VERITAS and multiwavelength observations of the Blazar B3 2247+381 in response to an IceCube neutrino alert

Abstract

While the sources of the diffuse astrophysical neutrino flux detected by the IceCube Neutrino Observatory are still largely unknown, one of the promising methods used towards understanding this is investigating the potential temporal and spatial correlations between neutrino alerts and the electromagnetic radiation from blazars. We report on the multiwavelength target-of-opportunity observations of the blazar B3 2247+381, taken in response to an IceCube multiplet alert for a cluster of muon neutrino events compatible with the source location between May 20, 2022 and November 10, 2022. B3 2247+381 was not detected with VERITAS during this time period. The source was found to be in a low-flux state in the optical, ultraviolet and gamma-ray bands for the time interval corresponding to the neutrino event, but was detected in the hard X-ray band with NuSTAR during this period. We find the multiwavelength spectral energy distribution is well described using a simple one-zone leptonic synchrotron self-Compton radiation model. Moreover, assuming the neutrinos originate from hadronic processes within the jet, the neutrino flux would be accompanied by a photon flux from the cascade emission, and the integrated photon flux required in such a case would significantly exceed the total multiwavelength fluxes and the VERITAS upper limits presented here. The lack of flaring activity observed with VERITAS, combined with the low multiwavelength flux levels, and given the significance of the neutrino excess is at 3$\sigma$ level (uncorrected for trials), makes B3 2247+381 an unlikely source of the IceCube multiplet. We conclude that the neutrino excess is likely a background fluctuation.

I Introduction

The IceCube Neutrino Observatory [3] has opened up a new window to extreme environments in the Universe, where cosmic rays are accelerated to ultra-high energies. Any detection of a TeV – PeV neutrino emitting source would help to directly answer some century-long open questions on the origin of cosmic rays. The neutrinos observed must be produced in cosmic-ray interactions [62]. Moreover, while both gamma rays and neutrinos point directly back to their sources, the latter are also not substantially attenuated as they travel through space. Higher energy gamma rays readily interact with background photons leading to pair-production and attenuation. The IceCube Observatory has detected a diffuse flux of astrophysical neutrinos [1], with the observed isotropic distribution of arrival directions suggesting an extragalactic origin. A number of extragalactic sources are proposed as candidate high-energy neutrino emitters. These include blazars, a subclass of active galactic nuclei (AGN) with relativistic jets pointing along our line-of-sight (for a review see [63]).

The IceCube Collaboration has reported evidence for neutrino emission from GeV gamma-ray sources, TXS 0506+056 [78] and NGC 1068 [49]. While the gamma-ray flux from the active galaxy NGC 1068 is attenuated and not variable, the TeV detected blazar TXS 0506+056 exhibits gamma-ray flaring episodes [11]. This, in addition to, the difference in AGN type of the sources and their respecitive redshifts, along the two sources representing different neutrino spectra [49], suggests the existence of at least two populations of extragalactic neutrino sources. The 2017 gamma-ray flare of TXS 0506+056 in spatial and temporal coincidence with the $\sim$ 290 TeV neutrino event IceCube-170922A [55, 78] also highlights the importance of flaring blazars in the search for neutrino sources.

More recently, IceCube detected a track-like event with an energy 171 TeV on December 8, 2021, which was found to be near the blazar PKS 0735+178 [47]. An association between the IceCube event and the PKS 0735+178 flare was investigated in multiwavelength studies [72, 14], but no firm connection was established. The above alerts are examples of single high energy ($>$ 100 TeV) neutrino events; other such studies can be found in [52], [36], [42], [58] and [5]. These alerts have been publicly distributed by IceCube since 2016. The information about these events is shared using the General Coordinates Network (GCN), an open-source platform created by NASA to receive and transmit low-latency alerts about astronomical transient phenomena. These single-event alerts have a typical localization uncertainty of $\sim 1^{\circ}$ and the region of interest (RoI) defined by the neutrino localization uncertainty often contains potential neutrino sources such as AGN or transient sources [4, 6].

While higher energy neutrinos are believed to have a higher probability of being of astrophysical origin, it should also be noted that the observation of multiple neutrinos from a single location is also a signature for neutrinos of astrophysical origin. In light of this, the IceCube Collaboration has been operating a complimentary approach known as the Gamma-ray Follow-Up (GFU) program [48, 31] since 2012. The goal of the GFU program is to enable prompt follow-up investigations of known gamma-ray sources for which IceCube has detected a cluster of candidate neutrino events, typically having energies $\sim$ 1 TeV, above a pre-defined significance, using imaging atmospheric Cherenkov telescopes (IACTs).

The GFU program utilizes a realtime neutrino event selection that selects events that are consistent with muon tracks arising from muon neutrino interactions in or near the IceCube detector volume. After significant data reduction, this sample is dominated by events arising from the interaction of cosmic rays in the atmosphere, with atmospheric neutrinos – an irreducible background - dominating the sample from the northern sky, and atmospheric muons dominating in the southern sky. The GFU searches time windows that span from seconds to as long as 180 days, looking for significant clustering in time and space of neutrinos in this sample. The alert thresholds are set so that random background fluctuations are suppressed to low levels, requiring a pre-trials statistical significance of $3\sigma$, and generate alerts at a frequency of about 10 per year.

The IceCube GFU program selected a list of 190 sources to be monitored for neutrino clusters and observable with VERITAS. These objects include all extragalactic TeV sources detected with IACTs, the Galactic Center, and the Crab Nebula, as well as sources from the 3FGL [12] or 3FHL [17] catalogs based upon the following criteria [31]:

• Extragalactic sources having a known redshift, z $\leq$ 1.0. This is done because sources at higher redshifts are difficult to detect with current generation IACTs due to absorbtion by extragalactic background light (EBL).

• 3FGL sources with variability index $>$ 77.2 and 3FHL sources for which the number of Bayesian blocks from variability analysis $>$ 1. This is done because the majority of detections with current generation IACTs occur during flares.

• Sources with a maximum elevation of $>45^{\circ}$ at the VERITAS site, in order to allow for optimal observations with VERITAS.

• Assuming the source produces gamma-ray flares with a 10-fold increase on the average Fermi-LAT flux, the extrapolated flux above 100 GeV exceeds the VERITAS $5\sigma$ sensitivity within 5 hours of observations, allowing for detection in short observing times.

This study investigates a GFU alert report privately shared with VERITAS by IceCube. It contained information on alerts of a cluster of muon neutrino candidate events from directions compatible with the source B3 2247+381, initially comprising 4 alerts received at a significance of $>3\sigma$ between August 10, 2022 and September 22, 2022. These initial alerts triggered the VERITAS target-of-opportunity (ToO) observations and initiated a multiwavelength campaign including NuSTAR observations. In total, there were 7 alerts over a duration of 174 days between May 20, 2022 and November 10, 2022, and are shown in Fig. 1. It should be noted that these 7 alerts were likely not independent and rather one primary alert developing in realtime, and retriggering multiple times as more events came in. The significance of these alerts, taking into account the event directions, angular uncertainties, event energies, and the time duration of the time window was found to be $3.2\sigma$. The corresponding false alert rate, quantifying how often the observed significance (or higher) is found at this location in a background-only scenario, is 0.0355 per year.

Figure 1: The evolution of alert significance at the location of B3 2247+381. In total, 7 alerts were received over a duration of 174 days between May 20, 2022 and November 10, 2022 and the significance of these alerts, taking into account their directions, angular uncertainties and energies, as well as the duration of the time window, was found to be $3.2\sigma$. It should be noted, this significance value is not fully corrected for trials as the GFU search includes trials corrections for the multiple time windows searched per source but not for the number of catalog entries searched.

The initial event on May 20, 2022 (MJD 59719) corresponded to a muon neutrino candidate event with an energy proxy of $E_{\mu}=~80.6$ TeV. The energy proxy is determined from the energy loss of the muon track and is described in Section 9.1 of [2]. This initial event pushed the significance of the neutrino cluster above 2.5$\sigma$, which then reached the alert threshold of 3$\sigma$ around August 9, 2022 (MJD 59800), resulting in a GFU alert to be triggered. Subsequent neutrinos with much lower energy proxies, $0.5$ TeV $<E_{\mu}<6.1$ TeV, continued to push the significance of the neutrino multiplets from this location to above 3$\sigma$, triggering six more GFU alerts, even though the cumulative significance of the neutrino multiplet did not continue to grow, as seen in Fig. 1.

The privately distributed alert stream was formed by clusters of neutrino events in time and space around B3 2247+381. It should be noted that the $3.2\sigma$ significance stated in the IceCube report is not fully corrected for trials as the GFU search includes trials corrections for the multiple time windows searched per source but not for the number of catalog entries searched. Furthermore, the GFU sample is dominated at this location by neutrino events with energies between 0.5 TeV and 6 TeV, hence these are highly likely to be atmospheric neutrino background. As the associated source of GFU alerts is already known, the goal of IACT and multiwavelength observations is to investigate possible changes to the state of the source, for example spectral changes or flaring states.

The BL Lacertae object B3 2247+381 ($z~=~0.119$; [32]) was first detected in the very high energy (VHE) regime with the MAGIC telescopes, based on 14.2 hours of good quality gamma-ray data collected between September 30, 2010 and October 30, 2010 [18]. These observations were triggered by a high optical state and yielded a significance level of 5.6$\sigma$. The observations revealed a relatively soft VHE spectrum with a photon index of $-3.2\pm 0.6$. Furthermore, no significant short-term flux variability was observed and the spectral energy distribution (SED) was successfully modelled using a one-zone synchrotron self-Compton (SSC) model. Moreover, the observed flux was found to be consistent with an upper limit obtained with MAGIC during a prior observation taken in 2006 during a low optical state, and no connection between the high optical state and the VHE gamma-ray emission could be established. The inclusion of B3 2247+381 in the IceCube GFU list was motivated by this initial detection at VHE energies.

The aim of this paper is to investigate the VERITAS and multiwavelength follow-up to the GFU alert. We study B3 2247+381 in other wavebands to characterize the transition between the low-energy and high-energy components of the broadband SED. Furthermore, we aim to present a proof-of-concept study of a follow-up to an alert within the framework of the GFU program and also investigate if there were any significant changes to the electromagnetic emission state of the source during the time interval corresponding to the alert.

II Observations and Data Analysis

II.1 VERITAS

The Very Energetic Radiation Imaging Telescope Array System (VERITAS) is an array comprising four 12 m IACTs. It is located at the Fred Lawrence Whipple Observatory (FLWO) in southern Arizona, USA ($30^{\circ}$ 40’ N, $110^{\circ}$ 57’ W, 1.3 km above sea level; Holder 46). Each telescope contains a camera comprising 499 photomultiplier tubes, and covering a field of view diameter of 3.5^∘. When the highest-energy photons enter the atmosphere, they produce an air shower containing particles which in turn produce a burst of blue light known as Cherenkov radiation. Gamma-ray observatories like VERITAS effectively use the whole atmosphere as their detector, and track the blue visible light produced from the air shower using optical telescopes.

VERITAS is capable of detecting gamma rays having energies in the range from 85 GeV to above 30 TeV, with an energy resolution of $\Delta E/E~\sim$ 15% (at 1 TeV) and an angular resolution of $\sim$ 0.1^∘ (68% containment at 1 TeV). In its current configuration, VERITAS typically detects a source with a flux of 1% of the steady state gamma-ray flux of the Crab Nebula at a statistical significance of $>5\sigma$ after 25 hours of observation. However, it should be noted that the time required to reach a $5\sigma$ detection is longer for sources with a spectrum softer than the Crab Nebula (photon index $\Gamma>2.49$) or those observed at large zenith angles ($\theta_{z}>40^{\circ}$) [66].

Following the IceCube alert, VERITAS started ToO observations of B3 2247+381 on September 22, 2022 (MJD 59844). VERITAS observed the source for 5 hours between September 23, 2022 and October 3, 2022 (MJD 59845 - MJD 59855) at a mean elevation of $73^{\circ}$ while hindered by poor weather on several nights during this period. These observations were performed using a standard “wobble” observing mode [34] with a $0.5^{\circ}$ offset in each of the four cardinal directions in order to simultaneously determine the background event rate. Furthermore, quality cuts were applied to the data set to remove events affected by bad weather. The VERITAS data were analyzed using the Eventdisplay analysis package [59], and independently confirmed with the VEGAS analysis package [27], yielding consistent results.

The VERITAS analysis parameterizes the principal moments of the elliptical shower images before applying a set of cuts to these parameters in order to reject cosmic-ray background events. The cuts are determined from a boosted decision tree algorithm [56], optimized for soft-spectrum sources ($\Gamma\sim 4$), and previously trained on gamma-ray shower simulations. During this process, we rejected events having fewer than two telescope images. Gamma-ray candidate events which fall within a squared angular distance, $\theta^{2}~\leq$ 0.008 deg², between the reconstructed event origin and B3 2247+381 are considered in the source (ON) region. Furthermore, the background is estimated using the reflected region model [24], where circular background OFF regions having the same size as the ON region are placed at the same radial distance from the center of the camera.

B3 2247+381 was not detected with VERITAS during the time interval investigated in this study. An excess of only 11 gamma-ray candidate events was recorded in the source region with $N_{\text{ON}}$ = 141 ON events, $N_{\text{OFF}}$ = 2384 OFF events, a background normalization factor, $\alpha=0.055$, and $N_{\text{Excess}}=N_{\text{ON}}-\alpha N_{\text{OFF}}=10.5$, corresponding to a statistical significance of 0.88$\sigma$, calculated following the method of [57]. The upper limit at 99% confidence level for the average integral flux above 200 GeV is $3.6\times 10^{-12}$ cm^-2 s^-1, or 1.5% of the Crab Nebula flux above the same energy threshold. All VERITAS upper limits were computed assuming a power-law spectrum with a photon index of 3 [69, following].

Furthermore, the analysis of archival observations of B3 2247+281 taken with VERITAS and comprising a total of 27 hours of quality-selected data accumulated from September 2009 to November 2020 found an excess of 65 gamma-ray candidate events in the source region with $N_{\text{ON}}$ = 836 ON events, $N_{\text{OFF}}$ = 14130 OFF events, a background normalization factor, $\alpha=0.055$, and corresponding to a statistical significance of 2.2$\sigma$ and a flux upper limit of $1.7\times 10^{-12}$ cm^-2 s^-1 above a threshold of 200 GeV. This can be compared to the time-averaged flux limit during the IceCube alert period of $3.6\times 10^{-12}$ cm^-2 s^-1. The differences between these limits likely originate from the much greater exposure time of the archival data (301 minutes of data during the IceCube alert period compared to 1683 minutes of archival data), which leads to a more constraining limit, rather than any intrinsic gamma-ray variability of the source. For comparison, the integral flux of B3 2247+381 above 200 GeV during the 2010 MAGIC detection was (5.0 $\pm$ $0.6_{\text{stat}}$ $\pm$ $1.1_{\text{sys}}$) $\times~10^{-12}$ cm^-2 s^-1 [18]. A lightcurve of the VERITAS flux above a threshold of 200 GeV and binned in monthly intervals is shown in Fig. 2. The VERITAS energy spectral limits, shown in the multiwavelength SED in Fig. 3, were derived using a spectral binning of 5 energy bins per decade in energy.

II.2 Fermi-LAT

The Fermi-Large Area Telescope (LAT; Atwood et al. 22) is a pair conversion telescope capable of detecting gamma-ray photons in the energy range from 20 MeV to above 500 GeV. The pair conversion process forms the basis for the operation of the Fermi-LAT by providing a unique signature for gamma rays. This distinguishes them from charged cosmic rays and allows a determination of the incident photon directions via the reconstruction of the trajectories of the electron positron pairs. Primarily operating in survey mode, the Fermi-LAT scans the entire sky every three hours. In this paper, we initially analyzed Fermi-LAT data during the IceCube neutrino alert period, between MJD 59719 and MJD 59893, corresponding to midnight on May 20, 2022 until midnight on November 10, 2022. Throughout the analysis, we use the Fermi Science Tools version 2.20, FERMIPY version 1.2 [83], in conjunction with the latest PASS 8 instrument response functions (IRFs; [21]).

The Fermi-LAT data were processed using a binned maximum likelihood analysis. Photons with energies between 100 MeV and 300 GeV detected within a RoI of radius 10^∘ centered on the location of B3 2247+381 were selected for the analysis. We selected only photon events from within a maximum zenith angle of 90^∘ in order to reduce contamination from background photons from the Earth’s limb, produced from the cosmic-ray interactions with the upper atmosphere.

A spatial bin size of 0.1^∘ per pixel and two energy bins per decade were used. All sources contained in the 4FGL-DR3 catalog [8] within 20^∘ of the RoI center, were included in the model with their spectral parameters fixed to their catalog values. This takes into account the gamma-ray emission from sources lying outside the RoI which might contribute photons to the data, especially at low energies, due to the size of the point spread function of the Fermi-LAT. Since the time interval considered in this work is beyond that covered in the 4FGL-DR3 catalog, the gtfindsrc routine was also applied to search for any additional point sources present in the data and not included in the catalog. No significant additional point sources, having a test statistic (TS; [60]) $\geq$ 9 (roughly corresponding to a significance of $\sim$ 3$\sigma$) were detected, indicating that all sources in the data had been accounted for. Moreover, the contributions from the isotropic and Galactic diffuse backgrounds were modeled using the most recent templates for isotropic and Galactic diffuse emission, iso_P8R3_SOURCE_V3_v1.txt and gll_iem_v07.fits, respectively.

The normalization factor for both the isotropic and Galactic diffuse emission templates were left free along with the spectral normalization of all modeled sources within the RoI. Moreover, the spectral shape parameters of all modeled sources within 3^∘ of B3 2247+381 were left free to vary while those of the remaining sources were fixed to the values reported in the 4FGL-DR3 catalog. B3 2247+381 was found to have an integral flux upper limit of $1.43\times 10^{-9}$ $\text{cm}^{-2}\text{s}^{-1}$ in the energy range 100 MeV – 300 GeV during the IceCube neutrino alert period.

We then repeated the same analysis routine over a wider time interval August 4, 2008, the start of the Fermi-LAT mission until midnight on June 1, 2023 (MJD 54683 – MJD 60096). B3 2247+381 was detected at a statistical significance of 23$\sigma$ (TS = 542.3), corresponding to an integral flux of $(2.84~\pm~0.42)~\times 10^{-9}$ $\text{cm}^{-2}\text{s}^{-1}$. The spectrum was found to be best modeled by a power law (PL):

$$\centering\frac{dN}{dE}=N_{0}\left(\frac{E}{E_{0}}\right)^{-\Gamma}\@add@centering$$ (1)

where $N_{0}$ is the normalization, $E_{0}$ is the pivot energy, and $\Gamma$ is the spectral index. The corresponding best-fit spectral parameters are $N_{0}=(2.25\pm 0.17)\times 10^{-14}$ $\text{cm}^{-2}\text{s}^{-1}\text{MeV}^{-1}$, and $\Gamma~=~1.71\pm 0.05$. These are, within uncertainties, compatible with the spectral values listed in the 4FGL-DR3 catalog [8] for this source, explicitly $N_{0}=(2.60\pm 0.19)\times 10^{-14}$ $\text{cm}^{-2}\text{s}^{-1}\text{MeV}^{-1}$ and $\Gamma=1.74\pm 0.05$.

In order to investigate the temporal variability of the gamma-ray flux of the source during the full time interval, the lightcurve of the integral flux in the energy range 100 MeV–300 GeV, shown in Fig. 2, is calculated with 90-day bins and keeping the spectral index fixed to the catalog value. We find no evidence of variability for B3 2247+381 in the long-term Fermi-LAT lightcurve (see Section III) or in the Fermi-LAT differential energy spectrum, shown in the multiwavelength SED in Fig. 3, which considers the wider time interval, from MJD 54683 to MJD 60096, for improved statistics.

Figure 2: Multiwavelength lightcurves of B3 2247+381 between August 4, 2008 until midnight on June 1, 2023 (MJD 54683 – MJD 60096). The blue narrow band corresponds to the IceCube neutrino excess period and the red narrow band corresponds to the time interval of the MAGIC detection [18]. (A) Monthly-binned VHE lightcurve for VERITAS observations above an energy threshold of 200 GeV. Also shown for comparison are MAGIC data from the 2010 detection [18]. (B) 90-day-binned Fermi-LAT lightcurve in the 0.1–300 GeV band. Upper limits are shown for individual bins having a significance of lower than $2\sigma$. (C) Swift-XRT and NuSTAR fluxes in the 0.3–10 keV band and 3–12.5 keV band, respectively. To aid visual comparison, the NuSTAR data have been scaled up by a factor of 10. (D) Optical lightcurve in all filters using data from the 48” optical telescope at the FLWO. Errors plotted are statistical uncertainties only. (E) Optical flux observations: Optical data from the ASAS-SN Sky Patrol in the g and V band (green and magenta, respectively) and ATLAS R band data (in red).

II.3 NuSTAR

The Nuclear Spectroscopic Telescope Array (NuSTAR) is a space-based telescope capable of detecting hard X-ray photons in the energy range between 3–79 keV. NuSTAR has an angular resolution of 18” (FWHM) and two focal plane modules (FPMA and FPMB) with a $13^{\prime}\times 13^{\prime}$ field of view [43]. The NuSTAR ToO observations of B3 2247+381 (ObsID = 80802636002, exposure = 40.5 ks) took place on September 26, 2022 (MJD 59848), following the IceCube alert received via MoU. The NuSTAR data were processed using NuSTAR data analysis software (NuSTARDAS) version 2.1.1 contained within HEASOFT version 6.29.

The source spectra were extracted from a circular region of radius $1^{\prime}$, chosen to maximize the source significance with respect to the background. The spectra were binned such that there was a 5$\sigma$ significance within each bin. The binned spectra, shown in the multiwavelength SED in Fig. 3, from FPMA and FPMB were simultaneously fitted using XSPEC [20] in the energy range between 3 keV and 12 keV, as beyond this the background begins to dominate. B3 2247+381 was strongly detected in the hard X-ray band with NuSTAR at a significance of $42\sigma$ during these observations. The NuSTAR spectra were well fit by an absorbed power law with the Galactic column density $N_{H}$ in the direction of B3 2247+381 held fixed at 1.57 $\times 10^{21}$ atoms cm^-2 [82]. The best-fit photon index of the observations was found to be $\Gamma=2.89\pm 0.12$ and a $\chi^{2}_{red}=70.57/63=1.12$ was obtained. The time-averaged integral flux of B3 2247+381 measured with NuSTAR during the period investigated is $(1.45\pm 0.05)\times 10^{-12}$  erg  cm^-2  s^-1.

II.4 Swift

The X-Ray Telescope (XRT), onboard the Neil Gehrels Swift Observatory, is a grazing-incidence focusing X-ray telescope sensitive to energies from 0.2 keV to 10 keV [41, 25]. Swift-XRT observations of B3 2247+381 were taken for the time interval between May 20, 2022 until midnight on November 10, 2022 (MJD 59719 – MJD 59893). A total observing time of $\sim$5.5 ks was accumulated during this period in photon-counting (PC) mode. The Swift-XRT lightcurve in the energy range 0.3–10 keV was retrieved from the public online tool “the Swift-XRT data products generator”  [29, 30] and is shown in Fig. 2.

The Ultraviolet/Optical Telescope (UVOT), also onboard the Neil Gehrels Swift Observatory, is a photon-counting telescope. It is sensitive to photons having energies ranging roughly between 1.9 eV and 7.3 eV [70]. Swift-UVOT observations are performed in parallel to the Swift-XRT measurements in six filters with central wavelengths of $V$ (5468 Å), $B$ (4392 Å), $U$ (3465 Å), $UVW1$ (2600 Å), $UVM2$ (2246 Å), and $UVW2$ (1928 Å). The Swift-UVOT data were analyzed with HEASOFT version 6.29. A source region with a radius of $5.0^{\prime\prime}$ centered on B3 2247+381, and a background region of the same size, away from the blazar and containing no known sources in any band, were used to extract signal and background counts, respectively. The magnitude of the source was then computed using uvotsource and converted to flux using the zero-point for each of the UVOT filters from Poole et al. [68].

The measured optical-UV fluxes, obtained from Swift-UVOT observations between September 23, 2022 and November 10, 2022 (MJD 59845 - MJD 59893) are shown in the multiwavelength SED in Fig. 3. This includes correcting for interstellar extinction using the approach of [71], assuming a reddening of $E(B-V)=0.0167$ [73]. To build a consistent SED, the contribution of the host-galaxy emission has been removed in the UVOT filters, using the R filter and the HIRAC-camera observations of TeV J2250+384 (B3 2247+381) with the Nordic Optical Telescope (NOT) [64]. Assuming the simplest model of [28], in conjunction with an elliptical-galaxy template from PEGASE.2 [33], the host estimation in the UVOT bands within a 5” aperture was estimated at V: $2.60\times 10^{-12}$, B: $1.63\times 10^{-12}$, U: $4.92\times 10^{-13}$, UVW1: $1.82\times 10^{-13}$, UVM2: $1.57\times 10^{-13}$, and UVW2: $1.74\times 10^{-13}$ erg cm^-2 s^-1, respectively.

II.5 Optical

The All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 75, Kochanek et al. 54) is an automated program aimed at routinely surveying the entire visible sky for bright transient sources with minimal observational bias. It can reach a depth of roughly 17 mag in V band and 18.5 mag in g band. While ASAS-SN originally comprised two stations, located at the Cerro Tololo International Observatory (CTIO, Chile) and the Haleakala Observatory (Hawaii), respectively, in 2017 three further stations were added, including a second unit at CTIO and one unit each at the McDonald Observatory (Texas) and the South African Astrophysical Observatory. ASAS-SN observations are available in two optical bands: the V band, centered at $\sim$551 nm and corresponding to observations from the two original stations, and the g band, centered at $\sim$480 nm and corresponding to observations performed with the three new stations. For the purpose of this study, we obtain B3 2247+381 observations from the ASAS-SN Sky Patrol for the time interval between August 4, 2008 until midnight on June 1, 2023 (MJD 54683 – MJD 60096), and shown in Fig. 2.

The Asteroid Terrestrial-impact Last Alert System (ATLAS; [79]) is a high cadence all-sky survey system comprising four independent units, one each at Haleakala and Mauna Loa in the Hawaiian islands in the Northern Hemisphere and one each at the El Sauce Observatory, Chile and the South African Astronomical Observatory in the Southern Hemisphere. ATLAS is optimized to be efficient for finding potentially dangerous asteroids, as well as in tracking and searching for highly variable and transient sources, such as AGN. Optical observations of B3 2247+381 with ATLAS in the R band, centered at $\sim$679 nm and having a typical cadence of one data point per two days, are shown in Fig. 2.

Images of B3 2247+381 have also been collected with the 48” optical telescope at the FLWO as part of an ongoing study of more than seventy blazars that aims to track and analyze their behavior over time. These data were collected using SDSS r^′, SDSS i^′, Harris V, and Harris B filters. Images are processed by an automated pipeline that uses aperture photometry to determine magnitudes. For each image taken at FLWO, a measurement is made of the blazar’s magnitude as well as the magnitude of its check star, which is used to establish and maintain the quality of photometric data. Additional quality is enforced by removing outliers. A total of 517 images were collected over 135 nights, between May 2015 and July 2023. An optical lightcurve of B3 2247+381, showing the resulting 453 good quality measurements in the different bands with the 48” optical telescope, is presented in Fig. 2. The data have not been corrected for optical extinction or any potential host galaxy contribution.

III Discussion

The SEDs of blazars typically display two broad emission features. While the first feature in the radio to X-ray waveband is commonly attributed to synchrotron emission from relativistic electrons and positrons within the jet, the origin of the second feature in the X-ray to gamma-ray band is less clear. In leptonic models (for example Jones et al. 51, Sikora et al. 76), this emission is attributed to the IC scattering between the energetic leptons in the jet and a field of lower energy photons. These may be the same synchrotron photons produced by the leptons (synchrotron self-Compton [SSC] model) in the case of BL Lacertae objects like B3 2247+381, or photons from an external source (external inverse Compton [EIC] model) in the case of Flat Spectrum Radio Quasars. The sources of external photons include the accretion disk, the broad-line region, and the molecular torus.

The connection between the parent cosmic-ray proton energy and the energies of the neutrinos and photons produced in their interactions is not straightforward nor fully understood. While it is expected that a neutrino would carry approximately $5\%$ of the parent proton energy, Doppler boosting in the jet can lead to the neutrino gaining back energy in the observer frame. This can result in a neutrino with an energy of  100 TeV in the observer frame being produced by protons with maximum energies in the jet frame as low as  100 TeV (depending on the Doppler boost factor). It is anyway possible for the maximum proton energy to be much higher, with some models (for example Cerruti et al. 26) deriving maximum proton energies above $10^{17}$ eV.

As the medium is expected to be opaque to the TeV-PeV pionic photons being produced, the electromagnetic cascading of these gamma rays results in some, or all, of their energy being redistributed into lower energy observable bands (keV to MeV), posing an additional challenge to deriving proton energy information from electromagnetic and neutrino observations. Additionally, disentangling the hadronic contribution to the SED is particularly challenging as it is likely subdominant when compared to the main leptonic component (for example [53]), unless the proton energies are well above the PeV range (for example Cerruti et al. 26). The modeling of these different contributions and their connections is an active area of research in the blazar modeling community.

[38] demonstrate that a moderate enhancement in the number of cosmic rays during a flare leads to a proportional rise in the neutrino flux, which is accompanied by hard X-rays and TeV gamma rays. Therefore, a multi-messenger approach combining both neutrino and electromagnetic observations is vital for neutrino source identification. Moreover, hard X-ray observations from NuSTAR, as well as gamma-ray observations at TeV energies have provided some of the strongest constraints on hadronic emission models in blazars [14]. In cases where gamma-rays and neutrinos cannot be linked with a simple model, or VHE radiation is not observed, other models can be used that associate neutrinos with hard X-rays cascading from gamma-ray interactions [7]. Such “hybrid models” combine a leptonic explanation of the observed emission with a hadronic component, the strength of which is dictated by the X-ray observations [38]. Finally, in addition to SEDs, [38] also highlight the importance of hard X-ray and gamma-ray variability since the strongest temporal correlations are expected between these two bands with neutrinos.

We characterize the gamma-ray variability of B3 2247+381 by evaluating the goodness-of-fit parameter for a constant value function fit to the lightcurve. As seen in panel (B) of Fig. 2, the Fermi-LAT lightcurve is very well fit by a constant function at a flux value of  ($3.2\pm 0.3)~\times 10^{-9}$ $\text{cm}^{-2}\text{s}^{-1}$, with a $\chi^{2}/\text{n.d.f}=40.8/52=0.8$, which suggests the absence of any GeV variability in the source for the time interval investigated. Furthermore, a $\chi^{2}/\text{n.d.f}=152/52=2.9$ was obtained when testing for variability using the $\text{ts}_{\text{var}}$ routine in FERMIPY, and also accounts for the upper limits in the lightcurve. The alternate hypothesis that a source’s lightcurve is not constant, $\text{ts}_{\text{var}}$, as defined in equation 4 of [65] is distributed as as a $\chi^{2}$ distribution with 23 degrees of freedom, where $\text{ts}_{\text{var}}$ $>41.6$ corresponds to a variable source with 99% confidence. For comparison, no variability was observed for the source during the MAGIC detection, where a fit of a constant function to the Fermi-LAT flux points yielded a flux value of  $(3.7\pm 0.5)~\times 10^{-9}$ $\text{cm}^{-2}\text{s}^{-1}$, with a $\chi^{2}/\text{n.d.f}=0.7$ [18].

There is evidence of variability on monthly timescales in the long-term X-ray data, as seen in the increased X-ray flux measured by Swift-XRT in Figure 2. This lightcurve is poorly fit by a constant function, corresponding to a $\chi^{2}/\text{n.d.f}=42.2$. It should be noted that there were no Swift-XRT observations of B3 2247+381 for the time period between the 2010 MAGIC detection and the IceCube alert, making it difficult for a detailed study of the long-term variability. Nevertheless, there was an increase in the X-ray flux state of B3 2247+381 corresponding to the MAGIC detection, with an average photon count rate of $0.58~\pm~0.11$ cts/s, compared to $0.23~\pm~0.03$ cts/s for the period corresponding to the neutrino alert. For comparison, no strong X-ray intra-night variability was found by [18] during the MAGIC detection.

The optical data of B3 2247+381 show stronger evidence of variability, with the ASAS-SN g observations in particular showing an increase in flux during the time interval corresponding to the neutrino alert as seen in Fig. 2. All three optical lightcurves are fit by a constant function, with $\chi^{2}/\text{n.d.f}$ of 2.04, 1.94 and 1.61 for the ATLAS, ASAS-SN g and ASAS-SN V data respectively. The time-averaged optical flux measured with ATLAS in the R band for the entire observation period was $0.95\pm 0.05$ mJy. For comparison, B3 2247+381 was observed to be in a steady state of $\sim 1.8$ mJy between 2006–2009 in the R band by the Tuorla group, using the 35 cm telescope at the KVA observatory on La Palma, Canary Islands, Spain before the optical flare in late summer in 2010, when it reached an average flux of 2.4 mJy [18]. The FLWO B,V,r^′ and i^′ data in Fig. 2 also support the conclusion that B3 2247+381 was not in a particularly bright optical state at the time of the neutrino alert.

It should be noted that we cannot conclusively claim nor rule out any associations between the neutrino multiplet and the blazar, based on these observations alone. In the case of previous IceCube neutrino alerts with possible blazar counterparts, TXS 0506+056 [11] and PKS 0735+178 [14], flaring states were detected through multiwavelength observations that were temporally coincident with the neutrino alerts, strengthening the statistical significance of the correlation between the neutrino and electromagnetic emission. Nevertheless, we still aim to shed light on the mechanisms responsible for the emission from B3 2247+381 during the time interval investigated by modeling the multiwavelength SED. The multiwavelength broadband SED of B3 2247+381, shown in Fig. 3, comprises the upper limits obtained with VERITAS, along with observations with NuSTAR and optical telescopes. Although these observations were not conducted exactly simultaneously, the lack of observed variability on nightly timescales in any of these wavelengths allows for the combination of these observations into a broadband SED. The Fermi-LAT SED points considered are time-averaged over a wider interval due to low statistics and no apparent variability in the gamma-ray flux from MJD 54683 to MJD 60096. The synchrotron component is loosely constrained by the relatively flat optical/UV spectra, while the hard X-ray spectra observed with NuSTAR is found to occur in the low-energy component of the SED, as opposed to characterizing the transition from low- to high-energy, as seen in the blazar PKS 0735+178 [14], with the transition occurring well above 10 keV.

In order to interpret the multiwavelength behavior of B3 2247+381, we model the broadband observations using Bjet_MCMC [44], a new public tool that can simply fit the broadband SEDs of blazars. Bjet_MCMC has previously been used to model several jetted AGN emitting in the VHE band, such as AP Librae [45], HESS J1943+213 [19], 1ES 1215+303 [81], and also PKS 1222+216 and TON 599 [15]. We consider a simple one-zone pure SSC radiation model consisting of a single blob with a broken power law differential electron energy density (EED), a scenario most commonly used to model BL Lacertae objects such as B3 2247+381.We expect a SSC contribution to blazar SEDs, regardless of whether or not the IceCube alert originates from astrophysical neutrinos associated with B3 2247+381. If the neutrino emission indeed originates from B3 2247+381, we would expect an additional hadronic component to the SED, therefore a simple one-zone SSC model would not be able to fully explain the observed SED, as was seen in the case of [14]. The SED fit is performed using the Markov-chain Monte-Carlo (MCMC) method.

The pure SSC model contains by default nine free parameters, including three describing the blob, namely the Doppler factor, $\delta$, the magnetic-field strength, $B$, and the blob radius, $R$, and six parameters which define the EED, namely the particle density factor, $N_{e}$, the two spectral indices of the spectrum, $n_{1}$ and $n_{2}$, the break Lorentz factor of the electron distribution, $\gamma_{\text{break}}$, and the minimum and maximum Lorentz factors of the EED, $\gamma_{\text{min}}$ and $\gamma_{\text{max}}$, respectively. The redshift is fixed $z~=~0.119$ [32] and a flat $\Lambda$CDM cosmology is assumed [23]. This choice of model also allows for comparison with the parameters obtained during the 2010 MAGIC detection [18]. We consider a wide range of parameters in order to make as few assumptions as possible, and many parameters are in log scale to allow the walkers to move through multiple orders of magnitude.

Figure 3: The broadband SED of B3 2247+381. The VERITAS spectral limits are time-averaged over the time interval of the ToO observations between September 23, 2022 and October 3, 2022 (MJD 59845 – 59855). The Fermi-LAT spectrum, time-averaged over a wider time interval between August 4, 2008 and June 1, 2023 (MJD 54683 – 60096), is shown using purple points while those corresponding to the IceCube neutrino alert window between May 20, 2022 and November 10, 2022 (MJD 59719 – 59893) is shown using pink points. The NuSTAR spectra were obtained on September 26, 2022 (MJD 59848). The Swift-UVOT and Swift-XRT observations correspond to the time interval between September 23, 2022 and November 10, 2022 (MJD 59845 – 59893). The grey points represent archival data, provided by the Space Science Data Center–ASI, taken with Wide-field Infrared Survey Explorer (WISE) and the Swift-XRT. It should be noted that none of these archival data are considered in the fitting, as we consider only contemporaneous datasets. The MAGIC data from the 2010 detection [18] are shown with gray square markers. A minimum significance of two standard deviations is required for each flux point and upper limits at 99% confidence level are quoted otherwise. The blue curve and shaded region represents the one-zone synchrotron and SSC radiation model fit and the corresponding $1\sigma$ confidence interval respectively, obtained using Bjet_MCMC, with the residuals for each SED point considered shown in the bottom panel, and the parameters shown in Table 1. The Bjet_MCMC shows the SED including the effect of the extragalactic background light (EBL) using the [35] model. The VERITAS upper limits have therefore not been corrected for EBL absorption.

Table 1: The fit parameters and the parameter range corresponding to the $1\sigma$ confidence level of a simple one-zone synchrotron and SSC radiation model fitted to the multiwavength broadband SED of B3 2247+381. Also tabulated is the redshift, $z~=~0.119$, which is a fixed value during model fitting [32]. The model is shown in Fig. 3. Also shown for comparison are the parameters obtained during the 2010 MAGIC detection [18].

Parameter	Description	Scale	Parameters	$1\sigma$ range	MAGIC
B [G]	magnetic field strength	$\text{log}_{10}$	-2.0	[-4, -0.7]	-1.2
$N_{e}[\text{cm}^{-3}]$	particle density factor	$\text{log}_{10}$	6.0	[3.5, 7.3]	3.3
$\gamma_{\text{min}}$	low-energy cutoff	$\text{log}_{10}$	0.8	[0, 2.4]	3.5
$\gamma_{\text{max}}$	high-energy cutoff	$\text{log}_{10}$	6.9	[5.1, 7.3]	5.8
$\gamma_{\text{break}}$	energy break	$\text{log}_{10}$	5.5	[4.9, 6.6]	4.9
$n_{1}$	first index	linear	2.75	[2.68, 2.79]	2.0
$n_{2}$	second index	linear	5.65	[4.86, 6.55]	4.4
$\delta$	Doppler factor	linear	73	[38, 98]	35
$R[\text{cm}]$	blob radius	linear	$9.4\times 10^{15}$	[$4.4\times 10^{14}$, $1.9\times 10^{18}$]	$8\times 10^{15}$
$z$	redshift	linear	0.119	–	0.119

For this model, we considered 100 walkers, 5000 steps, and a burn-in phase of 1000 steps and the MCMC walkers display a good general $\chi^{2}$ convergence. Fig. 4 (see Appendix A) shows the Bjet_MCMC $\chi^{2}$ values of walkers at each step for the SED fit for all walkers, smallest $\chi^{2}$ and median $\chi^{2}$ respectively. Fig. 5 (see Appendix A) shows a corner plot of the posterior probability distribution of the free parameters from the SED fit. The values of the parameters used in the SED fit are reported in Table 1, and also shown for comparison are the parameters used to fit the SED during the MAGIC detection [18]. The fitted model returned a $\chi^{2}$ = 95.5 with 84 degrees of freedom, corresponding to $\chi^{2}_{\text{red}}=1.14$.

Most of the parameters are found to have good agreement with the values obtained in [18] and typical values seen for BL Lacertae objects [77]. The Bjet_MCMC SED fit investigates wide ranges for all parameters [44], and probes outside the parameter range of [18]. Moreover, the [18] fitting algorithm acts primarily on the electron normalization, source radius and Doppler factor, with only slight changes to the other parameters. However, B3 2247+381 is relatively weak in the Fermi-LAT energy range, and as a result no statistically significant flux variability could be observed in timescales of days to months. It should be noted that large uncertainties obtained for the parameters prevent us from drawing any strong conclusions. Finally, the TeV non-detection severely limits a detailed study of the IC component of the SED.

As there is no evidence of association between the neutrino flare and the blazar, we remain agnostic about the origin of the gamma-ray emission. The one-zone leptonic scenario seems to reasonably fit the SED, which would justify the absence of an association of an increase in multi-wavelength emission during the time of the neutrino alert. In this paper, we briefly explore the viability of the neutrino emission, assuming the alerts are associated with B3 2247+381. In a hadronic emission scenario, a significant number of accelerated protons or ions may exist in the jet and neutrino emission may occur through two different channels, of which p$\gamma$ interactions [37] are more likely. Neutrino emission through the $pp$ channel would require very large particle density in the jet and hence a large jet power [67].

As seen in Fig. 3, the VERITAS upper limits are not constraining and since there is no direct evidence for a cut-off in the GeV gamma-ray spectrum observed by Fermi-LAT, we can infer the maximum neutrino flux that provides a VHE gamma-ray non-detection in a hadronic scenario. Assuming no absorption of the gamma rays due to pair production in the broad and narrow line regions, the VERITAS upper limits can be used to make an approximation of the possibility of neutrino production, as the EBL absorption is relatively low at the redshift $z~=~0.119$ [32]. Using the simulated response of IceCube, with GFU event selection to a neutrino flux, the observed event excess (12.7 events with a best fit neutrino spectral index of -2.33 over a 174.0 day time period) would correspond a neutrino flux of

$$EF_{\nu}(E)\approx 5.8\times 10^{-11}\ \mathrm{erg\,cm^{-2}\,s^{-1}}.$$ (2)

Furthermore, the target photons that are required for pion production would be observed at $\sim$ 100 keV, above the energies covered by NuSTAR observations, but the flux in the SSC model is low, which mandates a high proton source power. Given the significance of the neutrino excess at 3$\sigma$ level (uncorrected for trials), and that B3 2247+381 is an unlikely source of the neutrino excess, we conclude that the neutrino excess is likely caused by a fluctuation in the irreducible atmospheric neutrino background. Furthermore, the cascade emission of charged and neutral pions in the 100 TeV energy band provides a gamma-ray flux constraint at a few hundred GeV, similar to the neutrino flux [39]. The upper limits obtained with VERITAS are a factor $\sim$ 20 lower and, assuming no gamma-ray absorption above 100 GeV, the expected rate of neutrinos would be around $\sim$ 0.05 per month and any association would likely be a chance coincidence.

IV Conclusions

In this paper, we report on the VERITAS and multiwavelength follow-up observations of B3 2247+381, triggered by a Gamma-ray Follow-Up alert received by private communication from IceCube. The report contained information on alerts for a cluster of muon neutrino candidate events from directions compatible with B3 2247+381 received at a significance of $>3\sigma$ between May 20, 2022 and November 10, 2022. While these observations did not yield a detection with VERITAS, B3 2247+381 was, for the first time, significantly detected in the hard X-ray band with NuSTAR.

The broadband SED of B3 2247+381 comprising the upper limits obtained with VERITAS, as well as extensive follow-up observations across the electromagnetic spectrum, in particular the time-averaged Fermi-LAT data, and contemporaneous Swift and NuSTAR observations were investigated. We considered a simple one-zone synchrotron and SSC radiation model consisting of a single blob accelerating electrons with a broken power law electron energy density. A good agreement was seen between the fitted model and the multiwavelength spectral data over the entire energy range considered, and this broadly agrees with the results presented in [18].

Furthermore, this investigation also serves as one of the first VERITAS follow-ups to an alert within the framework of the GFU program. It should be noted that while the primary objective of the GFU program is to notify VERITAS and other instruments in the collective observatory network of potential astrophysical neutrino flares from directions compatible with known gamma-ray emitters, the duration of the flares is not predetermined and can span from seconds to the order of months. In this study, we report no significant changes to the gamma-ray emission state of the source despite the significant detection in the hard X-ray regime for the period investigated. Having just a spatial correlation alone makes this source different from previously studied sources with VERITAS, namely TXS 0506+056 [9] and PKS 0735+178 [14], both of which coincided both spatially and temporally with a gamma-ray flaring state. In looking for spatial coincidence between the neutrino multiplet and gamma-ray sources, we also note that, unlike TXS 0506+056 and PKS 0735+178, which were independently associated with IceCube alerts, B3 2247+381 was a-priori selected for GFU monitoring, which includes already known or likely VHE gamma-ray emitters. The two situations are therefore different.

Moreover, assuming the neutrinos originate from hadronic processes within the jet, the neutrino flux would be accompanied by a photon flux from the cascade emission, and the integrated photon flux required in such a case would significantly exceed the total multiwavelength fluxes and the VERITAS upper limits presented here. Therefore, due to no evidence of flaring activity with VERITAS, combined with the low multiwavelength fluxes, B3 2247+381 is disfavored as the origin of the IceCube multiplet. The neutrino hotspot seen with IceCube may still be of astrophysical origin and further multiwavelength observations of B3 2247+381 and a better understanding of instrument uncertainties might allow for a potential study of multi-zone contributions along the jet and confirm whether or not hadronic contributions and variability can be entirely ruled out. More generally, future study of AGN potentially associated with neutrino alerts, in particular those also coincident with flaring periods that have improved statistics, are required to draw further conclusions. For example, the viability of neutrino emission and its possible implications on the multiwavelength SED are investigated for PKS 0735+178 in [14].

This study also highlights some of the challenges in searching for neutrino-emitting blazars, including the limited localization precision of the IceCube Observatory and the effect of poor weather conditions on IACT observations. Furthermore, the large number of gamma-ray blazars as potential counterparts makes the association between neutrino events and a gamma-ray blazar difficult. As in this study, the electromagnetic emission can often be explained by invoking leptonic models alone without any further need for an hadronic component. Moreover, the jet power and the proton luminosity parameters required in lepto-hadronic models are often too high, and exceed the Eddington limit, for a short period of activity.

These challenges can only be addressed with further multiwavelength follow-up investigations of flaring blazars in temporal and spatial coincidence with astrophysical neutrino alerts, both single high energy neutrino events and those within the framework of the GFU program. Furthermore, as seen in [74], such studies are even more effective when we combine observations from multiple IACTs in order to provide a more complete coverage of the entire sky. This also helps to account for cases where the visibility of a source from a single observatory site is adversely affected due to factors such as bad weather, the presence of the Sun or Moon, or technical problems.

Despite a non-detection at VHE wavelengths, this observation campaign will inform future IceCube GFU follow-up campaigns with VERITAS. We have raised our GFU alert threshold to trigger on neutrino events with a significance $>3.5\sigma$ and will continue to monitor the alerts for a further increase in significance. A rising neutrino alert significance would indicate that the alert location is continuing to emit high energy neutrinos. In contrast, in the case of the alert at the B3 2247+381 location, the initial rise in significance is due to a single event, which then continues to boost the cumulative significance of subsequent events of lower energy and significance, without further rise in significance. Continuing to observe GFU alerts with IACTs, such as VERITAS, provides important opportunities to explore the viability of blazars as neutrino sources. Since we expect hadronic emission mechanisms to produce VHE gamma rays along with IceCube-detected neutrinos, VHE SED points or upper limits are important for constraining whether or not a hadronic component is necessary to model the multiwavelength blazar emission, which would motivate a physical connection with the observed neutrinos. VERITAS receives 1-2 GFU events per observing season, which allows for sufficient observations to be conducted to detect the levels of flaring activity expected from the blazars that have been selected for GFU targets.

Furthermore, future neutrino detectors such as IceCube-Gen2 [50], KM3NeT (for example [16]), P-ONE [80], and others will strongly improve the all-sky sensitivity to high energy neutrino events, as their acceptances are complementary to IceCube’s. Finally, neutrino follow-up efforts remain a major focus of the wider multi-messenger community (for example [10, 40, 61, 13]) and will hopefully result in directly constraining the century-old puzzle of the origin of cosmic rays in the future.

Appendix A Additional Bjet_MCMC OUTPUTS FOR B3 2247+381

Fig. 4 shows the Bjet_MCMC $\chi^{2}$ values of walkers at each step for the SED fit for all walkers, smallest $\chi^{2}$ and median $\chi^{2}$ respectively. Fig. 5 shows a corner plot of the posterior probability distribution of the free parameters from the SED fit.

Figure 4: The Bjet_MCMC $\chi^{2}$ values of walkers at each step for the SED fit of B3 2247+381. Top left panel: all walkers, Top right panel: smallest $\chi^{2}$, Bottom panel: median $\chi^{2}$.

Figure 5: The corner plot of the posterior probability distribution of the free parameters from the SED fit of B3 2247+381.