© The Authors 2023

1. Introduction

According to observations to date, the energetics of diffuse gamma rays (Ackermann et al. 2015), neutrinos (IceCube Collaboration 2021), and cosmic rays (Fenu & Pierre Auger Collaboration 2017) share similar power-law energetics. At lower energies, though, there is more flux in neutrinos than could be expected from gamma rays, if they are produced by the same sources. Results suggest that there could be a population of neutrino sources being obscured in gamma rays (e.g., IceCube Collaboration 2021). The tension at lower energies (E < 100 TeV) between the diffuse gamma-ray background measured by the Fermi Large Area Telescope (Fermi-LAT) and the diffuse neutrino background measured by the IceCube Neutrino Observatory also suggests the existence of a population of gamma-ray dark neutrino sources (Murase et al. 2016; Murase 2019), for which rather an anti-correlation holds between the neutrino and gamma-ray fluxes, at least in the epochs of efficient neutrino emission.

After a significant 3σ connection was found between a ∼290 TeV IceCube neutrino alert (IC-170922A) and the gamma-ray flaring blazar TXS 0506+056, this blazar quickly came into the focus of multimessenger astronomy. Follow-up searches in the offline IceCube data revealed a neutrino flare from the direction of TXS 0506+056 at the turn of the year 2014–2015 with an even higher significance (3.5σ). The individual events of this neutrino flare contained energies one order of magnitude lower (∼10 TeV) than the 2017 event, and they do not satisfy the strict selection cuts for single high-energy alerts.

As it turned out, the 2014–2015 neutrino flare (IceCube Collaboration 2018) happened during a low gamma-state of the source candidate blazar TXS 0506+056 (Garrappa et al. 2019a). An explanation suggests that if the medium is ideal for generating neutrinos it will absorb the pionic gamma photons that are also generated in the hadronic process (Halzen et al. 2019; Halzen & Kheirandish 2020). Observations suggest this might actually be true, at least for a population of high-energy (HE) neutrino sources (e.g., Kun et al. 2021). The PeV–neutrino event called Big Bird arrived in the local minimum of a ∼85-day-long dip of the source candidate blazar PKS B1424-418 (Kadler et al. 2016). The neutrino event IC-190730A was detected from the direction of PKS 1502+106 when this blazar was in the middle of a deep, about 180-day-long gamma-ray minimum (Kun et al. 2021). The neutrino event IC-170922A was observed in a ∼28-day-long gamma-ray dip of TXS 0506+056.

As is discussed in Kun et al. (2021), temporal gamma-suppression potentially resolves the apparent contradiction of the blazar models simultaneously producing a detectable neutrino flux and a gamma flare, since at the time of efficient neutrino production the observed gamma photon flux drops, because the energy of pionic gamma rays cascades down to X-ray energies.

In this paper we examine eight Fermi-LAT blazars with good photon statistics (Kun et al. 2022) that we found within the 90% containment area of the 70 IceCube track-type neutrino events. In Sect. 2 we introduce our previous search to find neutrino-source candidates in the Fermi-LAT 4FGL-DR2 catalog, and summarize the source properties of individual neutrino-source candidates. In Sect. 3 we present our Fermi-LAT light curve analysis. In Sect. 4 we present our results, while in Sect. 5 we discuss them and give our conclusions.

2. Source selection

2.1. Previous catalog search to find Fermi-LAT 4FGL-DR2 neutrino-source candidates

In Kun et al. (2022), we searched gamma-ray, X-ray, and radio point sources in the 90% containment area of the 70 IceCube track-type neutrino events recorded between October 13, 2009 and September 30, 2019 and listed by Giommi et al. (2020). We found 29 gamma-ray neutrino-source candidates in the Fermi-LAT 4FGL-DR2 catalog (Abdollahi et al. 2020), 61 such X-ray point sources in the Swift-XRT Point Source catalog (keeping only the AGN-type objects, Evans et al. 2020), and 87(96) neutrino counterparts in the Combined Radio All-Sky Targeted Eight GHz Survey (Healey et al. 2007) at 4.8 GHz (8.4 GHz). These are the so-called ν-samples.

After the construction of the ν-samples, we conducted a randomness test of them to see if our selection method chooses random ν-samples or ν-samples with flux properties that are significantly different from the full catalogs. We note that we did not scramble the declinations of the neutrinos because of the strong variation in the sensitivity of IceCube with declination. We found that the constructed ν-samples represent special samples of sources taken from the full Fermi-LAT 4FGL-DR2/Swift-XRT 2SXPS/CRATES catalogs with similar significance (2.1σ, 1.2σ, 2σ at 4.8 GHz and 2.1σ at 8.4 GHz, respectively). In this sense, special means that the neutrinos do not randomly select samples. We already analyzed the neutrino-radio connection in Kun et al. (2022). After collecting redshifts and deriving subsamples of the CRATES catalog complete in the redshift–luminosity plane, we found that the 4.8 GHz (8.4 GHz) subsample can explain (90% C.L.) between 4% and 53% (3% and 42%) of the 70 IceCube track-type neutrino events listed by Giommi et al. (2020).

In this paper we focus on the possible neutrino-gamma ray connection. We list the name, the sky coordinates, the flux properties, the class, and the association of the 29 neutrino-source candidate Fermi-LAT blazars in Table 1. We generated the gamma-ray light curves plus and minus half a year around the detection time of the associated neutrinos. For sources with integrated flux below 10⁻⁹ph cm⁻² s⁻¹ in the energy range from 0.1 to 100 GeV, and not experiencing bright flaring activity in the time window, there are not enough photon statistics to obtain good resolution on the light curve. We therefore applied the above gamma-ray photon flux threshold to obtain a sample of eight candidate sources (associated with seven neutrino events) from this selection that we study in the following. In the next subsection we introduce the eight neutrino-source candidate Fermi-LAT blazars that constitute our sample of interest in this study.

2.2. Individual neutrino-source candidates

RBS 0958 is a BL Lac object at redshift z = 0.13793 ± 0.00004 (Ahn et al. 2012), associated with the Fermi-LAT source 4FGL J1117.0+2013. RBS 0958 is a high-synchrotron peaked (e.g., Paiano et al. 2017) VHE candidate (Furniss et al. 2015) and is considered to be a candidate TeV emitter based on the position of its synchrotron peak (Costamante et al. 2018). RBS 0958 was observed by the High Energy Stereoscopic System (H.E.S.S.) as having an integral flux upper limit of I(> 610 GeV) = 1.44 × 10⁻¹² ph cm⁻² s⁻¹ (Aharonian 2015).

GB6 J1040+0617 is a low-synchrotron peaked BL Lac object at redshift z = 0.7351 ± 0.0045 (Ahn et al. 2012; Maselli et al. 2015), associated with the Fermi-LAT source 4FGL J1040.5+0617. Garrappa et al. (2019a) highlights the increased gamma-ray activity of this source for a period of ∼90 days around the time of the neutrino detection. Gabányi et al. (2019) analyzed 3.4 μm and 4.6 μm band GB6 J1040+0617 data of the Wide-field Infrared Survey Explorer (WISE). They find the source was brightest, increasing its flux density by almost 90%, in the mission phase preceding the neutrino event IceCube-141209A by 4.5 days.

PKS 1313-333 is a flat-spectrum radio quasar (FSRQ) at redshift z ≈ 1.210 (Jauncey et al. 1982), associated with the Fermi-LAT source 4FGL J1316.1-3338. Fermi-LAT detected an increasing gamma-ray flux from a source positionally consistent with PKS 1313-333¹, which was confirmed by Astro-rivelatore Gamma a Immagini Leggero (AGILE)². Integrating photons from 02:00 UT on January 10, 2016 to 02:00 UT on January 12, 2016, a maximum likelihood analysis of AGILE data yields the detection of the source at a significance level of about 5σ, with a flux of (1.6 ± 0.6)×10⁻⁶ ph cm⁻² s⁻¹ (E > 100 MeV), in agreement with the Fermi-LAT measurement. An increasing gamma-ray flux was observed again in 2021³.

TXS 0506+056 is a BL Lac object (suggested to have at least temporarily intrinsic properties of a FSRQ, Padovani et al. 2019) at redshift z = 0.3365 ± 0.0010 (Paiano et al. 2018), associated with the Fermi-LAT source 4FGL J0509.4+0542. This source is a TeV⁴ emitter. A multimessenger campaign, triggered by the coincident observation of a gamma-ray flare and the 290 TeV IceCube neutrino IC-170922A (assuming a neutrino spectrum of ${\epsilon }_{\nu }^{-2.19}$), identified TXS 0506+056 as the possible source of high-energy neutrinos at a 3σ significance level (IceCube Collaboration et al. 2018). With this, TXS 0506+056 became the first identified source of high-energy cosmic neutrinos at such a high significance. The blazar TXS 0506+056 is 0.1° away from the best-fit neutrino location of IC-170922A (56.5% signalness). The quantity signalness is the probability that an observed neutrino has an astrophysical (s = 1) or atmospheric (s = 0) origin (see its definition in Aartsen et al. 2017). Going back to the archival data, the IceCube Collaboration found a 160-day-long neutrino flare in 2014–2015 from the direction of TXS 0506+056 at a significance level of 3.5σ (IceCube Collaboration 2018). Despite being more dominant than IC-170922A, this neutrino flare happened during a low-gamma state of TXS 0506+056 (Garrappa et al. 2019a). Putting this behavior into a common physical picture with IC-170922A and the coincident gamma-flare still challenges the current blazar models.

The Major Atmospheric Gamma Imaging Cherenkov Telescopes (MAGIC) measured the upper limit on the > 90 GeV flux of TXS 0506+056 as < 3.56 × 10⁻¹¹ ph cm⁻² s⁻¹ in the epoch MJD 58020.18 (Ansoldi et al. 2018), about 1 day after the detection of IC-170922A, and MAGIC measured a significant increase in > 90 GeV flux, (4.7 ± 1.4)×10⁻¹¹ ph cm⁻² s⁻¹ at MJD 58029.22, and (8.7 ± 2.0)×10⁻¹¹ ph cm⁻² s⁻¹ at MJD 58030.24, 10–11 days after the detection of IC-170922A. Confirming this elevated high-energy state, the Very Energetic Radiation Imaging Telescope Array System (VERITAS) Collaboration reported the integral flux (> 110 GeV) of TXS 0506+056 as (8.9 ± 1.6)×10⁻¹² ph cm⁻² s⁻¹ (Abeysekara et al. 2018). H.E.S.S. performed follow-up observations of the high-energy neutrino on September 23, 2017 and the following night (September 24, 2017 at 03:10 UTC). A preliminary on-site calibration and analysis searching for a point-like gamma-ray source from within the 90% uncertainty region of the neutrino event IC-170922A revealed no significant detection (upper limit < 1.0 × 10⁻¹¹ ph cm⁻² s⁻¹ for > 175 GeV, Schüssler et al. 2019). Employing the ANTARES neutrino telescope data recorded in a cone of 3°, centered on the best-fit coordinates of IC-170922A and within a ±1 h time-window centered on the event time, Albert et al. (2018) constrained the neutrino fluence to < 15 GeV cm⁻² for γ = 2 (integrated between 3.3 TeV and 3.4 PeV) and to < 34 GeV cm⁻² for γ = 2.5 (integrated between 450 GeV and 280 TeV), assuming a point source with power-law spectrum dN/dE ∼ E^−γ.

The radio flux density of TXS 0506+056 measured by the Owens Valley 40m Radio telescope started to increase in the beginning of 2016, and was still increasing after the neutrino detection. During the 2014–2015 neutrino flare, the radio flux density was average. With very long baseline interferometry (VLBI), the radio jet of TXS 0506+056 shows strong signs of spine-sheath structure within the inner 1 mas, corresponding to about 70–140 pc in the deprojected distance from the millimeter-VLBI core (Ros et al. 2020). It has been shown that at 15 GHz (e.g., Kun et al. 2019) and 43 GHz (Ros et al. 2020) the VLBI core is responsible for the radio brightening. Several scenarios have been suggested to explain the neutrino emission (Britzen et al. 2019; Garrappa et al. 2019a; Reimer et al. 2019; Petropoulou et al. 2020; Zhang et al. 2020).

PKS 1454-354 is a low-synchrotron peaked FSRQ type blazar at redshift z ≈ 1.424 (Jackson et al. 2002), associated with the Fermi-LAT source 4FGL J1457.4-3539. This source is located within the 90% error ellipse of the neutrino IC-181014A. The Fermi-LAT collaboration reported an increase in gamma-ray activity from the source in September 2008 (Abdo et al. 2009). PKS 1454-354 has a weak one-sided jet.

NVSS J042025-374443 is a radio source at the (photometric) redshift z ≈ 0.3 (Malkin 2018), associated with the Fermi-LAT source 4FGL J0420.3-3745. This is another > 100 MeV gamma-ray source located within the 90% IC-190504A localization error, at a distance of roughly 37 arcmin from the best-fit neutrino localization (PKS 0426-380 is located at a distance of roughly 72 arcmin).

PKS 0426-380 is a reclassified FSRQ at redshift z ≈ 1.111 (Heidt et al. 2004; Sbarufatti et al. 2005), associated with the Fermi-LAT source 4FGL J0428.6-3756. This source is located within the 90% error ellipse of the neutrino IC-190504A. Tanaka et al. (2013) reported the Fermi-LAT detection of two very-high-energy (E ∼ 100 GeV) gamma-ray photons from the directional vicinity of PKS 0426-380. Zhang et al. (2017) find a possible gamma-ray quasi-periodicity of this blazar with a period of 3.35 ± 0.68 yr at the significance level of approximately 3.6σ.

PKS 1502+106 is an FSRQ, the fifteenth brightest gamma-ray source at > 100 MeV in terms of energy flux among 2, 863 sources in the fourth catalog of AGN detected by Fermi-LAT (4LAC Ajello et al. 2020). The redshift of PKS 1502+106 was estimated to be z = 1.839 (Smith et al. 1977), later confirmed by the good signal-to-noise ratio spectrum of the Sloan Digital Sky Survey (SDSS; z = 1.8385 ± 0.0024 at high confidence, Schneider et al. 2010). Given its distance and still-enormous observed gamma-ray photon flux, PKS 1502+106 must have extremely high intrinsic luminosity. This source is highly variable in the Fermi-LAT gamma-ray band (Abdo et al. 2010a). It showed strong gamma-ray flares in the Fermi-LAT energy range in 2008 (Abdo et al. 2010b). The flaring activity of PKS 1502+106 was renewed between 2009 and 2010. From 2015, the source showed more strong gamma-ray flares (see e.g., in Kun et al. 2021).

VLBI observations reveal a one-sided, curved radio jet and a bright core structure at parsec scales (e.g., Karamanavis et al. 2016). At radio frequencies, the source is also highly variable.

The IceCube detector recorded a high-energy neutrino event on June 30, 2019 (with 67% signalness) that was distributed through the Gold alert stream (IceCube Collaboration 2019). The most probable neutrino energy was approximately 300 TeV, assuming a ${\epsilon }_{\nu }^{-2.19}$ neutrino spectrum. PKS 1502+106 lies with an offset of 0.31° from the best-fit neutrino location of IC-190730A. At the time of the corresponding neutrino detection, PKS 1502+106 was in a quiet state in the optical (Stein et al. 2019; Lipunov et al. 2019), UV (Franckowiak et al. 2020), X-ray (Santander et al. 2019; Kun et al. 2021), and gamma-ray regimes (Garrappa et al. 2019b). A long-term radio outburst that started in 2014 had just reached its peak flux when the IceCube neutrino was recorded (Kiehlmann et al. 2019). This source became part of the Icecube and ANTARES a priori defined monitored source list (34 gamma-ray bright AGN).

After the detection of IC-190730A, the nearby candidate PKS 1502+106 quickly became the center of focus of many studies. Kun et al. (2021) studied cosmic neutrinos from temporarily gamma-suppressed blazars, such as PKS 1502+106. As it is shown, the neutrino event IC-190730A (IceCube Collaboration 2019) arrived during an exceptionally low gamma-ray state of PKS 1502+106 that lasted for several months. Rodrigues et al. (2021) identify one quiescent state and two flaring states with hard and soft gamma-ray spectra. They find two hadronic models that can describe the multiwavelength emission during all three states: a leptohadronic model with a contribution from photohadronic processes to X-rays and gamma rays, and a proton synchrotron model where the emission from keV to 10 GeV comes from proton synchrotron radiation. Oikonomou et al. (2021) analyzed UV, optical, and X-ray observations of PKS 1502+106, as well as collecting other multiwavelength data. They find that the leptohadronic modeling, in which the multiwavelength emission of PKS 1502+106 originates beyond the broad-line region (BLR) and inside the dust torus, is most consistent with the observations. Analyzing 15 GHz Very Long Baseline Array (VLBA) and astrometric 8 GHz VLBA data, Britzen et al. (2021) present evidence for a radio ring structure that develops with time. They suggest that the neutrino IC-190730A was most likely produced by proton-proton interaction in the blazar zone beyond the BLR, enabled by episodic encounters of the jet with dense clouds, that is to say some molecular clouds in the narrow line region.

3. Analysis of Fermi-LAT data of eight blazars

The LAT instrument on board the Fermi Gamma-ray Space Telescope is designed to cover the energy band from 20 MeV to greater than 300 GeV. In this study we employed one year of Pass8 Fermi-LAT⁵ data⁶ for each of the eight Fermi-LAT neutrino-source candidates, in a time range plus and minus half a year centered on the detection time of the corresponding neutrinos listed in Table 2. The region of interest (ROI) was centered at the J2000 sky coordinates of the 4FGL-DR2 sources, encompassing an area of the sky within 15° ROI radius. We selected event type “front+back” (evtype = 3), which is the recommended type for a point source analysis.

We generated the likelihood light curve of the eight brightest 4FGL-DR2 neutrino-source candidates. We performed a binned likelihood analysis of the data utilizing the fermipy v1.0.1 and ScienceTools v2.0.8 packages, both built in the FermiBottle Docker container and analysis environment provided by the Fermi Science Support Center⁷. The analysis was carried out on a cluster of 16 vCPU-s (Intel Skylake 16 × 2.2 GHz) in the ELKH Cloud⁸. The instrument response functions +P8R3SOURCEV2+ were employed together with templates of the Galactic interstellar emission model +glliemv07.fits+ and of the isotropic diffuse emission +isoP8R3SOURCEV2v1.txt+⁹. We applied the nominal data quality cut +(DATAQUAL > 0) (LATCONFIG==1)+, and a zenith angle cut θ < 90° to eliminate Earth limb events. Time intervals when the Sun was closer than 15° to the target sources were filtered out.

The detection level of each source was measured by the definition of test statistic (TS, Mattox et al. 1996), and the detection limit of new sources was set to TS_min = 25 (∼5σ). The test statistic is defined as TS = −2ln(L_max, 1/L_max, 0), where L_max, 0 is the maximum likelihood value for a model without an additional source (the “null hypothesis”) and L_max, 1 is the maximum likelihood value for a model with the additional source at a specified location.

Since the two models are nested, according to Wilks’ theorem (Wilks 1938) the TS-distribution will asymptotically follow a χ² distribution. Then the TS distribution should be drawn from a ${\chi }_n^2$ distribution, where n is the difference in the degree of freedom between the models with and without the additional source. In our Fermi-LAT analysis the normalization factor of the Galactic and isotropic diffuse components and sources with TS > 25 were freed. We associated power-law profiles to the additional sources that are detected in our analysis but that do not appear in the 4FGL-DR2. We consider new point sources to be separate sources if the angular distance between them is larger then 0.3 degrees. We present a summary of additional sources found in the ROI of seven target sources in Table A.1. To get more information about variability in the light curves, we applied the adaptive binning algorithm (Lott et al. 2012), setting the relative flux uncertainty to 15%. With adaptive binning, our analysis is more sensitive to the light variations. To achieve the same photon statistics throughout the light curve, the bins are shorter when the source is bright and the bins are wider when the source is fainter.

We repeated the analysis with 3–4 days fixed binning for two sources that have the best photon statistics in our sample, namely TXS 0506+056 and PKS 1502+106. We show the light curve for these sources with fixed binning overplotted on their adaptively binned light curve in Appendix C (see Figs. C.1 and C.2, respectively). The shapes of the light curves with fixed binning are similar to the adaptively binned light curves, but with fixed binning the light curves are more noisy when the sources were faint. In comparing the results from the two light curve generation methods, we do not see different structures that would alter the results of this paper. Moreover, fixed binning of the light curve of the other sources in our sample, being fainter than the two shown here, would lead to a large amount of non-detections. This further justifies the choice of the adaptive binning method employed in our work.

4. Results

Below, we refer to the centers and widths of the bins in the form of center ±width/2 in order to help understand the time characteristics of the bins. We note here that the ± does not mean an error. We modeled the light curves of the neutrino-source candidates between 100 MeV and 300 GeV, except for GB6 J1040+0617, PKS 1313-333, and TXS 0506+056, where, to avoid contamination by the neighboring bright gamma-ray sources, we set the minimum energy of the analysis to 300 MeV. We fit the normalization factor and the spectral index of the source of interest in each bin. Normalization of sources with TS-value > 25 and sources within 5 degrees of the source of interest was also free (regardless of their test statistic). To better understand significant variations in the light curves, we applied the Bayesian block algorithm of Scargle et al. (2013) to them. Following Garrappa et al. (2019a), we set the probability to p = 0.05. The Bayesian blocks are shown on the light curve of the gamma-ray sources in the appendix, as stated below, from source to source. We summarize the characteristic timing and flux values of minimum and maximum photon flux bins compared to the neutrino bins in Table 3, and elaborate more on the results below.

RBS 0958 is located 1.19° away from the best-fit localization of the neutrino event IC-130408A, within the 90% error ellipse of the neutrino. We show the one-year-long likelihood light curve of RBS 0958 in Fig. B.1 along with the photon index as function of time. The corresponding neutrino was detected in T_ν = MJD 56390.1887627, with reconstructed neutrino energy of $E_{\nu }=30.8_{-3.5}^{+3.3}$ TeV. The neutrino came in the Fermi-LAT bin MJD 56417.0 ± 29.4 days, with an average flux of (1.46 ± 0.58)×10⁻⁸ ph cm⁻² s⁻¹. This bin is situated on the decreasing branch of a gamma flare, which is centered to MJD 56303.2 ± 10.0 days (86.97 days before the neutrino detection). Though the Bayesian blocks (with p = 0.05) do not indicate a significantly variable flux during the one-year-long time window (see Fig. B.1), the peak photon flux of the suspected flare ((8.0 ± 2.4)×10⁻⁸ ph cm⁻² s⁻¹) is about 45–83 times higher than the average photon flux of RBS 0958 in the 4FGL-DR2 catalog (1.24 × 10⁻⁹ ph cm⁻² s⁻¹). This indicates that there might indeed be a gamma-ray flare going on in this window. Based on its photon indices, RBS 0958 shows signs of spectral softening during the gamma-ray flare. The peak flux of this flare was (7.97 ± 2.35)×10⁻⁸ ph cm⁻² s⁻¹, which is about 5.5 times the flux of the bin in the gamma-ray light curve containing the detected neutrino (from here on we refer to this specific bin of the gamma-ray light curve as “neutrino bin”). The bins before (MJD 56370.6 ± 17.0) and after (MJD 56477.0 ± 30.6) the neutrino bin have similar fluxes to the flux of the neutrino bin ((1.65 ± 1.02)×10⁻⁸ ph cm⁻² s⁻¹ and (1.08 ± 0.59)×10⁻⁸ ph cm⁻² s⁻¹, respectively).

GB6 J1040+0617 is located within the 90% error ellipse of IC-141209A, 0.45° away from the best-fit localization of the neutrino. We show the one-year-long likelihood light curve of GB6 J1040+0617 in Fig. B.2 along with the photon index as a function of time. The corresponding neutrino was detected in T_ν = MJD 57000.14311, with reconstructed neutrino energy of $E_{\nu }=97.4_{-9.6}^{+9.6}$ TeV. The neutrino came in the Fermi-LAT bin MJD 56994.9 ± 22.8 days, with an average flux of (1.57 ± 0.46)×10⁻⁸ ph cm⁻² s⁻¹. This is the peak flux of the gamma-ray flare shown by this source during the one-year-long light curve centered on the time of the neutrino detection. The flux in the previous bin (MJD 56949.6 ± 22.5) is about half the flux of the neutrino bin ((7.9 ± 4.1)×10⁻⁹ ph cm⁻² s⁻¹). The flux in the next bin (MJD 57038.2 ± 20.5) is about 2/3 of the flux in the neutrino bin ((1.11 ± 0.38)×10⁻⁸ ph cm⁻² s⁻¹).

PKS 1313-333 is located within the 90% error ellipse of IC-160814A (T_ν = MJD 57614.91), 1.29° away from its best-fit sky coordinates. We show the one-year-long likelihood light curve of PKS 1313-333 in Fig. B.3 along with the photon index as a function of time. The neutrino came in the Fermi-LAT bin MJD 57632.6 ± 51.0, with flux of (5.9 ± 2.4)×10⁻⁹ ph cm⁻² s⁻¹. The previous bin, which is the peak of the light curve, is situated at MJD 57571.6 ± 10.0 with a flux of (2.57 ± 0.91)×10⁻⁸ ph cm⁻² s⁻¹, which is more than four times the flux of the neutrino bin. The next bin (MJD 57734.7 ± 51.0) has an average flux of (2.3 ± 2.0)×10⁻⁹ ph cm⁻² s⁻¹, about half the flux of the neutrino bin. This is also the faintest bin of the light curve.

TXS 0506+056 is located within the 90% error ellipse of IC-170922A (T_ν = MJD 58018.87), 0.1° away from its best-fit sky coordinates. We show the one-year-long likelihood light curve of TXS 0506+056 in Fig. B.4 along with the photon index as a function of time. The neutrino came in the 2.3-day-long Fermi-LAT bin centered to MJD 58018.0, which has a flux of (1.10  ±  0.32)×10⁻⁷ ph cm⁻² s⁻¹. This bin is on the decreasing side of a gamma-ray dip. The highest-flux bin (MJD 57938.0 ± 0.8) has a flux of (3.23 ± 0.64)×10⁻⁷ ph cm⁻² s⁻¹, which is about three times the flux in the neutrino bin. The highest-flux bin is situated 80.0 days before the detection time of the neutrino. The minimum flux of the light curve (2.2 ± 0.1)×10⁻⁸ ph cm⁻² s⁻¹ is reached in the bin MJD 57962.8 ± 3.3, centered 155 days prior to the detection.

PKS 1454-354 is located within the 90% error ellipse of IC-181014A (T_ν = MJD 58405.5), 0.85° away from best-fit sky coordinates of this neutrino event. We show the one-year-long likelihood light curve of PKS 1454-354 in Fig. B.5 along with the photon index as a function of time. The neutrino came in the Fermi-LAT bin MJD 58396.7 ± 22.8, with the photon flux of (2.7 ± 1.8)×10⁻⁸ ph cm⁻² s⁻¹, indicating the source was in a faint phase. We note that the photon flux in the neutrino bin is about an order of magnitude larger than the average photon flux of this source in the 4FGL-DR2 catalog (3.49 × 10⁻⁹ ph cm⁻² s⁻¹), and that it is about an order of magnitude smaller than the peak photon flux of this source in the same catalog (2.56 × 10⁻⁷ ph cm⁻² s⁻¹).

NVSS J042025-374443 is located within the 90% error ellipse of IC-190504A (T_ν = MJD 58607.77), 0.76° away from its best-fit sky coordinates. We show the one-year-long likelihood light curve of NVSS J042025-374443 in Fig. B.6 along with the photon index as a function of time. Though the source is rather faint in the studied observational window, permitting us to get only three bins of the light curve, we can note the following. The neutrino came in the Fermi-LAT bin MJD 58575.1 ± 60, with a flux of (9.5 ± 4.9)×10⁻⁹ ph cm⁻² s⁻¹. The previous bin is situated as MJD 58470.1 ± 45.0, with a flux of (1.36 ± 0.52)×10⁻⁸ ph cm⁻² s⁻¹. The next bin, MJD 58710.1 ± 75, has an average flux of (7.4 ± 3.9)×10⁻⁹ ph cm⁻² s⁻¹. The highest-flux bin, MJD 58470.1 ± 45.0, is centered 137.63 days before the detection time of the corresponding neutrino and has an average flux of (1.36 ± 0.52)×10⁻⁸ ph cm⁻² s⁻¹. We note that the photon flux of the neutrino bin is 4–13 times higher than the average photon flux of this source in the 4FGL-DR2 catalog (1.09 × 10⁻⁹ ph cm⁻² s⁻¹), and about 2–8 times smaller than the peak photon flux of the source in the 4FGL-DR2 catalog (3.56 × 10⁻⁸ ph cm⁻² s⁻¹).

PKS 0426-380 is also located within the 90% error ellipse of IC-190504A (T_ν = MJD 58607.77), 1.47° away from its best-fit sky positions. We show the one-year-long likelihood light curve of PKS 0426-380 in Fig. B.7 along with the photon index as a function of time. The neutrino came in the Fermi-LAT bin MJD 58607.8 ± 2.7 that has an average flux of (1.97 ± 0.41)×10⁻⁷ ph cm⁻² s⁻¹. The highest-flux gamma-ray bin MJD 58651.3 ± 1.6 is centered 43.6 days after the neutrino, and the peak of a gamma-ray flare shown by the source. This light curve bin has an average flux of (3.76 ± 1.02)×10⁻⁷ ph cm⁻² s⁻¹, which is about two times higher than the flux of the neutrino bin. PKS 0426-380 is a significantly stronger gamma source than the other gamma source, NVSS J042025-374443, that also lies within the 90% error ellipse of IC-190504A.

PKS 1502+106 is located within the 90% error ellipse of IC-190730A (T_ν = MJD 58694.87), 0.58° away from the best-fit sky coordinates of the corresponding neutrino. This 299 TeV event was distributed through the Gold channel of the IceCube real-time alert system. Its signalness is 0.67 and the false alarm rate is 0.677 yr⁻¹. We show the one-year likelihood light curve of PKS 1502+106 in Fig. B.8 along with the photon index as a function of time. The neutrino came in the Fermi-LAT bin MJD 58703.2 ± 25.5, which has a flux of (5.4 ± 1.0)×10⁻⁸ ph cm⁻² s⁻¹. This bin marks the global minimum of the one-year light curve of PKS 1502-106, the lowest-flux bin of a prolonged deep minimum lasting about 178 days. The highest-flux bin (MJD 58590.9 ± 7.0), centered 112 days before the neutrino detection, has an average flux of (3.4 ± 0.8)×10⁻⁷ ph cm⁻² s⁻¹), which is more than six times the flux of the neutrino bin.

5. Discussion and conclusions

High-energy, low-optical-depth neutrino formation processes predict that the neutrino flux is proportional to the gamma-ray photon flux and one therefore expects the emission of high-energy neutrinos during gamma flares.

We calculated the chance probability of detecting a neutrino in a bin with a gamma-ray photon flux higher then the observed gamma-ray photon flux of the neutrino bin:

$$\begin{array}{c}\hfill p=P(F_{\gamma ,max}>F_{\gamma ,\nu })=\frac{{\mathrm{\Sigma }}_i{\int }_{F_{\gamma ,\nu }}^{\infty }\mathcal{N}(x,F_i,{\sigma }_i)\mathrm{d}Fx}{{\mathrm{\Sigma }}_i{\int }_{-\infty }^{\infty }\mathcal{N}(x,F_i,{\sigma }_i)\mathrm{d}Fx}\end{array}$$(1)

where 𝒩(x, F_i, σ_i) is the Gaussian evaluated at x, when we assume the uncertainty of the flux is normally distributed. When the p-value is close to zero, it means that the neutrino came in gamma maximum and there is no chance of getting higher gamma-ray photon flux during the studied time range. This is the case when the medium is transparent to the gamma photons and we see a gamma peak during the detection of high-energy neutrinos from the source.

When the p-value is close to unity, it means that the neutrino came in a gamma minimum and therefore all bins of the light curve are brighter than the flux of the bin in which the neutrino arrived. This is the case when the neutrino flux and gamma-ray photon flux are anti-correlated, for example due to the high in-source gamma-gamma optical depth, and at the time of maximum neutrino emissivity the pionic gamma photons are maximally absorbed and cascade down to X-ray energies. This is the case of beam dampening. Using the above considerations we calculated the p-values (see in Table 3).

It seems that the corresponding neutrino came in the global maximum of the one-year-long gamma light curve of GB6 J1040+0617, and in global minimum of the one-year long gamma light curve of PKS 1502+106. For the other neutrino-source candidates, PKS 1313-333, PKS 1454-354 NVSS J042025-374443, RBS 0958, TXS 0506+056, and PKS 0426-380, the neutrino arrival time seems to not correlate with the gamma-ray photon flux, even though fainter phases are more characteristic of the neutrino arrivals.

What we can conclude is that, apart from GB6 J1040+0617, none of the neutrino-source candidates were significantly bright in a time window of ±6 months around the corresponding IceCube neutrinos. We note here that our findings are based on a one-year-long gamma-ray photon flux curve of the respective sources. Taking into account the whole Fermi-LAT light curves might alter the outcomes of this analysis.

The fact that, out of the eight neutrino candidate sources (associated with seven neutrino events) studied, we find only one neutrino-source candidate blazar at peak gamma-ray photon flux during the neutrino emission means that either the sources in our sample are not the sources of the corresponding IceCube neutrinos, or an in-source effect (e.g., the suppression of gamma rays due to high gamma-gamma opacity, Kun et al. 2021) complicates the multimessenger scenario of neutrino emission, at least for these blazars. It is also possible that a mix of these two scenarios is at play. Future studies conducted at even higher energies, at which gamma-ray absorption should be more evident, or alternative measurements, such as HI column density measurements, might be helpful to disambiguate between these scenarios.

Acknowledgments

We thank Simone Garrappa, Deirdre Horan, Josefa Becerra Gonzalez, Melissa Pesce Rollins and Patrik Milán Veres for comments and discussions. We acknowledge support from the Deutsche Forschungsgemeinschaft DFG, within the Collaborative Research Center SFB1491 “Cosmic Interacting Matters – From Source to Signal” (project No. 445052434). E.K. thanks the Alexander von Humboldt Foundation for its Postdoctoral Fellowship. I.B. acknowledges the support of the Alfred P. Sloan Foundation and NSF grants PHY-1911796 and PHY-2110060. This paper makes use of publicly available Fermi-LAT data provided online by the https://fermi.gsfc.nasa.gov/ssc/data/access/Fermi Science Support Center. On behalf of Project ‘fermi-agn’ we thank for the usage of the ELKH Cloud that significantly helped us achieving the results published in this paper.

References

Appendix A: Additional sources found in the ROI of the target sources

Appendix B: Light curves

Fig. B.1. Gamma-ray flux variation of RBS 0958. Upper panel: likelihood light curve of RBS 0958 with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.2. Gamma-ray flux variation of Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.3. Gamma-ray flux variation of PKS 1313-333. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.4. Gamma-ray flux variation of TXS 0506+056. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.5. Gamma-ray flux variation of PKS 1454-354. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.6. Gamma-ray flux variation of NVSS J042025-374443. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.7. Gamma-ray flux variation of PKS 0426-380. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.

Fig. B.8. Gamma-ray flux variation of PKS 1502+106. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as function of time.

Appendix C: Comparison of adaptive and fixed binning of the light curves of two bright sources

Fig. C.1. Light curve of TXS 0506+056 with fixed binning (green crosses) overplotted on its adaptively binned light curve (orange crosses). For very faint bins of the fixed binning light curve we plotted only the upper limits (down-pointing green arrows).

Fig. C.2. Light curve of PKS 1502+106 with fixed binning (green crosses) overplotted on its adaptively binned light curve (orange crosses). For very faint bins of the fixed binning light curve we plotted only the upper limits (down-pointing green arrows).

All Tables

All Figures

	Fig. B.1. Gamma-ray flux variation of RBS 0958. Upper panel: likelihood light curve of RBS 0958 with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.2. Gamma-ray flux variation of Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.3. Gamma-ray flux variation of PKS 1313-333. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.4. Gamma-ray flux variation of TXS 0506+056. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.5. Gamma-ray flux variation of PKS 1454-354. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.6. Gamma-ray flux variation of NVSS J042025-374443. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.7. Gamma-ray flux variation of PKS 0426-380. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as a function of time.
In the text

	Fig. B.8. Gamma-ray flux variation of PKS 1502+106. Upper panel: likelihood light curve with Bayesian blocks (p = 0.05). The continuous green line represents the Bayesian block binning. Lower panel: photon index as function of time.
In the text

	Fig. C.1. Light curve of TXS 0506+056 with fixed binning (green crosses) overplotted on its adaptively binned light curve (orange crosses). For very faint bins of the fixed binning light curve we plotted only the upper limits (down-pointing green arrows).
In the text

	Fig. C.2. Light curve of PKS 1502+106 with fixed binning (green crosses) overplotted on its adaptively binned light curve (orange crosses). For very faint bins of the fixed binning light curve we plotted only the upper limits (down-pointing green arrows).
In the text