TeV Neutrinos and Hard X-Rays from Relativistic Reconnection in the Corona of NGC 1068

Active galactic nuclei (AGN) are among the most energetic sources in the Universe. They are powered by accreting supermassive black holes, and their main emission features (see, e.g., Marconi et al. 2004; Ghisellini 2013; Padovani et al. 2017) can be summarized as (i) an optical-ultraviolet (OUV) thermal component, often referred to as the "big blue bump," which is produced by the accretion disk; (ii) a second thermal component radiated by the dusty torus in the infrared; and (iii) an X-ray power law extending up to 100 keV, which is believed to be associated with the AGN "coronal" region.

The coronal region of AGN has been proposed as a production site of multimessenger emission since the late 1970s (Eichler 1979; Silberberg & Shapiro 1979; Berezinsky & Ginzburg 1981). Very recently, the IceCube collaboration has reported compelling evidence of a high-energy neutrino flux from NGC 1068 (Aartsen et al. 2020; Abbasi et al. 2022), a nearby Seyfert II galaxy located at a distance of ∼10.1 Mpc (Tully et al. 2008; Padovani et al. 2023, in preparation). This neutrino flux has been measured with a significance of 4.2σ in the energy band 1.5–15 TeV, showing a soft neutrino spectrum $d{{\rm{\Phi }}}_{\nu }/{{dE}}_{\nu }\propto {E}_{\nu }^{-3.2}$. A comparable TeV gamma-ray flux, which would be expected for an optically thin source, is not observed, with stringent upper limits set by the MAGIC telescope(Acciari et al. 2019). This suggests that the neutrino source must be located close to the AGN innermost regions, where the dense optical and X-ray radiation would reprocess gamma rays to lower energies (Berezinsky & Ginzburg 1981). Thus, the GeV gamma rays observed by Fermi-LAT (Abdo et al. 2010) are hard to explain within the coronal model and generally attributed to a different region, e.g., the weak jet observed in the core of the galaxy (Lenain et al. 2010), the circumnuclear starburst (Yoast-Hull et al. 2014; Ambrosone et al. 2021; Eichmann et al. 2022), a large-scale AGN-driven outflow (Lamastra et al. 2016), a combination of successful and failed line-driven wind (Inoue et al. 2022), or an ultrafast outflow (Peretti et al. 2023).

Neutrino emission from AGN cores directly points to the presence of a relativistic hadronic population, but the specific particle acceleration mechanism at work remains an open question. Inoue et al. (2020) and Murase et al. (2020) provided the first phenomenological modeling of NGC 1068, prescribing that protons are energized via diffusive shock acceleration or stochastic acceleration, respectively. Both models highlighted the importance of Bethe–Heitler (BH) interactions on optical photons as the main mechanism limiting proton acceleration. Kheirandish et al. (2021) discussed magnetic reconnection as a potential mechanism for proton acceleration in a weakly magnetized corona. Assuming a proton power-law spectrum with a slope of −2, they computed the neutrino emission from pp and p γ inelastic collisions and applied their analysis to NGC 1068. Recently, Murase (2022) provided model-independent multimessenger constraints on the acceleration processes, as well as on the production mechanism of high-energy neutrinos in AGN coronae. This study highlighted the importance of the electromagnetic cascade for reprocessing the GeV–TeV gamma rays into the MeV band.

The aim of this Letter is to explore the role of reconnection in magnetically dominated plasmas for powering the X-ray and neutrino emission in AGN coronae. In this scenario, rough energy equipartition between magnetic fields, X-ray radiation, and relativistic protons is established in the reconnection region. Motivated by recent 3D particle-in-cell (PIC) simulations of relativistic reconnection (Zhang et al. 2021, 2023; Chernoglazov et al. 2023), we assume that the spectrum of accelerated protons is a broken power law, with the break energy E_p,br being constrained by energy conservation (i.e., the energy density of accelerated protons is at most comparable to the magnetic energy density). ⁷ The proton spectrum is hard below the break (slope around −1), whereas it is softer above the break, which could explain the soft neutrino spectrum of NGC 1068. Our model has only two free parameters, namely, the size of the reconnection region and the number density of nonthermal protons in the coronal region, which directly determines E_p,br. All other parameters can be either inferred by observations or benchmarked with PIC simulations.

In this Letter, we demonstrate that the neutrino emission expected from relativistic reconnection in the coronal region of NGC 1068 is consistent with IceCube observations, and we show that the peak neutrino luminosity scales quadratically with the X-ray luminosity and is inversely proportional to the black hole mass. For large X-ray luminosities, protons cool so fast that they lose most of their energy inside the reconnection layer, which effectively becomes a calorimeter; in this regime, the scaling saturates into a linear proportionality with the X-ray luminosity, independently of the black hole mass and layer size. We find that BH and photohadronic interactions initiate an electromagnetic cascade that ultimately establishes and maintains a population of electron–positron pairs, which naturally makes the corona moderately optically thick for the Comptonization of low-energy photons up to the hard X-ray band.

In this section, we identify the main properties of the reconnection layer. Assuming that the hard X-ray flux comes from Comptonization within the reconnection layer, we can infer the electron number density and also provide an estimate for the magnetic field strength in the reconnection layer.

We envision the formation of a reconnection layer in the vicinity of a black hole with mass M = 10⁷ M₇ M_⊙. In order to maintain our conclusions as generally as possible, we do not specify the exact origin and location of this reconnection layer. One possible scenario is the equatorial current sheet that naturally arises in the innermost region of magnetically arrested disks (Avara et al. 2016; Dexter et al. 2020; Ripperda et al. 2020, 2022; Porth et al. 2021; Nathanail et al. 2022; Scepi et al. 2022). Another possibility is the reconnection layer adjacent to the jet boundary (Nathanail et al. 2020; Chashkina et al. 2021; Davelaar et al. 2023) or beyond the Y-point of field lines connected to both the disk and the black hole (El Mellah et al. 2022, 2023). Regardless of its specific origin, the presence of this macroscopic reconnection layer facilitates the rapid conversion of the available magnetic energy into particle energy. A fraction of this energy can be reprocessed into hard X-rays. Traditionally, this has been associated with the Comptonization of low-energy photons from the accretion disk by a population of hot coronal electrons with ∼100 keV temperature (Sunyaev & Truemper 1979; Sunyaev & Titarchuk1980; Haardt & Maraschi 1991, 1993). More recently, alternative scenarios have been proposed, in which X-rays are Comptonized by the bulk motions of a turbulent plasma (Groselj et al. 2023) or by the transrelativistic bulk motions of the stochastic chain of reconnection plasmoids (Beloborodov 2017; Sironi & Beloborodov 2020; Sridhar et al. 2021, 2023). We refrain from adopting a specific scenario to describe the X-ray coronal emission of NGC 1068 and instead use some general considerations applicable to accreting systems with magnetically powered hard X-ray emission.

We consider a reconnection layer of length L, which we scale with the black hole gravitational radius r_g = GM/c² ≃ 1.5 × 10¹² M₇ cm. The cross-sectional area of the "upstream" plasma flowing into the reconnection layer is estimated as A_in ∼ L². The characteristic thickness of the layer of reconnected plasma (the reconnection "downstream") is determined by the size of the largest plasmoids, which can be estimated as w ∼ β_rec L. Here, β_rec represents the plasma inflow speed into the reconnection layer normalized to the speed of light. Its value is also representative of the ratio between the reconnection electric field and the reconnecting magnetic field B. The choice of β_rec ∼ 0.1 (Comisso & Bhattacharjee 2016; Cassak et al. 2017) is appropriate for the collisionless relativistic regime in which the magnetization σ = B²/4π ρ c² exceeds unity, where ρ is the plasma mass density. As shown below, this is indeed our regime of interest.

In both X-ray binaries (Zdziarski 2000; Malzac et al. 2001; Reig et al. 2018) and luminous AGN (Wilkins & Gallo 2015; Petrucci et al. 2020; Tripathi & Dewangan 2022), the shape of the X-ray spectrum suggests that the optical depth for Comptonization is moderate (i.e., ∼0.1–10). Motivated by radiative PIC simulations (Sridhar et al. 2021, 2023), we use the condition τ_T = n_e σ_T w ∼ 0.5, where σ_T is the Thomson cross section and n_e is the cold pair number density (hereafter the electron density), to infer

$$\begin{eqnarray}&&{n}_{e}\simeq {10}^{12}\,{\beta }_{-1}^{-1}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-1}{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 1 }$$

where β₋₁ = β_rec/10⁻¹.

The strength B of the reconnecting field can be obtained by assuming that the coronal hard X-ray emission comes from dissipation of magnetic energy in reconnection, as argued in Beloborodov (2017). We assume that a fraction η_X ∼ 0.5 of the Poynting flux crossing the reconnection plane is dissipated into Comptonized X-rays (the rest stays in electromagnetic fields; see, e.g., Sironi et al. 2015), namely,

$$\begin{eqnarray}&&{F}_{{\rm{X}}}={\eta }_{{\rm{X}}}\displaystyle \frac{{\beta }_{\mathrm{rec}}c}{4\pi }{B}^{2}.\end{eqnarray} \tag{ 2 }$$

Noting that the X-ray flux is related to the intrinsic X-ray luminosity as L_X = 2F_X L², where the factor of 2 accounts for emission from both sides of the layer, the magnetic energy density is written as

$$\begin{eqnarray}&&{u}_{B}=\displaystyle \frac{{B}^{2}}{8\pi }\simeq 3\times {10}^{7}{L}_{{\rm{X}},43}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-2}{\beta }_{-1}^{-1}\,\mathrm{erg}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 3 }$$

where L_X = 10⁴³ L_X,43 erg s⁻¹ is the bolometric intrinsic X-ray luminosity of NGC 1068. This implies a magnetic field

$$\begin{eqnarray}&&B\simeq 3\times {10}^{4}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-1}{L}_{{\rm{X}},43}^{1/2}{\beta }_{-1}^{-1/2}\,{\rm{G}}.\end{eqnarray} \tag{ 4 }$$

Taking the pair number density from the optical depth argument presented above, we can estimate the magnetization of the system as

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{{B}^{2}}{4\pi {n}_{e}{m}_{e}{c}^{2}}\simeq 70\,{L}_{{\rm{X}},43}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-1},\end{eqnarray} \tag{ 5 }$$

where we have assumed that ions do not appreciably contribute to the mass density; namely, the proton number density is n_p ≪ n_e m_e /m_p . As we discuss in Section 3.4, this condition is supported by neutrino observations. As anticipated above, reconnection in a system with σ ≫ 1 is identified as "relativistic."

In this section, we first discuss the physics of proton acceleration in relativistic reconnection and the expected proton spectrum. We then describe the hierarchy of timescales (or, equivalently, of energy scales), accounting for proton acceleration, escape from the reconnection layer, and various cooling losses. We determine the steady-state proton spectrum and use it to predict the TeV neutrino flux.

3.1. Spectrum of Accelerated Protons

3.2. Proton Acceleration, Escape, and Cooling

Zhang et al. (2021, 2023) elucidated the physics of particle acceleration in relativistic 3D reconnection. They showed that particles gain most of their energy in the inflow (upstream) region while meandering between the two sides of the reconnection layer. Here, the timescale for proton acceleration t_acc ≃ E_p /β_rec eBc is of the order of

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}\simeq {10}^{-3}\displaystyle \frac{{E}_{p}}{25\,\mathrm{TeV}}\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}{L}_{{\rm{X}},43}^{-1/2}{\beta }_{-1}^{-1/2}\,{\rm{s}}.\end{eqnarray} \tag{ 6 }$$

Particles leave the region of active acceleration when they get captured by one of the flux ropes/plasmoids of reconnected plasma. At that point, they are no longer actively accelerated, and they either advect out of the system (on a timescale t_esc) or cool due to radiative losses. In the discussion that follows, we focus on the proton spectrum in the downstream region of reconnected plasma, assuming—based on PIC simulations (Werner & Uzdensky 2017; Zhang et al. 2021, 2023; Chernoglazov et al. 2023)—that the prior phase of active energization injects a broken power-law spectrum scaling as $\propto {E}_{p}^{-1}$ below the break and $\propto {E}_{p}^{-s}$ above the break. Protons escape/advect out of the reconnection region on a typical timescale

$$\begin{eqnarray}&&{t}_{\mathrm{esc}}\simeq \displaystyle \frac{L}{c}\simeq 250\,\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\,{\rm{s}}.\end{eqnarray} \tag{ 7 }$$

The dominant channel of proton energy loss is photohadronic (p γ) inelastic scattering on the X-ray photon field; we discuss in Appendix A additional cooling processes that turn out to be subdominant. To compute the p γ cooling rate, we parameterize the X-ray coronal emission as follows. The X-ray energy density is given by u_X ≃ L_X/L² c. In the optically thick regime, we expect this estimate to be revised by a corrective factor proportional to τ_T , but since τ_T ∼ 1, we neglect the exact numerical factors in what follows. Motivated by spectral modeling of NGC 1068 (Bauer et al. 2015; Marinucci et al. 2016), we adopt a power-law spectrum of photon index −2 between ${E}_{{\rm{X}},\min }=100\,\mathrm{eV}$ and ${E}_{{\rm{X}},\max }=100\,\mathrm{keV}$,

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{{\rm{X}}}}{{{dE}}_{{\rm{X}}}}=\displaystyle \frac{{L}_{{\rm{X}}}}{3\mathrm{log}(10){E}_{{\rm{X}}}^{2}{L}^{2}c}.\end{eqnarray} \tag{ 8 }$$

The frequency range covered by the coronal spectrum is consistent with a scenario in which photons are initially injected into the OUV band from the accretion disk, typically at energies even below 100 eV, and subsequently undergo Comptonization in the corona up to ∼100 keV. We emphasize, though, that the X-ray spectrum below about 1–2 keV is not motivated by observations; on the other hand, as we discuss in more detail in Section 3.3, a different value of ${E}_{{\rm{X}},\min }$ would not alter our results.

We assume an effective cross section that is constant in energy, $\hat{\sigma }={\sigma }_{p\gamma }{K}_{p\gamma }=70$ μb, for interaction energies above the center-of-mass energy threshold _r = 390 (in units of m_e c²), where σ_{p γ} is the p γ total cross section, and K_{p γ} is the proton inelasticity (Atoyan & Dermer 2003; Dermer & Atoyan 2003). Using the X-ray spectrum in Equation (8), we estimate the p γ energy-loss timescale as

$$\begin{eqnarray}&&\begin{array}{l}{t}_{p\gamma }=\displaystyle \frac{9\mathrm{log}(10){L}^{2}\ {\epsilon }_{r}\,{m}_{e}{c}^{2}}{4\hat{\sigma }{L}_{{\rm{X}}}({E}_{p}/{m}_{p}{c}^{2})}\\ =\,4.8\times {10}^{3}\,{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{2}\displaystyle \frac{25\,\mathrm{TeV}}{{E}_{p}}{L}_{{\rm{X}},43}^{-1}\,{\rm{s}}\end{array}\end{eqnarray} \tag{ 9 }$$

for proton energy ${E}_{p}\lt {E}_{p}^{* }$, which is defined as

$$\begin{eqnarray}&&{E}_{p}^{* }=\displaystyle \frac{{\epsilon }_{r}\,{m}_{p}{m}_{e}{c}^{4}}{2{E}_{{\rm{X}},\min }}\simeq 930\,\mathrm{TeV}.\end{eqnarray} \tag{ 10 }$$

In contrast, for ${E}_{p}\gt {E}_{p}^{* }$, the cooling timescale flattens to the asymptotic energy-independent value

$$\begin{eqnarray}&&\begin{array}{l}{t}_{p\gamma }=\displaystyle \frac{3\mathrm{log}(10){L}^{2}{\epsilon }_{r}\,{E}_{{\rm{X}},\min }}{\hat{\sigma }{L}_{{\rm{X}}}}\\ \,=\,87\,{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{2}{L}_{{\rm{X}},43}^{-1}\,{\rm{s}},\end{array}\end{eqnarray} \tag{ 11 }$$

since the scattering does not happen anymore at threshold for pion production and becomes dominated by multi-pion interactions with the whole X-ray spectrum.

Another relevant energy scale is obtained by equating the cooling timescale with the escape/advection timescale

$$\begin{eqnarray}&&{E}_{p,\mathrm{cool}}=490\ \mathrm{TeV}\ {L}_{{\rm{X}},43}^{-1}\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}.\end{eqnarray} \tag{ 12 }$$

Above E_p,cool, protons cool much faster than they can escape the reconnection layer and therefore lose most of their energy while inside of it.

The maximum energy to which protons can be accelerated is set by the equality between the acceleration and the dominant energy-loss timescale. At very high energies, this can either be synchrotron (estimated in Appendix A) or p γ, yielding (first and second term in the square bracket, respectively)

$$\begin{eqnarray}&&\begin{array}{l}{E}_{p,\mathrm{rad}}\simeq 380\ {\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{1/2}{L}_{{\rm{X}},43}^{-1/2}{\beta }_{-1}^{3/4}\\ \min \left[1,5.6{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{1/2}{L}_{{\rm{X}},43}^{-1/2}{\beta }_{-1}^{-1/4}\right]\ \mathrm{PeV}.\end{array}\end{eqnarray} \tag{ 13 }$$

Figure 1 collects the energy-dependent timescales we have introduced for a benchmark choice of parameters representative of the likely conditions in NGC 1068 (see Section 3.6). We also show the synchrotron energy-loss timescale and the BH energy-loss timescales due to X-ray and OUV photons separately, as discussed in Appendix A. We can directly grasp the main features discussed in the text: the acceleration timescale is much faster than all the other timescales up to hundreds of PeV. However, protons cool faster than the escape timescale due to p γ interactions, already very close to the break energy.

Zoom In Zoom Out Reset image size

3.3. Steady-state Proton Spectrum

Here we consider the proton spectrum in the downstream region of reconnected plasma, assuming that the prior phase of active acceleration in the upstream injects a broken power law scaling as $\propto {E}_{p}^{-1}$ below the break and $\propto {E}_{p}^{-s}$ above the break. When accounting for cooling losses, the steady-state proton spectrum will be composed of several power-law segments that depend on the hierarchy of the proton energy scales. For ${E}_{p,\mathrm{br}}\lt {E}_{p,\mathrm{cool}}\lt {E}_{p}^{* }\lt {E}_{p,\mathrm{rad}}$, as shown in Figure 1, we expect

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{p}}{{{dE}}_{p}}\propto \left\{\begin{array}{l}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-1},{E}_{p}\lt {E}_{p,\mathrm{br}}\\ {\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-s},{E}_{p,\mathrm{br}}\lt {E}_{p}\lt {E}_{p,\mathrm{cool}}\\ {\left(\displaystyle \frac{{E}_{p,\mathrm{cool}}}{{E}_{p,\mathrm{br}}}\right)}^{-s}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{cool}}}\right)}^{-s-1},{E}_{p,\mathrm{cool}}\lt {E}_{p}\lt {E}_{p}^{* }\\ \displaystyle \frac{{E}_{p,\mathrm{br}}^{s}{E}_{p,\mathrm{cool}}}{{({E}_{p}^{* })}^{s+1}}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p}^{* }}\right)}^{-s},{E}_{p}^{* }\lt {E}_{p}\lt {E}_{p,\mathrm{rad}}.\end{array}\right.\end{eqnarray} \tag{ 14 }$$

Notice that the break at ${E}_{p}^{* }$, which is determined by the lower energy of the X-ray target ${E}_{{\rm{X}},\min }$, is typically much above the peak of the proton spectrum and therefore does not impact the neutrino production close to the peak. Therefore, a lower value for ${E}_{{\rm{X}},\min }$ would not significantly affect our results.

The overall normalization of the proton spectrum is determined by the condition of near-equipartition between magnetic field energy density and proton energy density. Specifically, we normalize the spectrum by assuming that the proton energy density is u_p = η_p u_B . In the absence of proton cooling losses, PIC simulations indicate that ${\eta }_{p}\,={({u}_{p}/{u}_{B})}_{\mathrm{uncool}}\gtrsim 0.1$ (Comisso & Jiang 2023; French et al. 2023; Totorica et al. 2023). In what follows, we consider the special case where the protons that carry most of the energy (i.e., those at E_p,br) are marginally fast-cooling, namely, E_p,cool ≳ E_p,br, and use η_p = 0.1η_p,−1 as a typical value. In Appendix B, we generalize our results to the fast-cooling regime where E_p,cool ≪ E_p,br.

If we take the energy hierarchy ${E}_{p,\mathrm{br}}\sim {E}_{p,\mathrm{cool}}\ll {E}_{p}^{* }\ll {E}_{p,\mathrm{rad}}$ with a postbreak slope of s = 3 (motivated by the neutrino observations), we can compute the total proton number density,

$$\begin{eqnarray}&&\begin{array}{l}{n}_{p}\simeq \displaystyle \frac{2{\eta }_{p}{u}_{B}}{3{E}_{p,\mathrm{br}}}\left[\displaystyle \frac{1}{3}+\mathrm{log}\left(\displaystyle \frac{{E}_{p,\mathrm{br}}}{{E}_{p,\min }}\right)\right]\approx \displaystyle \frac{20{\eta }_{p}{u}_{B}}{3{E}_{p,\mathrm{br}}}\\ \simeq 5\times {10}^{5}\displaystyle \frac{{L}_{{\rm{X}},43}{\eta }_{p,-1}}{{\beta }_{-1}}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-2}\displaystyle \frac{25\,\mathrm{TeV}}{{E}_{p,\mathrm{br}}}\ {\mathrm{cm}}^{-3},\end{array}\end{eqnarray} \tag{ 15 }$$

where we have replaced the square bracket by its typical order-of-magnitude value of 10 for E_p,br ∼ 25 TeV (i.e., the typical value needed to reproduce the IceCube neutrino measurements). Equation (15) confirms our initial assumption that n_p m_p ≪ n_e m_e ; see Equation (1). Charge neutrality therefore requires the largest part of the lepton population in the coronal region to be made of electron–positron pairs.

3.4. Neutrino Emission

Photohadronic collisions inject charged pions, which subsequently decay to neutrinos. We estimate the production rate of neutrinos and antineutrinos of all flavors as

$$\begin{eqnarray}&&{E}_{\nu }^{2}{Q}_{\nu }={L}^{3}{\beta }_{\mathrm{rec}}{\left[{K}_{\nu }\displaystyle \frac{{E}_{p}^{2}{{dn}}_{p}/{{dE}}_{p}}{{t}_{p\gamma }}\right]}_{{E}_{p}=20{E}_{\nu }},\end{eqnarray} \tag{ 16 }$$

where Q_ν is the rate of neutrino production differential in neutrino energy, β_rec L³ is the volume of the reconnection region, and the factor K_ν is the fraction of proton energy converted to neutrinos. For interactions at Δ resonance, K_ν = 1/4, while for multi-pion interactions, K_ν = 1/2. For the target X-ray spectrum we are considering, using the two-bin parameterization for the p γ interaction rate of Dermer & Atoyan (2003), we find that around E_p ≃ E_p,br, a reasonable approximation is K_ν ≃ 0.35. We assume here that each neutrino carries on average an energy of E_ν = E_p /20.

Since ${t}_{p\gamma }\propto {E}_{p}^{-1}$ for ${E}_{p}\lt {E}_{p}^{* }$, the neutrino spectrum is hardened by one power of energy as compared to the proton spectrum, whereas it has the same slope as the proton spectrum for ${E}_{p}\gt {E}_{p}^{* }$. Taking the postbreak proton slope to be s = 3 as in Equation (14), we expect a broken power law with ${E}_{\nu }^{2}{Q}_{\nu }\propto {E}_{\nu }^{2}$ for E_ν < 0.05E_p,br; ${E}_{\nu }^{2}{Q}_{\nu }\propto {E}_{\nu }^{0}$ for 0.05E_p,br < E_ν < 0.05E_p,cool; and ${E}_{\nu }^{2}{Q}_{\nu }\propto {E}_{\nu }^{-1}$ for E_ν > 0.05E_p,cool (since the proton spectrum is steepened by photohadronic losses to ${{dn}}_{p}/{{dE}}_{p}\propto {E}_{p}^{-4}$; see also Equation (14)). Thus, for our benchmark value s = 3 of the postbreak proton slope, the neutrino spectrum is consistent with IceCube observations. The peak value of ${E}_{\nu }^{2}{Q}_{\nu }$, which lies in the energy range 0.05E_p,br–0.05E_p,cool, is

$$\begin{eqnarray}&&\begin{array}{l}{({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}=\displaystyle \frac{2{K}_{\nu }{\eta }_{p}{E}_{p,\mathrm{br}}{L}_{{\rm{X}}}^{2}\hat{\sigma }}{27\mathrm{log}(10){\eta }_{{\rm{X}}}{\epsilon }_{r}{m}_{e}{m}_{p}{c}^{5}L}\\ \simeq 5.9\times {10}^{39}{\eta }_{p,-1}{L}_{{\rm{X}},43}^{2}\displaystyle \frac{5{r}_{g}}{{{LM}}_{7}}\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\,\mathrm{TeV}}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{array}\end{eqnarray} \tag{ 17 }$$

For L_X,43 ≳ 5–10, this is comparable with the required luminosity needed to explain the IceCube flux, provided that E_p,br is in the range 50–100 TeV.

3.5. Proton-mediated Pair Enrichment

So far, we have left unspecified the origin of the leptonic population, whose number density we can infer from the Compton opacity. A thrilling possibility is that the proton population itself is responsible for maintaining a steady density of cold pairs in the coronal region. MeV photons, which are primarily emitted by BH relativistic pairs (as we will show in Section 3.6), as well as cascaded from p γ-produced gamma rays, inject pairs almost at rest after their attenuation. Because of the hadronic origin of MeV photons, we argue that the typical luminosity in the MeV band is comparable with the neutrino luminosity (see also Section 3.6). An order-of-magnitude estimate of the MeV photon energy density is

$$\begin{eqnarray}&&{u}_{\mathrm{MeV}}\sim \displaystyle \frac{{({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}}{{L}_{{\rm{X}}}}{u}_{{\rm{X}}}\simeq \displaystyle \frac{{({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}}{{L}^{2}c},\end{eqnarray} \tag{ 18 }$$

and their number density is n_MeV ∼ u_MeV/1 MeV. While this is only an order-of-magnitude estimate, it captures the expected scaling with the energy injected by hadrons; we validate it by comparing with our numerical simulations in Section 3.6. Finally, we can estimate the number of pairs produced by MeV photon–photon attenuation as

$$\begin{eqnarray}&&\begin{array}{l}{n}_{e}^{\gamma \gamma }\sim 0.1{\sigma }_{T}{n}_{\mathrm{MeV}}^{2}L\\ \simeq 2.5\times {10}^{6}\,{L}_{{\rm{X}},43}^{4}{\eta }_{p,-1}^{2}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-5}{\left(\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\,\mathrm{TeV}}\right)}^{2}\,{\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ 19 }$$

For L_X,43 ≳ 5–10, the attenuation of ∼MeV photons can therefore lead to a pair population that constitutes a sizable fraction of the pair density inferred from the Compton opacity in Equation (1). The above estimate should be taken as a lower limit, because we have only considered pair production via attenuation of ∼MeV photons. We have not included, e.g., the contribution of relativistic pairs—produced in photohadronic interactions and from the attenuation of GeV–TeV gamma-ray photons—that cool down to Lorentz factors γ_e ∼ 1.

3.6. Neutrino Emission from NGC 1068

We propose that reconnection-powered coronal neutrino emission is the dominant source of the neutrino flux observed by IceCube in coincidence with NGC 1068. Results are shown for a luminosity distance d_L = 10.1 Mpc (see Tully et al. 2008; Padovani et al. 2023, in preparation) and a black hole mass M₇ = 0.67 (see Greenhill et al. 1996; Padovani et al. 2023, in preparation).

The first requirement is that the neutrino luminosity matches the observed all-flavor luminosity observed by IceCube. Assuming a peak neutrino energy of 1.25 TeV, corresponding to a proton break energy E_p,br ≃ 25 TeV, the all-flavor neutrino luminosity is (0.8–4) × 10⁴² erg s⁻¹ (based on the best-fit value and 95% confidence region shown in Figure 4 of Abbasi et al. 2022 after accounting for the difference in the luminosity distance used). The second requirement is that the neutrino spectrum observed by IceCube must be sufficiently soft, which can be accommodated in our model if the postbreak proton slope is s ∼ 3 and, simultaneously, the cooling break E_p,cool is not much larger than the break energy E_p,br ≃ 25 TeV. In this way, the neutrino spectrum is not hardened significantly by the energy dependence of the p γ efficiency (${t}_{p\gamma }\propto {E}_{p}^{-1}$ in Equation (16)), and the reconnection layer acts almost as a proton calorimeter. Therefore, our second constraint is that E_p,cool ≳ E_p,br ≃ 25 TeV. We take E_p,cool to lie between 25 and 100 TeV (the fast-cooling case E_p,cool ≪ E_p,br is discussed in more detail in Appendix B.)

Figure 2 shows the regions of parameter space where both requirements are simultaneously satisfied according to our analytical estimates. We consider as independent parameters the size L of the reconnection layer and the X-ray luminosity in the 2–10 keV range; while the latter is, in principle, observed, there is a relatively large uncertainty on its value, ${L}_{{\rm{X}},[2\mbox{--}10]\,\mathrm{keV}}\,={3}_{-2}^{+3}\times {10}^{43}$ erg s⁻¹ (see Marinucci et al. 2016 and Padovani et al. 2023, in preparation; we rescaled their luminosity to account for the different choice of luminosity distance), motivating the range chosen for the figure. The axis at the top shows the bolometric (0.1–100 keV) intrinsic X-ray luminosity L_X, which we have used in our analytical estimates. We also show in green the region where the pairs produced by the attenuation of MeV photons according to Equation (19) can contribute a sizable fraction, if not the entirety, of the pair population inferred from the Compton opacity argument. The allowed parameter space provides a coherent explanation for both the neutrino signal and the moderate optical depth of the coronal regions.

Zoom In Zoom Out Reset image size

To confirm our analytical estimates, we have also performed numerical calculations of the neutrino and gamma-ray production within a spherical coronal region of size R; we describe in Appendices C and D the relation between the spherical geometry of our numerical calculations and the planar one used in our analytical estimates, as well as the main properties of the numerical code we used. We show the resulting neutrino spectral energy distribution in Figure 3 (orange solid line), whose peak is in reasonable agreement with the analytical expectation (orange circle). ⁸ Our calculations do not account for neutron–photon interactions; we have separately evaluated them numerically, finding that they can account for a 50% increase in the neutrino flux shown in Figure 3.

Zoom In Zoom Out Reset image size

The large opacity to photon–photon pair production of the coronal region suppresses the gamma-ray spectrum at energies ≳1 MeV, making it consistent with the observations. Synchrotron radiation of BH pairs is an important source of ∼MeV photons, which produce pairs with γ_e ∼ 1. Attenuation of more energetic photons arising from p γ interactions leads to the production of relativistic secondary electrons and positrons, which cool down to γ_e ∼ 1 within a dynamical time due to the strong synchrotron losses. As a result, attenuation of gamma-ray photons produced from photohadronic interactions provides another channel for injecting cold pairs into the system. With our numerical code, we find that the density of pairs resulting from the attenuation of hadronically produced photons is n_e,cold ∼ 3 × 10¹¹ cm⁻³, which leads to τ_T = σ_T R n_e,cold ∼ 0.3. This numerical value is also consistent with the analytical estimate computed for a spherical geometry; see Equation(C11).

We note that the radiative calculations we presented have some limitations. First, cold pairs are "passive" (i.e., nonemitting) particles in our code, since the radiative processes of nonrelativistic electrons (e.g., cyclotron radiation and Compton scattering) are not implemented. Second, we do not account for the energization of leptons due to magnetic reconnection or their synchrotron emission, which might be important for the low-energy end of the photon spectrum. For this reason, we only show the photon spectrum down to 1 eV and consider the X-ray coronal emission given (i.e., no interaction between cold pairs and X-ray photons is considered).

We propose that protons energized by magnetic reconnection in a highly magnetized, pair-dominated plasma interact with hard X-rays from the corona and produce the high-energy neutrinos observed by IceCube from the Seyfert galaxy NGC 1068. The key ingredient of this study—grounded on recent 3D PIC simulations of relativistic reconnection—is that the X-ray coronal emission, the energy density of the reconnecting magnetic fields, and the relativistic proton population are all connected with each other.

If hard X-ray photons are produced via Comptonization of lower-energy photons by pairs in the reconnection region, as previously proposed, then the energy density of up-scattered photons will be a fraction (of order ∼0.1) of the energy density of the reconnecting magnetic fields. Protons accelerated in the reconnection region will also carry a fraction of ∼0.1 of the magnetic energy density. To avoid an "energy crisis," the hard $\propto {E}_{p}^{-1}$ proton spectrum established at low energies must soften ($\propto {E}_{p}^{-s}$, with s > 2) above a characteristic break energy E_p,br. From the peak of the observed neutrino flux, we infer E_p,br ≃ 25 TeV. Under these conditions, the neutrino emission can be fully determined by only two parameters: the size of the reconnection region and the X-ray luminosity of the corona (which is observationally constrained but with large error bars).

In our scenario, the proton break energy happens to be just below the energy at which photohadronic cooling becomes faster than the proton escape/advection from the reconnection layer. This implies that most of the energy in the accelerated protons is lost to neutrinos and gamma rays. While the latter are not able to escape the optically thick reconnection layer and are cascaded down to lower energies, neutrinos carry away a significant fraction of the proton energy, explaining the relatively large flux observed by IceCube. Our scenario naturally leads to a robust scaling of the peak neutrino luminosity with black hole mass and X-ray luminosity, i.e., ${({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}\propto {M}^{-1}{L}_{{\rm{X}}}^{2}$, which comes from the direct connection between X-ray, magnetic, and nonthermal proton energy density. In the case of fast cooling, with E_p,cool ≪ E_p,br, the scaling with L_X becomes linear with no dependence on black hole mass and layer size, since the reconnection layer becomes calorimetric.

Intriguingly, we also find that the pair population induced by the cascade initiated by p γ interactions and the attenuation of BH-produced MeV photons is comparable with the one inferred from the Compton optical thickness of the corona. Even though this result is not a strict requirement, i.e., the bulk of the pairs might have a different origin, it is a natural outcome of our model, at least for the parameters of NGC 1068. We are therefore prompted to speculate that the nonthermal proton population may also play an important role in maintaining the large lepton density in AGN coronae. A more thorough dynamic study of the interplay between protons and X-ray coronal photons is yet to be carried out.

Based on our findings, reconnection layers in AGN coronae are magnetically dominated (i.e., the magnetic energy density is much larger than the plasma rest-mass energy density) and significantly pair-enriched, to the point that the upstream mass density is mostly contributed by electron–positron pairs. Equivalently, if the magnetization σ is normalized to the pair mass density, we have 1 ≪ σ ≪ σ_p , where the "proton magnetization" is (see Equation (15))

$$\begin{eqnarray}&&{\sigma }_{p}=\displaystyle \frac{{B}^{2}}{4\pi {n}_{p}{m}_{p}{c}^{2}}\sim \displaystyle \frac{3{E}_{p,\mathrm{br}}}{10{\eta }_{p}{m}_{p}{c}^{2}}\simeq 8\times {10}^{4}\,{\eta }_{p,-1}^{-1}\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\,\mathrm{TeV}}.\end{eqnarray} \tag{ 20 }$$

The low value of the proton number density (equivalently, the high σ_p ) seems to be somewhat needed for the observability of the neutrino signal. A much higher proton density would correspond to a much lower break energy (at fixed magnetic energy density), so the neutrino signal would peak in an energy range unobservable by IceCube.

Interestingly, the three scenarios presented in Section 2 for the formation of a large-scale coronal reconnection layer are likely consistent with a pair-dominated, proton-poor composition (in GRMHD simulations, the particle density in the layer is found to be dominated by artificially injected density floors, rather than the physical density of the initialized electron–proton torus). ⁹ To the best of our knowledge, 3D PIC simulations in this regime (1 ≪ σ ≪ σ_p ) have yet to properly quantify how the postbreak proton slope s depends on the guide field strength. Werner & Uzdensky (2017) demonstrated that the high-energy spectral slope is steeper for stronger guide fields—with s ≃ 3 for a guide field comparable to the reconnecting field—but their simulations adopted an electron–positron composition. If our observationally motivated choice of s = 3 is indeed to be attributed to a guide field comparable to the reconnecting field, then the reconnection layer at the jet boundary is the most likely scenario (e.g., Sironi et al. 2021). But in general, whether a relatively strong guide field is attained in the three scenarios presented in Section 2 is yet to be robustly assessed. As far as current simulations are concerned (Ripperda et al. 2022), equatorial current sheets do not seem to host a guide field, although this statement has been verified only for a relatively large black hole spin. Finally, for the Y-point of field lines connected to the disk and the black hole, El Mellah et al. (2022) exhibits a low guide field, but the field geometry here is constrained by the boundary conditions to be anchored to the accretion disk. Future simulations, where the disk is not merely a boundary condition but is self-consistently evolved, will be needed to study the properties of the reconnection layer in this case. Such studies will be essential to support our observationally motivated choice of s ≃ 3 for NGC 1068.

A direct extrapolation of our results for NGC 1068 to other AGN might suggest that neutrino emission from coronal regions requires the presence of a jet. However, the requirement of a significant guide field, which could be realized in a jet/disk boundary layer, only holds for NGC 1068, for which a soft neutrino spectrum is observed. Harder neutrino spectra can be produced in equatorial current sheets at AGN without jets. In conclusion, steep neutrino spectra (s > 2) from coronal regions of AGN would typically require reconnection layers in the jet/disk boundary where stronger guide fields may be found; otherwise, a jet is not needed. Future discoveries of other neutrino-emitting Seyfert galaxies with IceCube and KM3Net will critically test our proposed scenario.

Note added.—While we were finalizing this work, the paper by Mbarek et al. (2023) appeared. The two works are independent and complementary. Some arguments are similar—most notably, by assuming that Comptonization happens in a magnetically dominated corona, one can infer the coronal density and field strength directly from X-ray observations—but there are also important differences. Mbarek et al. (2023) invoked acceleration in reconnection as the injection stage for further acceleration by turbulence. Instead, we propose a simpler, better-constrained model based only on reconnection. We also supplement our analytical estimates with numerical calculations of the full electromagnetic cascade, highlighting the importance of reconnection-accelerated protons for sustaining the population of electron–positron pairs that makes the corona moderately Compton thick.

This work was supported by a grant from the Simons Foundation (00001470; to L.S.). We are grateful to Paolo Padovani for useful comments on the draft. L.S. acknowledges support from DoE Early Career Award DE-SC0023015. This research was facilitated by the Multimessenger Plasma Physics Center (MPPC), NSF grant PHY-2206609, to L.S. M.P. acknowledges support from the MERAC Fondation through the project THRILL and the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the "2nd call for H.F.R.I. Research Projects to support Faculty members and Researchers" through the project UNTRAPHOB. D.F.G.F. is supported by the Villum Fonden under project No. 29388 and the European Union's Horizon 2020 Research and Innovation Program under Marie Skłodowska-Curie grant agreement No. 847523 "INTERACTIONS." L.C. acknowledges support from NSF grant PHY-2308944 and NASA ATP grant 80NSSC22K0667. E.P. acknowledges support from the Villum Fonden (No. 18994) and the European Union's Horizon 2020 research and innovation program under Marie Sklodowska-Curie grant agreement No. 847523 INTERACTIONS'. E.P. was also supported by Agence Nationale de la Recherche (grant ANR-21-CE31-0028).

In the main text, we have considered as a dominant loss channel the photohadronic interaction of reconnection-accelerated protons onto the X-ray photon field. Here we consider additional channels of energy loss.

Other works on the coronal emission, e.g., Murase et al. (2020), have pointed at BH interactions with the OUV photon field as the dominant energy-loss channel. The reason for our different conclusion is the location of the neutrino production site, which, in our scenario, is a reconnection region formed much closer to the central black hole. While the X-ray field is directly produced within the reconnection region, the OUV field is thought to be produced in a thin accretion disk lying at larger radii, ∼10–50 r_g ; therefore, the energy density of optical photons is geometrically diluted within the region of interest. In this section, we explicitly verify these statements.

Let us first consider BH scattering on the X-ray target to show that it is generally subdominant as compared to p γ interactions. The BH energy-loss timescale of a proton with Lorentz factor γ_p due to interactions with an isotropic photon field of differential number density distribution ${{dn}}_{\gamma }({E}_{\gamma })/{{dE}}_{\gamma }$ reads

$$\begin{eqnarray}&&{t}_{\mathrm{BH}}^{-1}({\gamma }_{p})={\sigma }_{T}\left(\displaystyle \frac{{m}_{e}}{{m}_{p}}\right)\displaystyle \frac{{m}_{e}{c}^{3}}{{\gamma }_{p}^{2}}{\int }_{{\gamma }_{p}^{-1}}^{+\infty }\displaystyle \frac{d\epsilon }{{\epsilon }^{2}}\displaystyle \frac{{{dn}}_{\gamma }(\epsilon {m}_{e}{c}^{2})}{{{dE}}_{\gamma }}{\rm{\Phi }}(2{\gamma }_{p}\epsilon ),\end{eqnarray} \tag{ A1 }$$

where

$$\begin{eqnarray}&&{\rm{\Phi }}(x)={\int }_{2}^{x}d\overline{\epsilon }\,\overline{\epsilon }\displaystyle \frac{{\sigma }_{\phi e}(\overline{\epsilon })}{{\sigma }_{T}\sqrt{1+2\overline{\epsilon }}},\end{eqnarray} \tag{ A2 }$$

and we take the BH cross section ${\sigma }_{\phi e}(\overline{\epsilon })$ from the numerical fits of Stepney & Guilbert (1983). For the X-ray photon target density given by Equation (8), which we report here for convenience,

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{{\rm{X}}}}{{{dE}}_{{\rm{X}}}}=\displaystyle \frac{{L}_{{\rm{X}}}}{3\mathrm{log}(10){E}_{{\rm{X}}}^{2}{L}^{2}c},{E}_{{\rm{X}},\min }\lt {E}_{{\rm{X}}}\lt {E}_{{\rm{X}},\max },\end{eqnarray} \tag{ A3 }$$

we obtain the energy-loss rate as

$$\begin{eqnarray}&&{t}_{\mathrm{BH},{\rm{X}}}^{-1}({\gamma }_{p})=\displaystyle \frac{{\sigma }_{T}{L}_{{\rm{X}}}}{3\mathrm{log}(10){m}_{p}{L}^{2}{c}^{2}}\displaystyle \frac{{\xi }_{{\rm{X}}}({\gamma }_{p})}{{\gamma }_{p}^{2}},\end{eqnarray} \tag{ A4 }$$

where ξ_X(γ_p ) is the dimensionless function defined by

$$\begin{eqnarray}&&{\xi }_{{\rm{X}}}({\gamma }_{p})={\int }_{\min \left({\gamma }_{p}^{-1},\displaystyle \frac{{E}_{{\rm{X}},\min }}{{m}_{e}{c}^{2}}\right)}^{\tfrac{{E}_{{\rm{X}},\max }}{{m}_{e}{c}^{2}}}\displaystyle \frac{d\epsilon }{{\epsilon }^{4}}{\rm{\Phi }}(2{\gamma }_{p}\epsilon ).\end{eqnarray} \tag{ A5 }$$

For comparison, the photohadronic timescale for ${\gamma }_{p}\lt {\gamma }_{p}^{* }\,={E}_{p}^{* }/{m}_{p}{c}^{2}$ is

$$\begin{eqnarray}&&{t}_{p\gamma }^{-1}({\gamma }_{p})=\displaystyle \frac{4\hat{\sigma }{L}_{{\rm{X}}}}{9\mathrm{log}(10){L}^{2}{\epsilon }_{r}{m}_{e}{c}^{2}}{\gamma }_{p}.\end{eqnarray} \tag{ A6 }$$

The ratio between these two timescales does not depend on the X-ray luminosity or the size of the region; noting that the function ${\xi }_{{\rm{X}}}({\gamma }_{p})/{\gamma }_{p}^{3}$ has a maximum value of about 2 × 10⁻⁴, we easily verify that the ratio between the two timescales is bounded by

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{p\gamma }}{{t}_{\mathrm{BH},{\rm{X}}}}\lesssim 0.3.\end{eqnarray} \tag{ A7 }$$

Therefore, BH energy losses from the X-ray field are at most comparable and slightly subdominant with respect to the photohadronic interactions off the same target field for protons with ${E}_{p}\leqslant {E}_{p}^{* }$. Nevertheless, the energy lost via p γ interactions is divided into three particle species (neutrinos, gamma rays, and pairs), while the proton energy lost via BH is channeled to pairs alone. Therefore, the emission of BH pairs and neutrinos may be comparable, as shown in Figure 3 (see also Petropoulou & Mastichiadis 2015).

The additional target field in the region of interest is the OUV field, which is produced mostly as a multitemperature spectrum from the geometrically thin, optically thick accretion disk at larger radii. First of all, we provide an estimate of how far away from the central black hole the OUV is mostly produced. We consider an OUV-integrated luminosity L_OUV = ϕ L_X, where the dimensionless factor is typically ϕ ≃ 10 (Murase et al. 2020, and references therein). Assuming an accretion efficiency η_a ≃ 0.1, this points to a mass accretion rate

$$\begin{eqnarray}&&\dot{M}\simeq \displaystyle \frac{\phi {L}_{{\rm{X}}}}{{\eta }_{\mathrm{acc}}{c}^{2}}.\end{eqnarray} \tag{ A8 }$$

The effective temperature of the produced OUV field as a function of radius is determined by the standard Shakura–Sunyaev theory of geometrically thin accretion disks. The typical energy of the OUV field relevant for BH interactions is of order of 10 eV, so we consider an effective temperature T_OUV = 3.33 eV; the corresponding radius is

$$\begin{eqnarray}&&{r}_{\mathrm{OUV}}={\left(\displaystyle \frac{3{GM}\dot{M}}{8\pi {\sigma }_{\mathrm{SB}}{T}_{\mathrm{OUV}}^{4}}\right)}^{1/3}\simeq 75{r}_{g}{L}_{{\rm{X}},43}^{1/3}{M}_{7}^{-2/3},\end{eqnarray} \tag{ A9 }$$

where σ_SB is the Stefan–Boltzmann constant. This distance is typically much larger than the one considered for the reconnection region in our scenario, which is already suggestive of the subdominant nature of the OUV field in our region of interest. To prove this conclusively, we explicitly estimate the BH losses on the OUV field.

Since we will conclude that energy losses on this target are not an important energy-loss channel, the precise spectral shape is not a decisive feature for our purposes. We model the OUV field as a thermal distribution with a typical temperature T_OUV = 3.33 eV,

$$\begin{eqnarray}&&{n}_{\mathrm{OUV}}({E}_{\gamma })=\displaystyle \frac{15\phi {L}_{{\rm{X}}}}{4{\pi }^{5}{r}_{\mathrm{OUV}}^{2}{T}_{\mathrm{OUV}}^{4}}\displaystyle \frac{{E}_{\gamma }^{2}}{{e}^{{E}_{\gamma }/{T}_{\mathrm{OUV}}}-1},\end{eqnarray} \tag{ A10 }$$

which is valid in the limit r_OUV ≫ L. We can now estimate the BH energy-loss timescale off the OUV field as

$$\begin{eqnarray}&&{t}_{\mathrm{BH},\mathrm{OUV}}^{-1}({\gamma }_{p})=\displaystyle \frac{15{\sigma }_{T}\phi {L}_{{\rm{X}}}{m}_{e}^{4}}{4{\pi }^{5}{r}_{\mathrm{OUV}}^{2}{T}_{\mathrm{OUV}}^{4}{m}_{p}}\displaystyle \frac{{\xi }_{\mathrm{OUV}}({\gamma }_{p})}{{\gamma }_{p}^{2}}.\end{eqnarray} \tag{ A11 }$$

We can again determine the ratio between this energy-loss timescale and the photohadronic timescale, noting that the function ${\xi }_{\mathrm{OUV}}({\gamma }_{p})/{\gamma }_{p}^{3}$ has a maximum value of about 9 × 10⁻²⁵:

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{p\gamma }}{{t}_{\mathrm{BH},\mathrm{OUV}}}\lesssim 3\times {10}^{-3}\,{(L/5{r}_{{\rm{g}}})}^{2}{M}_{7}^{4/3}{L}_{{\rm{X}},43}^{-2/3}.\end{eqnarray} \tag{ A12 }$$

We see that for the typical values considered in this work, BH energy losses on the OUV photons are indeed negligible compared to the p γ energy losses. The geometric reduction of the OUV signal inside the corona may raise tensions with the general explanation for the coronal X-rays as coming from the Comptonization of OUV photons. However, in principle, X-rays could also originate from the Comptonization of the synchrotron photons from the accelerated electrons or the low-energy photons injected into the photohadronic cascade that we have examined in the main text.

An additional source of cooling is synchrotron in the strong magnetic field. The typical energy-loss timescale is

$$\begin{eqnarray}&&\begin{array}{l}{t}_{\mathrm{syn}}=\displaystyle \frac{3{m}_{p}^{4}{c}^{3}}{4{\sigma }_{T}{u}_{B}{m}_{e}^{2}{E}_{p}}\\ \simeq 2.4\times {10}^{5}\,{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{2}\displaystyle \frac{25\,\mathrm{TeV}}{{E}_{p}}{L}_{{\rm{X}},43}^{-1}{\beta }_{-1}\,{\rm{s}},\end{array}\end{eqnarray} \tag{ A13 }$$

which is again much longer that the t_{p γ} loss timescale.

In the main text, we have restricted ourselves to the regime where E_p,cool ≳ E_p,br, so that the protons carrying the most energy are not cooled significantly before escaping the layer. Here, we generalize our treatment to the fast-cooling case as well. For clarity, we denote by ${\eta }_{p}^{* }$ the fraction of magnetic energy density that would go into nonthermal protons in the absence of cooling. Thus, in the absence of cooling effects on spectral slopes, we would have

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{p}^{(\mathrm{uncool})}}{{{dE}}_{p}}=\displaystyle \frac{{\eta }_{p}^{* }{u}_{B}}{2{E}_{p,\mathrm{br}}^{2}}\left\{\begin{array}{l}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-1},{E}_{p}\lt {E}_{p,\mathrm{br}}\\ {\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-s},{E}_{p,\mathrm{br}}\lt {E}_{p}\lt {E}_{p,\mathrm{rad}}\end{array}\right.,\end{eqnarray} \tag{ B1 }$$

where s ≃ 3 as in the main text. Notice that, as compared to Equations (14) and (15) that refer to the marginally fast-cooling case with E_p,cool ≃ E_p,br, the numerical factor is different here, since the spectrum after the break decreases as ${E}_{p}^{-3}$ rather than ${E}_{p}^{-4}$. Therefore, in the marginally fast-cooling case with E_p,cool ≃ E_p,br, the relation between the fraction η_p we used there and the fraction ${\eta }_{p}^{* }$ defined here is ${\eta }_{p}=3{\eta }_{p}^{* }/4$.

Accounting for cooling, the steady-state spectrum can be approximated as

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{p}}{{{dE}}_{p}}=\displaystyle \frac{{{dn}}_{p}^{(\mathrm{uncool})}}{{{dE}}_{p}}\min \left[1,\displaystyle \frac{{t}_{p\gamma }}{{t}_{\mathrm{esc}}}\right].\end{eqnarray} \tag{ B2 }$$

Therefore, the spectral shape of the protons depends on the hierarchy between E_p,br and E_p,cool. The case E_p,cool ≳ E_p,br was already discussed in the main text. In the case E_p,cool < E_p,br, protons cool significantly even below the break energy, leading to a typical spectrum

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{p}}{{{dE}}_{p}}\propto \left\{\begin{array}{l}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-1},{E}_{p}\lt {E}_{p,\mathrm{cool}}\\ {\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-2},{E}_{p,\mathrm{cool}}\lt {E}_{p}\lt {E}_{p,\mathrm{br}}\\ {\left(\displaystyle \frac{{E}_{p,\mathrm{cool}}}{{E}_{p,\mathrm{br}}}\right)}^{-2}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{cool}}}\right)}^{-s-1},{E}_{p,\mathrm{cool}}\lt {E}_{p}\lt {E}_{p}^{* }\\ \displaystyle \frac{{E}_{p,\mathrm{br}}^{2}{E}_{p,\mathrm{cool}}^{s-1}}{{({E}_{p}^{* })}^{s+1}}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p}^{* }}\right)}^{-s},{E}_{p}^{* }\lt {E}_{p}\lt {E}_{p,\mathrm{rad}}.\end{array}\right.\end{eqnarray} \tag{ B3 }$$

Therefore, the peak neutrino spectrum, which lies in the range 0.05E_p,cool < E_p < 0.05E_p,br, is

$$\begin{eqnarray}&&{({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}=8.8\times {10}^{40}\,{\eta }_{p,-1}^{* }{L}_{{\rm{X}},43}\mathrm{erg}\,{{\rm{s}}}^{-1};\end{eqnarray} \tag{ B4 }$$

namely, the neutrino luminosity is a fixed fraction of the X-ray luminosity, independent of the size of the reconnection region. This happens because the region becomes fully calorimetric, converting all of the proton energy via p γ interactions.

The hadronically produced pair number density would also be affected by the different form of the proton spectrum. Using the proton distribution accounting for cooling, we find

$$\begin{eqnarray}&&\begin{array}{l}{n}_{e}^{\gamma \gamma }=1.9\times {10}^{6}\ {L}_{{\rm{X}},43}^{4}{\eta }_{p,-1}^{* 2}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-5}{\left(\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\,\mathrm{TeV}}\right)}^{2}\ {\mathrm{cm}}^{-3}\\ \min \left[1,{\left(\displaystyle \frac{{E}_{p,\mathrm{cool}}}{{E}_{p,\mathrm{br}}}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ B5 }$$

In the fast-cooling regime E_p,cool < E_p,br, we obtain

$$\begin{eqnarray}&&{n}_{e}^{\gamma \gamma }=7.2\times {10}^{8}{L}_{{\rm{X}},43}^{2}{\eta }_{p,-1}^{* 2}{\left(\displaystyle \frac{{{LM}}_{7}}{5{r}_{g}}\right)}^{-3}\,{\mathrm{cm}}^{-3}.\end{eqnarray} \tag{ B6 }$$

These results allow us to extend the regions in parameter space shown by Figure 2 to the regime with E_p,cool ≲ E_p,br ≃ 25 TeV. We show the corresponding regions in Figure 4. For large X-ray luminosities, to the right of the thin black dashed line, we enter the fast-cooling regime, where the peak neutrino luminosity scales in proportion to the X-ray luminosity and becomes independent of the size of the reconnection region. This shows that the IceCube neutrino luminosity can also be explained in this fast-cooling regime, and that in fact, an even larger number of pairs can be attained in this regime compared to the marginally fast cooling.

Zoom In Zoom Out Reset image size

The numerical code we use for the radiative calculations (see Appendix D) is set up for a spherical geometry. Here we discuss the main changes induced by considering a spherical, rather than plane, region. To clearly differentiate between the two cases, we will define the size of the spherical region as R​​​​​. In relating the planar to the spherical case, we will choose an effective size for the spherical region by equating the volumes for the sphere 4π R³/3 and the reconnection region β_rec L³. Therefore, R is related to the full length of the current sheet as $R=L{(3{\beta }_{\mathrm{rec}}/(4\pi ))}^{1/3}\simeq 0.28\ L{\beta }_{-1}^{1/3}$.

The energy density of X-rays is now connected to the luminosity by

$$\begin{eqnarray}&&{u}_{{\rm{X}}}=\displaystyle \frac{3{L}_{{\rm{X}}}}{4\pi {R}^{2}c}.\end{eqnarray} \tag{ C1 }$$

From the equality

$$\begin{eqnarray}&&{F}_{{\rm{X}}}=\displaystyle \frac{{L}_{{\rm{X}}}}{4\pi {R}^{2}}={\eta }_{{\rm{X}}}\displaystyle \frac{c}{4\pi }{\beta }_{\mathrm{rec}}{B}^{2},\end{eqnarray} \tag{ C2 }$$

we obtain the magnetic energy density

$$\begin{eqnarray}&&\begin{array}{l}{u}_{B}=\displaystyle \frac{{u}_{{\rm{X}}}}{6{\eta }_{{\rm{X}}}{\beta }_{\mathrm{rec}}}\\ \simeq 4.8\times {10}^{6}{L}_{{\rm{X}},43}{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\right)}^{-2}{\beta }_{-1}^{-1}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ C3 }$$

The corresponding magnetic field is

$$\begin{eqnarray}&&B\simeq 1.2\times {10}^{4}{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\right)}^{-1}{L}_{{\rm{X}},43}^{1/2}{\beta }_{-1}^{-1/2}\,{\rm{G}}.\end{eqnarray} \tag{ C4 }$$

The acceleration timescale is also changed as

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}\simeq 2.5\times {10}^{-3}\displaystyle \frac{{E}_{p}}{25\,\mathrm{TeV}}\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}{L}_{{\rm{X}},43}^{-1/2}{\beta }_{-1}^{-1/2}\,{\rm{s}}.\end{eqnarray} \tag{ C5 }$$

The escape timescale is left unchanged and is

$$\begin{eqnarray}&&{t}_{\mathrm{esc}}\simeq \displaystyle \frac{R}{c}\simeq 250\,\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\,{\rm{s}}.\end{eqnarray} \tag{ C6 }$$

For the photohadronic timescale, the change in the X-ray energy density leads to the revised estimate

$$\begin{eqnarray}&&{t}_{p\gamma }\simeq 2\times {10}^{4}\,{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\right)}^{2}\displaystyle \frac{25\,\mathrm{TeV}}{{E}_{p}}{L}_{{\rm{X}},43}^{-1}\,{\rm{s}}.\end{eqnarray} \tag{ C7 }$$

Equating this with the escape timescale, we find the critical energy

$$\begin{eqnarray}&&{E}_{p,\mathrm{cool}}=2.1\,{L}_{{\rm{X}},43}^{-1}\,\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\,\mathrm{PeV}.\end{eqnarray} \tag{ C8 }$$

We similarly find the maximum cooling-limited proton energy as

$$\begin{eqnarray}&&\begin{array}{l}{E}_{p,\mathrm{rad}}\simeq 601\,{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}{L}_{{\rm{X}},43}}\right)}^{1/2}{\beta }_{-1}^{3/4}\\ \min \left[1,5.9{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}{L}_{{\rm{X}},43}}\right)}^{1/2}{\beta }_{-1}^{-1/4}\right]\,\mathrm{PeV}.\end{array}\end{eqnarray} \tag{ C9 }$$

In the peak neutrino spectrum, we have to account for an additional factor (2π) × (4π/3) in the denominator coming from the geometrical factors in u_p ∝ u_B and in the target photon density and a correction factor 4π/3β_rec to account for the relative volume factor between the spherical volume 4π R³/3 and the volume of the layer L³ β_rec, leading to the estimate

$$\begin{eqnarray}&&{({E}_{\nu }^{2}{Q}_{\nu })}^{\mathrm{pk}}\simeq {10}^{40}\ \displaystyle \frac{5{r}_{g}{L}_{{\rm{X}},43}^{2}{\eta }_{p,-1}}{{{RM}}_{7}{\beta }_{-1}}\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\ \mathrm{TeV}}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ C10 }$$

Finally, in the estimated number density of pairs, assuming that the fraction of MeV photons can be estimated just as in the planar case, we find

$$\begin{eqnarray}&&\begin{array}{l}{n}_{e}^{\gamma \gamma }\simeq 3.5\times {10}^{5}\,{L}_{{\rm{X}},43}^{4}{\eta }_{p,-1}^{2}{\left(\displaystyle \frac{{{RM}}_{7}}{5{r}_{g}}\right)}^{-5}\\ \qquad \times \ {\left(\displaystyle \frac{{E}_{p,\mathrm{br}}}{25\,\mathrm{TeV}}\right)}^{2}{\beta }_{-1}^{-2}\,{\mathrm{cm}}^{-3}.\end{array}\end{eqnarray} \tag{ C11 }$$

Notice that in the spherical case, with R ∼ 0.28L as discussed above, the pair number density can be much larger compared to the planar case (see Equation (19)) due to the spherical region being significantly smaller in all directions. Nonetheless, differences found in ${n}_{e}^{\gamma \gamma }$ between the planar and spherical geometry do not affect proton cooling and, in turn, our neutrino estimates.

For completeness, we show in Figure 5 the regions of R − L_{X,[2–10] keV} leading to a pair production comparable with the pair population inferred from Compton opacity, a peak neutrino flux consistent with the IceCube measurements, and 25 TeV < E_p,cool < 100 TeV. Compared to the planar case of Figure 2, for a compact enough region (R ∼ 1–2 r_g ), one can simultaneously produce a sizable number of pairs and saturate the observed neutrino flux without entering the fast-cooling regime (see also Figures 3 and 6). This is due to the more compact spherical geometry, which allows for a larger production even at lower luminosities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We supplement our analytical calculations of the neutrino spectrum with numerical results obtained with the proprietary code athe ν a (Dimitrakoudis et al. 2012). The code solves the kinetic equations for relativistic protons, secondary electrons and positrons, photons, neutrons, and neutrinos. The emitting region is assumed to be spherical with radius R. The magnetic field inside the source has strength B, and all particle populations inside the source are assumed to be isotropically distributed. The physical processes that are included in athe ν a and couple the various particle species are electron and proton synchrotron emission, synchrotron self-absorption, electron inverse Compton scattering, γ γ pair production, proton–photon (BH) pair production, and proton–photon and neutron–photon pion production. The photomeson interactions are modeled based on the results of the Monte Carlo event generator sophia (Mücke et al. 2000), while BH production is modeled using the Monte Carlo results of Protheroe & Johnson (1996); see also Mastichiadis et al. (2005). For the modeling of other processes, we refer the reader to Mastichiadis & Kirk (1995) and Dimitrakoudis et al. (2012). The code also returns the Thomson optical depth of the source based on the density of electron–positron pairs that cool due to radiative losses down to γ_e = 1. We note, however, that these cold pairs are "passive" because radiative processes of nonrelativistic electrons (e.g., cyclotron radiation and Compton scattering) are not included in the numerical code. With the adopted numerical scheme, energy is conserved in a self-consistent way, since all the energy gained by one particle species has to come from an equal amount of energy lost by another particle species. The adopted numerical scheme is ideal for studying the development of in-source electromagnetic cascades in both linear (Zhang et al. 2020) and nonlinear regimes, where the targets for photohadronic interactions are themselves produced by the proton population (e.g., Petropoulou et al. 2014; Mastichiadis & Petropoulou 2021).

As discussed in Appendix C, we consider a spherical region with a radius $R=L{(3{\beta }_{\mathrm{rec}}/(4\pi ))}^{1/3}\simeq 0.28\ L{\beta }_{-1}^{1/3}$. All particles escape from the source on an energy-independent timescale t_esc = R/c. Relativistic protons, after being accelerated into a broken power-law distribution by reconnection (see also Section 3.3), are then "injected" into the spherical region at a rate given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dq}}_{p}}{{{dE}}_{p}}={q}_{0}\left\{\begin{array}{l}{\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-1},{E}_{p}\leqslant {E}_{p,\mathrm{br}}\\ {\left(\displaystyle \frac{{E}_{p}}{{E}_{p,\mathrm{br}}}\right)}^{-s},{E}_{p,\mathrm{br}}\lt {E}_{p}\leqslant {E}_{p,\mathrm{rad}}\end{array}\right.,\end{eqnarray} \tag{ D1 }$$

where ${q}_{0}=(s-2)/(s-1){\eta }_{p}^{* }{{cu}}_{B}{R}^{-1}{E}_{p,\mathrm{br}}^{-2}$, assuming ${E}_{p,\mathrm{rad}}\,\gg {E}_{p,\mathrm{br}}\gg {E}_{p,\min }$. Here, we introduce ${\eta }_{p}^{* }$ to refer to the fraction of energy carried by the proton population in the absence of cooling, as explained in Appendix B. To compute the steady-state proton distribution and the emerging photon and neutrino emission, we study the system for 10R/c; this is also a typical lifetime of the reconnection region for β_rec ∼ 0.1. Results are shown for M = 6.7 × 10⁶ M_⊙, ${E}_{{\rm{X}},\min }=100\,\mathrm{eV}$, ${E}_{{\rm{X}},\max }\,=100\,\mathrm{keV}$, X-ray photon index −2, L_X = 5 × 10⁴³ erg s⁻¹, β_rec = 0.1, η_X = 0.5, ${\eta }_{p}^{* }=0.1$, s = 3, ${E}_{p,\min }={m}_{p}{c}^{2}$, E_p,br = 25 TeV, E_p,rad = 100 PeV, B = 10^5.1 G, and R ≃ 1.4r_g ≃ 1.4 × 10¹² cm.

Figure 6 shows the steady-state proton distribution we obtain numerically (thick blue line) and analytically as explained in Section 3.3 (cyan markers). The analytical expression describes well the spectral break at ${E}_{p}^{* }$ caused by the change in the energy dependence of the photopion loss timescale.

Some key findings from the numerical analysis are as follows.

• 1.  
  Protons with energies E_p ≳ E_p,cool cool efficiently through photomeson interactions (equivalently, the reconnection region is optically thick to photomeson interactions for protons beyond this energy).
• 2.  
  About 50% of the energy initially carried by the proton population into the radiative zone (i.e., the reconnection downstream) is transferred to secondary particles (resulting in η_p ≃ 0.05), with neutrinos taking about 12% of the energy lost by protons.
• 3.  
  We find that the density of pairs with γ_e ∼ 1 is n_e,cold ∼ 3 × 10¹¹ cm⁻³. This value is consistent with the analytical estimate in Equation (C11) for η_p = 0.05 and R = 1.4r_g . Interestingly, it is also close to the pair density needed to establish a Thomson optical depth of order unity in the corona (see Equation (1)).

Protons escaping the reconnection layer can provide an additional contribution to the neutrino emission. Photohadronic neutrino production is generally dominated by the reconnection layer due to the rapid geometric dilution of both the emitted proton population and the X-ray photon field. On the other hand, a reasonable concern is whether proton–proton scattering on the material in the accretion flow can significantly contribute to the neutrino emission. We now show that this is not the case by estimating the energy-loss timescale due to pp scattering on the accretion flow gas.

The proton density in the accretion disk can be estimated as

$$\begin{eqnarray}&&{n}_{p}^{\mathrm{acc}}\simeq \displaystyle \frac{\dot{M}}{4\pi {{rHv}}_{r}{m}_{p}},\end{eqnarray} \tag{ E1 }$$

where H denotes the disk's scale height, and v_r is the radial inflow velocity. Considering that the inner disk is thick, with H ∼ r, due to the disk expansion associated with radiation pressure, and taking v_ϕ ∼ (GM/r)^1/2 and v_r ∼ α v_ϕ as the azimuthal and radial velocity, with α ∼ 0.1 being the standard α parameter from the Shakura and Sunyaev theory of accretion disks, we use Equation (A8) and obtain

$$\begin{eqnarray}&&{n}_{p}^{\mathrm{acc}}\simeq 7\times {10}^{10}\displaystyle \frac{{L}_{{\rm{X}},43}{\phi }_{1}}{\alpha \,{\eta }_{\mathrm{acc},-1}{M}_{7}^{2}}{\left(\displaystyle \frac{L}{5{r}_{g}}\right)}^{-3/2}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ E2 }$$

where ϕ = 10 ϕ₁. We can now estimate the pp energy-loss timescale as (see, e.g., Murase 2018)

$$\begin{eqnarray}&&\begin{array}{l}{t}_{{pp}}={({n}_{p}^{\mathrm{acc}}{\sigma }_{{pp}}{\kappa }_{{pp}}c)}^{-1}\\ \simeq 3\times {10}^{4}\,\displaystyle \frac{\alpha {\eta }_{\mathrm{acc},-1}{M}_{7}^{2}}{{L}_{{\rm{X}},43}{\phi }_{1}}{\left(\displaystyle \frac{L}{5{r}_{g}}\right)}^{3/2}\,{\rm{s}},\end{array}\end{eqnarray} \tag{ E3 }$$

where σ_pp = 3 × 10⁻²⁶ cm² is the pp cross section and κ_pp = 0.5 is the inelasticity. In order to see whether pp interactions can be fast enough to significantly produce neutrinos, we should compare this timescale with the typical residence time of the protons in the inner region of the accretion disk after they escape the reconnection layer. If protons are able to freely stream out of the region, such a timescale is of the same order of magnitude as the escape timescale from the layer; therefore, the pp production is negligible. On the other hand, if proton escape is diffusive, their residence time in the disk would be larger, possibly leading to enhanced pp production. In any case, since only protons below the cooling energy E_p,cool are able to escape the reconnection layer while still retaining their energy, this additional channel of neutrino production would mostly be relevant at energies near or below the neutrino peak and therefore would not be observable at IceCube.

A final possibility that one could consider is that neutrons, produced in the reconnection layer by p γ interactions, could escape from it and decay outside, producing neutrinos. However, since the fraction of neutron energy carried by a neutrino in beta decay is typically of the order of 10⁻⁴ (see, e.g., Lipari et al. 2007; Fiorillo et al. 2021), these neutrinos would appear at much lower energies, in the GeV range. Therefore, we do not account for this contribution in our work.