Cosmic Neutrinos from Temporarily Gamma-suppressed Blazars

The origin of the isotropic cosmic flux of high-energy neutrinos, discovered by IceCube in 2013 (IceCube Collaboration 2013a, 2013b, 2020), is still unknown. The most promising strategy to identify the underlying sources has been multimessenger observations (IceCube Collaboration et al. 2018; Murase & Bartos 2019). Studies aiming to find sources of cosmic neutrinos have recently focused on the blazar subclass of active galactic nuclei (AGNs; e.g., Kadler et al. 2016; Krauß et al. 2018; Franckowiak et al. 2020; Giommi et al. 2020; Hoerbe et al. 2020), extreme particle accelerators directing their relativistic γ-ray jet toward Earth.

The spectral energy distribution of blazars is characterized by a low- and a high-energy bump (e.g., Abdo et al. 2010). The low-energy component is due to synchrotron radiation of charged particles gyrating about the magnetic field of the AGN, which is the primary process generating the radio continuum of jetted AGN (e.g., Boettcher et al. 2012). The origin of the high-energy emission is less clear, as the observed spectrum can be explained by leptonic, hadronic, or hybrid lepto-hadronic models (e.g., Murase et al. 2018; Gao et al. 2019; Rodrigues et al. 2019). In leptonic models the synchrotron self-Compton process or external inverse Compton fields, while in hadronic models synchrotron emission of high-energy protons and the decay of neutral pions can lead to the enhanced γ-ray emission of blazars (e.g., Biermann et al. 2011; Aharonian et al. 2013). Above 100 GeV, the emitted radiation is attenuated by extragalactic background light (e.g., Gilmore et al. 2012).

A multimessenger campaign initiated by the observation of a muon neutrino of 290 TeV energy (IC-170922A) identified the γ-ray flaring blazar TXS 0506+056 (IceCube Collaboration et al. 2018) as the source of the neutrino. It was argued that the rapid flux density variation of TXS 0506+056 at the arrival time of the neutrino came from the VLBI core at 15 GHz (Kun et al. 2019) and at 43 GHz (Ros et al. 2020). Knowing where to look, IceCube retrieved a second neutrino flare from the archival data producing 13 ± 5 neutrinos over 158 days between 2014 and 2015 (IceCube Collaboration 2018). The flare dominates the integrated flux of 10 yr of IceCube data, leaving the flare containing IC-170922A as a distant second. In contrast to IC-170922A, the γ-activity of the source was low at the time of the 2014/2015 neutrino flare (Garrappa et al. 2019). Halzen et al. (2019) suggested this could be the result of large optical depths for photomeson production required by efficient high-energy neutrino production, which in turn prevents the escape of high-energy pionic γ-rays. A population of γ-ray dark cosmic-ray accelerators as a neutrino origin was already suspected (e.g., Murase et al. 2016).

IceCube recorded its highest-energy cosmic neutrino alert ever on 2019 July 30, IC-190730A (IceCube Collaboration 2019), triggered by a muon neutrino with an estimated energy of 300 TeV. The Fermi-LAT strong point source 4FGL J1504.4+1029, associated with the active galaxy PKS 1502+106 (z = 1.838; Pâris et al. 2017), is located within 0.31 degrees from the best-fit neutrino location. 4FGL J1504.4+1029 is one of the brightest γ-ray sources. Given the large redshift of 1.844, the source possesses an extremely high intrinsic luminosity. Kiehlmann et al. (2019) reported that the flux density of PKS 1502+106 at 15 GHz measured with the OVRO 40 m Telescope (Richards et al. 2011) shows a long-term outburst peaking at the detection time of the neutrino. The VLBI flux density of PKS 1502+106 is dominated by the core at 15, 43, and 86 GHz observing frequencies (Karamanavis et al. 2016). Rodrigues et al. (2020) found three different activity states, one quiescent, and two flaring states distinguished by the hardness of the X-ray spectra. Recently, Britzen et al. (2021) found evidence for an evolving radio ring structure in PKS 1502+106. They also found that the radio emission is correlated with the γ-ray emission, with radio lagging the γ-rays. They suggested that the neutrino is most likely produced by pp interaction in the blazar zone.

In this Letter we examine the temporal structure of γ-ray emission for blazars TXS 0506+056, PKS 1502+106, and PKS B1424-418, which have associated high-energy neutrino counterparts. For 4FGL J1504.4+1029 (PKS 1502+106) we analyze the Fermi-LAT light curve and its correlation with the single dish radio flux density curve obtained by the OVRO 40 m Telescope (Richards et al. 2011), and emission properties at the time of the neutrino detected by IceCube. For each of the three blazars we examine, we find that γ-emission is locally suppressed at the time of the detection of their neutrino counterparts, and rapidly rises in the following weeks. Based on neutrino, γ-ray, and radio observations we present a scenario in which the neutrino emission is expected in temporary γ-suppressed states of jetted AGN.

We analyzed 12 yr of Fermi-LAT ⁸ Pass8 data ⁹ around the sky position of 4FGL J1504.4+1029 (PKS 1502+106; see also Franckowiak et al. 2020), obtained between 2008 August 4 and 2020 June 26 in the energy range from 100 MeV to 300 GeV. The search window was centered at R.A. _J2000 = 22610 and decl. _J2000 = 1049 encompassing an area of the sky within °15 ROI radius around the central coordinates. We selected event type "front+back" (evtype = 3), which is the recommended type for a point source analysis. We performed an unbinned likelihood analysis chain of the data utilizing the Fermi-LAT ScienceTools package v11r5p3 with fermipy package 0.17.3. The instrument response function P8R3_SOURCE_V2 was employed altogether with templates of the Galactic interstellar emission model gll_iem_v07.fits and of the isotropic diffuse emission iso_P8R3_SOURCE_V2_v1.txt. ¹⁰ We applied the nominal data quality cut (DATA_QUAL > 0) && (LAT_CONFIG==1), and a zenith angle cut θ < 90° to eliminate Earth limb events.

We applied a 30 day wide bin width to generate the γ light curve of PKS 1502+106, which is an unusually bright γ-ray source dominating its environment throughout the 12 yr. Therefore the input model for the light-curve generation included only PKS 1502+106 (both the prefactor and index were free parameters of the likelihood models) and diffuse components.

We show the γ-ray light curve of PKS 1502+106 in the top panel of Figure 1. The blazar exhibited larger and smaller flares throughout the Fermi observations. Notably, at the time of the detection of a coincident neutrino, its γ-flux was at a local minimum, showing rapid rise during the weeks following the minimum.

Zoom In Zoom Out Reset image size

We show the flux density curve of PKS 1502+106 obtained by the OVRO 40 m Telescope (Richards et al. 2011) at 15 GHz in the second panel of Figure 1, such that the radio flux densities are binned to the Fermi light curve. It is apparent that the γ flux and radio flux density changed simultaneously until about the end of 2014, both for high and for low flux states, and their derivatives changed similarly (mode1). After 2014, however, the radio flux density increased and stayed high while the γ flux generally decreased after a sharp rise in 2015 (mode2), such that both γ and radio changed more violently and less synchronized.

We calculated the Pearson correlation index R between the γ-ray light curve and the binned radio flux density curve in mode1. Its resulted value R = 0.85 reveals a strong, positive linear correlation. We cross-correlated the γ-ray light curve with the unbinned radio flux density curve in mode 1 by applying a a python implementation (Robertson et al. 2015) of the Discrete Correlation Function method (DCF; Edelson & Krolik 1988). The DCFs with 15, 30, 60, 90, and 120 day wide bins followed similarly shaped distributions centered at about τ ≈ 60 days. We estimated the centroid of the DCF, τ_c = (Σ_i τ_i DCF_i )/(Σ_i DCF_i ) (e.g., Raiteri et al. 2011) to be τ_c ≈ 59 days for the 60 days binning. The uncertainty of τ_c was estimated by performing Monte Carlo simulations following the "flux redistribution/random subset selection" technique (FR/RSS; Peterson et al. 1998). Time lag with uncertainties emerged at τ_c = 59 ± 32 days, pointing toward the possibility that the radio flux density curve is delayed with respect to the γ-ray light curve. Taking into account that all of the flux curves (Britzen et al. 2021) found a lag of ${151}_{-95}^{+40}$ days, similarly the radio lagging behind the γ-flux.

In mode 2 the γ-ray flux was generally decreasing after a rapid increase in the first half of 2015, while the radio flux density was at high state, suggesting a more complex behavior.

4.1. A Multimessenger Scenario

In hadronic blazar models the photomeson production is typically assumed to generate high-energy neutrinos in p γ interactions because the low particle density in the jet renders the pp process subdominant (e.g., Boettcher et al. 2012). The interaction of cosmic rays in the jet with stellar envelopes or gas clouds (e.g., Barkov et al. 2010) might, however, generate high-energy neutrinos through the pp process.

In the following, we assume that cosmic rays are proton dominated, since photodisintegration of heavier nuclei is unable to provide an effective mechanism for very-high-energy γ-ray emission from astrophysical sources (Aharonian & Taylor 2010).

The flavor-averaged diffuse high-energy cosmic neutrino flux is related to cosmic-ray flux by (Halzen et al. 2019):

$$\begin{eqnarray}&&\displaystyle \frac{1}{3}\displaystyle \sum _{\alpha }{E}_{\nu }^{2}\displaystyle \frac{{{dN}}_{\nu }}{{{dE}}_{\nu }}\simeq \displaystyle \frac{c}{4\pi }\left(\displaystyle \frac{1}{2}(1-{e}^{-{\tau }_{p\gamma }}){\xi }_{z}{t}_{H}\displaystyle \frac{{dE}}{{dt}}\right),\end{eqnarray} \tag{ 1 }$$

where the factor ξ_z accounts for the redshift evolution of sources (Ahlers & Halzen 2018). The IceCube all-flavor diffuse neutrino flux ∼ 3 × 10⁻¹¹ TeV cm⁻² s⁻¹ sr⁻¹ and the injection rate of cosmic rays above 10¹⁶ eV dE/dt ∼ 1–2 × 10⁴⁴ erg Mpc⁻³ yr⁻¹ lead to optical depth for p γ interactions as τ_{p γ} ≳ 0.4, assuming that muon neutrino emission rate follows a power law E^−Γ, where Γ ≃ 2.19 (IceCube Collaboration 2017). Using the connection between τ_{p γ} and τ_{γ γ} as (Murase et al. 2016; Halzen et al. 2019)

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }\approx \displaystyle \frac{{\eta }_{\gamma \gamma }{\sigma }_{\gamma \gamma }}{{\eta }_{p\gamma }{\hat{\sigma }}_{p\gamma }}{\tau }_{p\gamma },\end{eqnarray} \tag{ 2 }$$

τ_{p γ} ≳ 0.4 gives ${\tau }_{\gamma \gamma }\simeq { \mathcal O }(100)$, where the corresponding values of η_{γ γ} , η_{p γ} , σ_{γ γ} , and ${\hat{\sigma }}_{p\gamma }$ are given for the Δ-resonance and direct pion production by, e.g., Halzen et al. (2019). High-energy pionic γ-rays are then unable to escape from an emission region that effectively generates high-energy neutrinos.

Astrophysical (median) muon neutrino energies (E_ν ) between 119 and 4.8 PeV in the eight years of IceCube data (IceCube Collaboration 2017) correspond to proton energies (E_p ) from 2.4 PeV to 96 PeV (since E_p ≈ 20E_ν ). Assuming the meson production is dominated by the Δ-resonance, these protons interact with X-ray and UV target photons (see, e.g., Murase et al. 2018), and constrain the two-photon annihilation depth at E_γ ∼ 5–200 GeV energies. This energy range is well within the Fermi-LAT energy range of 0.1–300 GeV. Considering head-on collisions, the 5–200 GeV Fermi photons produce e⁺ e⁻ pairs with 10 keV(γ/10)²–3 eV(γ/10)² photons. In blazar environments the electron Lorentz factor is relativistic, rendering the energy of photons absorbing pionic γ-rays to the keV, MeV, or even GeV regime.

As inverse Compton scattering is Klein–Nishina suppressed at the highest energies (e.g., Moderski et al. 2005), we consider that e⁺ e⁻ pairs produced in the interaction of GeV pionic γ-rays with X-ray target photons losing their energy primarily via synchrotron cooling. The characteristic frequency of synchrotron emission (f_ν ) by Bethe–Heitler pairs is expected in the keV regime, while for photomeson production (including pairs from two-photon annihilation) f_ν is expected in the MeV regime (Murase et al. 2018). Above 10¹⁸ eV the Bethe–Heitler, above 5 × 10¹⁹ the photomeson pair production dominates, therefore f_ν shifts from keV to MeV as cosmic rays gain ultra-high energies and begin to produce neutrinos.

This scenario is consistent with multimessenger observations of the extragalactic, diffuse high-energy emission: in order to explain the diffuse high-energy neutrino flux together with the Fermi detection of the extragalactic γ-ray flux can only be done by either assuming very flat neutrino spectra or taking into account strong γ-ray absorption in the sources (Ahlers & Halzen 2015).

4.2. A Case Study: PKS 1502+106

Observations combined with the above scenario suggest that in mode 1 of PKS 1502+106 γ-rays originate in a transparent source with insufficient target photon density to produce a detectable neutrino flux (γ-transparent mode). The strong correlation between the γ and radio flux suggests the leptonic origin of γ-rays in this state. In mode 2 a dense photon (or even proton) target field interacting with the jet provides the conditions to produce neutrinos (γ-suppressed mode). The pionic γ-rays accompanying cosmic neutrinos lose energy in the source due to large τ_{γ γ} as shown above, and emerge at Earth below the Fermi-LAT detection threshold. We show the γ-transparent and the γ-suppressed blazar modes in Figure 2.

Zoom In Zoom Out Reset image size

The powerful radio flare started in the beginning of 2014 reveals a very active outflow. We calculated the γ to radio ratio in 56 time bins of the γ-suppressed mode. This ratio was the smallest in the time bin of the neutrino, yielding a p-value of 0.018. This suggests the suppression was most effective about the IceCube neutrino detection.

A high cascading contribution from the pionic γ-rays is expected in the γ-suppressed mode, leading to a high X-ray flux. We show the Swift-XRT 0.3–10 keV light curve of PKS 1502+106 ¹¹ in Figure 1 (Stroh & Falcone 2013). We see the 0.1–300 GeV γ-ray flux drops after a sharp peak at about MJD 57000 (Figure 1). The X-ray flux was peaking at about MJD 58100, with the γ-flux still above average, with no neutrino detected. The scenario suggests Swift observed the synchroton emission of Bethe–Heitler pairs these times. The increasing radio flux might reveal particle acceleration, with protons probably coaccelerating with the electrons/positrons. IceCube detected a neutrino while the γ-flux was maximally suppressed, and the 0.3–10 keV X-ray flux was about halved. The lowered X-ray flux might be due to that at sufficiently high energies photomeson production dominates the pair production instead of the Bethe–Heitler, and a large part of the synchrotron emission from pairs appears at MeV energies, outside the 0.3–10 keV range of Swift.

This is a scenario that can work as already discussed by Halzen & Kheirandish (2020), but to prove this in detail for PKS 1502+106, detailed theoretical modeling of the multiwavelength data is necessary, which is beyond the scope of this paper.

TXS 0506+056. IceCube recorded a 290 TeV muon neutrino in temporal and directional coincidence with a γ-flare of TXS 0506+056, the correlation being statistically significant at the 3σ level (IceCube Collaboration et al. 2018). We show the blazar's γ-flux measured by Fermi-LAT and the Major Atmospheric Gamma Imaging Cerenkov (MAGIC) Telescopes in the middle panel of Figure 3. The Fermi-LAT flux from TXS 0506+056 drops by about 30% from the week prior to the neutrino, and starts rising again within 2 weeks after the neutrino, in the middle of a large γ-ray flare. The MAGIC flux is undetectable at the time of the neutrino emission, but then rapidly increases by at least a factor of 10 within 2 weeks (IceCube Collaboration et al. 2018). Interestingly, a similar rise, although on the timescale of hours was also detected by the MASTER optical telescope in the aftermath of the neutrino detection from the blazar (Lipunov et al. 2020). We note MAGIC observed TXS 0506+056 for 2 hr on 2017 September 24 (MJD 58020) under nonoptimal weather conditions (IceCube Collaboration et al. 2018). Ansoldi et al. (2018) interpreted the γ-ray variability between MJD 58024 and MJD 58030 with a one-zone lepto-hadronic model, based on interactions of electrons and protons coaccelerated in a structured jet.

Zoom In Zoom Out Reset image size

PKS B1424-418. A temporal and directional coincidence was found between the outburst of blazar PKS B1424-418 and a PeV-energy neutrino event (HESE 35, dubbed as "Big Bird"; Kadler et al. 2016). In the first panel of Figure 3 we show the γ-ray light curve of PKS B1424-418 around the time of the neutrino detection (Kadler et al. 2016). While there is a longer outburst, at the actual time of the neutrino there is a sudden suppression in the γ-flux, which then quickly rises again afterwards. Therefore, we find that the neutrino detection coincided with a local γ minimum. Kadler et al. (2016) determined 5% for a chance coincidence. Our scenario suggests the neutrino should come in a γ-ray dip, which significantly lowers the probability of the chance coincidence.

For comparison, we also show the γ-ray light curve of PKS 1502+106 in the third panel of Figure 3, binned in 7 days wide bins. The temporary appearance of an X-ray/UV target photon field can temporarily suppress the γ-flux, even if the source is flaring. The timescale of the suppression potentially depends on the characteristics of target field, which is probably more massive in the case of PKS 1502+106. Investigating the origin of the target field source by source is beyond the scope of the paper.

We showed that the γ-ray flux of the three blazars to which a high-energy neutrino counterpart was associated, experienced a temporary γ-suppression at the time of the detected neutrino. If our scenario is correct, then it is substantially easier to discover neutrino–blazar connections due to the short time coincidence allowed between γ-suppressed periods and neutrinos. With this we suggest revising the strategy to find the origin of the cosmic high-energy neutrinos. It is possible that previously recorded IceCube and Fermi-LAT observations are already sufficient to identify the origin of the bulk of the high-energy neutrinos detected in the universe.

The emerging picture from our results is that efficient neutrino production is achieved by the combination of (i) an energetic outflow from the central supermassive black hole in the blazar and (ii) the temporary increase in the target photon density. The presence of the energetic outflows are indicated by the high radio flux and the high γ-flux before and after the detected neutrinos. Therefore, γ-suppression is likely not caused by changes in the outflow from the black hole, but by the appearance of target photons at greater radii. This is a temporary condition as the γ-flux rapidly increases afterwards.

The main message of this Letter is that the temporal γ-suppression potentially resolves the apparent contradiction of the blazar models simultaneously producing a detectable neutrino flux and a γ flare, since at the time of efficient neutrino production the observed γ-flux drops. More observations and a larger sample of sources behaving in a similar way to TXS 0506+056 and PKS 1502+106 would establish a more stable ground to this scenario, and we encourage studies in this direction.

The authors would like to thank Roger Clay and Marcos Santander for their valuable comments on the manuscript. They are grateful to the anonymous referee for the comments that improved the quality of this manuscript. E.K. was supported by the Premium Post-doctoral Research Program of the Hungarian Academy of Sciences. Support of NKFIH grant No. 123996 in the early stage of this work is acknowledged. I.B. acknowledges the support of the National Science Foundation under grant No. #1911796 and of the Alfred P. Sloan Foundation. J.B.T. acknowledges support from the DFG, project TJ 62/8-1. F.H.'s research is supported by the U.S. National Science Foundation under grants PLR-1600823 and PHY-1913607. On behalf of Project "fermi-agn" we are thankful for the usage of ELKH Cloud (formerly MTA Cloud, https://cloud.mta.hu/), which significantly helped us achieve the results published in this paper. This research has made use of data from the OVRO 40 m monitoring program (Richards et al. 2011) which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.