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

The recent discovery of astrophysical neutrinos from the Seyfert galaxy NGC 1068 suggests the presence of nonthermal protons within a compact “coronal” region close to the central black hole. The acceleration mechanism of these nonthermal protons remains elusive. We show that a large-scale magnetic reconnection layer, of the order of a few gravitational radii, may provide such a mechanism. In such a scenario, rough energy equipartition between magnetic fields, X-ray photons, and nonthermal protons is established in the reconnection region. Motivated by recent 3D particle-in-cell simulations of relativistic reconnection, we assume that the spectrum of accelerated protons is a broken power law, with the break energy 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 dnp/dEp∝Ep−1 below the break and dnp/dEp∝Ep−s above the break, with IceCube neutrino observations suggesting s ≃ 3. Protons above the break lose most of their energy within the reconnection layer via photohadronic collisions with the coronal X-rays, producing a neutrino signal in good agreement with the recent observations. Gamma rays injected in photohadronic collisions are cascaded to lower energies, sustaining the population of electron–positron pairs that makes the corona moderately Compton thick.


INTRODUCTION
Active galactic nuclei (AGN) are among the most energetic sources in the Universe.They are powered by accreting supermassive black holes (SMBHs), 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 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; (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 multi-messenger emission since the late seventies (Berezinsky 1977;Silberberg & Shapiro 1979;Eichler 1979).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, and Padovani et al., in prep.).This neu-Corresponding author: damiano.fiorillo@nbi.ku.dk trino flux has been measured with a significance of 4.2σ in the energy band 1.5-15 TeV, showing a soft neutrino spectrum dΦ ν /dE ν ∝ E −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 are 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 ultra-fast 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) pro-vided the first phenomenological modeling of NGC 1068, prescribing that protons are energized via diffusive shock acceleration (DSA) 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 slope −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 multi-messenger 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 three-dimensional particle-in-cell (PIC) simulations of relativistic reconnection (Zhang et al. 2021(Zhang et al. , 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)1 .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 non-thermal 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 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 Bethe-Heitler 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.

PROPERTIES OF THE RECONNECTION LAYER
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 7 M 7 M ⊙ .In order to maintain our conclusions as general 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 magneticallyarrested disks (Avara et al. 2016;Dexter et al. 2020;Ripperda et al. 2020;Porth et al. 2021;Scepi et al. 2022;Nathanail et al. 2022;Ripperda 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 & Titarchuk 1980;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 trans-relativistic bulk motions of the stochastic chain of reconnection plasmoids (Beloborodov 2017;Sironi & Beloborodov 2020;Sridhar et al. 2021;Sridhar et al. 2023).We refrain from adopting a specific scenario to describe the X-ray coronal emission of NGC 1068, and use instead 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 2 ≃ 1.5 × 10 12 M 7 cm.The cross-sectional area of the "upstream" plasma flowing into the reconnection layer is estimated as A in ∼ L 2 .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 2 /4πρc 2 exceeds unity, where ρ is the plasma mass density.As shown below, this is indeed our regime of interest.
In both X-ray binaries (Reig et al. 2018;Malzac et al. 2001;Zdziarski 2000) 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;Sridhar et al. 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 (henceforth, the electron density), to infer where 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 Noting that the X-ray flux is related to the intrinsic X-ray luminosity as L X = 2F X L 2 , where the factor of two accounts for emission from both sides of the layer, the magnetic energy density is written as where L X = 10 43 L X,43 erg/s is the bolometric intrinsic X-ray luminosity of NGC 1068.This implies a magnetic field Taking the pair number density from the optical depth argument presented above, we can estimate the magnetization of the system as 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".

RELATIVISTIC PROTONS IN THE RECONNECTION LAYER
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 time scales (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.

Spectrum of accelerated protons
Our understanding of the physics of particle acceleration in relativistic reconnection has greatly advanced in recent years thanks to first-principles particle-in-cell (PIC) simulations (Werner & Uzdensky 2017;Zhang et al. 2021;Zhang et al. 2023;Chernoglazov et al. 2023).Motivated by 3D PIC simulations of reconnection in the relativistic regime, we assume that, in the absence of proton cooling losses: (i) protons are accelerated linearly in time with an energization rate of approximately dE p /dt ≃ β rec eBc, (ii) nearly all protons flowing into the reconnection layer are injected into the acceleration process, and (iii) the proton spectrum can be modeled as a broken power law with dn p /dE p ∝ E −1 p below some break energy E p,br , and dn p /dE p ∝ E −s p with s ≳ 2 above the break energy, up to a cutoff energy E p,max ≃ β rec eBL corresponding to the maximum energy that protons can attain in the absence of cooling losses.
The value of the post-break slope depends on various physical conditions, including the strength of the "guide" field (i.e., the non-reconnecting component orthogonal to the oppositelydirected fields).In fact, Werner & Uzdensky (2017) used 3D PIC simulations to show that, as the guide field strength B g increases from B g /B = 0 to B g /B ∼ 1, the high-energy spectral slope increases from s ≃ 2 to s ≃ 3 (so, the spectrum gets steeper).In the following, we deduce the post-break slope from IceCube neutrino observations, which suggest s ∼ 3.
The characteristic energy at the spectral break is associated with the so-called proton magnetization σ p = B 2 /4πn p m p c 2 , such that the break energy is of order E p,br ∼ σ p m p c 2 .This is equivalent to stating that the energy density in the nonthermal proton population, u p ∼ n p E p,br , is comparable to the magnetic energy density, u B .We will later show that 25 TeV is the ballpark value for the break energy needed to explain the IceCube neutrino observations.3.2.Proton acceleration, escape and cooling Zhang et al. (2021) and Zhang et al. (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 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 they 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. 2021Zhang et al. , 2023;;Chernoglazov et al. 2023)-that the prior phase of active energization injects a broken power-law spectrum scaling as ∝ E −1 p below the break and as ∝ E −s p above the break.Pro- tons escape/advect out of the reconnection region on a typical timescale The dominant channel of proton energy loss is photohadronic (pγ) inelastic scattering on the X-ray photon field; we discuss in App.A additional cooling processes which 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 2 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 X,min = 100 eV and E X,max = 100 keV, The frequency range covered by the coronal spectrum is consistent with a scenario in which photons are initially injected in the optical-ultraviolet (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 Sec.3.3, a different value of E X,min would not alter our results.
We assume an effective cross section that is constant in energy, σ = σ pγ K pγ = 70 µb, for interaction energies above the center-of-mass energy threshold ϵ r = 390 (in units of m e c 2 ), 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 Eq. 8, we estimate the pγ energy-loss timescale as for proton energies E p < E * p , which is defined as In contrast, for E p > E * p the cooling timescale flattens to the asymptotic energy-independent value 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 Above E p,cool , protons cool much faster than they can escape the reconnection layer, and therefore lose most of their energy while being 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 App.A) or pγ, yielding (first and second term in the square bracket, respectively) PeV.
Fig. 1 collects the energy-dependent timescales we have introduced, for a benchmark choice of parameters representative of the likely conditions in NGC 1068 (see Sec. 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 App. A. We can grasp directly 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.

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 ∝ E −1 p below the break and as ∝ E −s p above the break.When accounting for cooling losses, the steadystate proton spectrum will be composed of several power-law segments that depend on the hierarchy of the proton energy scales.For E p,br < E p,cool < E * p < E p,rad , as shown in Fig. 1, we expect (14) Notice that the break at E * p , which is determined by the lower energy of the X-ray target E 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 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 η p = (u p /u B ) uncool ≳ 0.1 (French et al. 2023;Totorica et al. 2023;Comisso & Jiang 2023).In what follows, we consider the special case where 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 App.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,br ∼ E p,cool ≪ E * p ≪ E p,rad , with post-break slope of s = 3 (motivated by the neutrino observations), we can compute the total proton number density, where we have replaced the square bracket by its typical order of magnitude value 10 for E p,br ∼ 25 TeV (i.e.typical value needed to reproduce the IceCube neutrino measurements).Eq. ( 15) confirms our initial assumption that n p m p ≪ n e m e , see Eq. ( 1).Charge neutrality therefore requires the largest part of the lepton population in the coronal region to be made of electron-positron pairs.

Neutrino emission
Photohadronic collisions inject charged pions, which subsequently decay to neutrinos.We estimate the production rate of neutrinos and anti-neutrinos of all flavors as where Q ν is the rate of neutrino production differential in neutrino energy, β rec L 3 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 twobin 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 , 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 > E * p .Taking the post-break proton slope to be s = 3 as in Eq. 14, we expect a broken power law with since the proton spectrum is steepened by photohadronic losses to dn p /dE p ∝ E −4 p ), see also Eq. 14.Thus, for our benchmark value s = 3 of the post-break proton slope, the neutrino spectrum is consistent with IceCube observations.The peak value of E 2 ν Q ν , which lies in the energy range 0.05E p,br − 0.05E p,cool , is 5.9 × 10 39 η p,−1 L 2 X,43 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.

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 Bethe-Heitler relativistic pairs (as we will show in Sec.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 Sec. 3.6).An order-ofmagnitude estimate of the MeV photon energy density is and their number density is n MeV ∼ u MeV /1 MeV.While this is only an order-of-magnitude estimate, it captures the  [2,10] keV parameter space leading to a neutrino signal from NGC 1068 consistent with IceCube observations.The axis at the top shows the bolometric intrinsic X-ray luminosity LX , which we use in our equations.We show in red the requirements on the cooling break (Eq.12) and in blue the constraints on the peak neutrino luminosity (Eq.17).We show in green the region where the pairs injected by the attenuation of hadronically-produced MeV photons have a number density from 5% to 100% of the pair density inferred from the Compton opacity.
expected scaling with the energy injected by hadrons; we validate it by comparing with our numerical simulations in Sec.3.6.Finally, we can estimate the amount of pairs produced by MeV photon-photon attenuation as For L X,43 ≳ 5 − 10, the attenuation of ∼MeV photons can therefore lead to a pair population that constitutes a sizeable fraction of the pair density inferred from the Compton opacity in Eq. 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.

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, and Padovani et al., in prep.) and a black hole mass M 7 = 0.67 (see Greenhill et al. 1996, and Padovani et al., in prep.).
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 42 erg/s (based on the best-fit value and 95% confidence region shown in Fig. 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 post-break 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γ ∝ E −1 p in Eq. 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 TeV and 100 TeV (the fast-cooling case E p,cool ≪ E p,br is discussed in more detail in App.B.) Fig. 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 X,[2−10] keV = 3 +3 −2 × 10 43 erg/s (see Marinucci et al. 2016, and Padovani et al., in prep; 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 keV -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 Eq. 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 coronal regions.
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 Fig. 3 (orange solid line), whose peak is in reasonable agreement with the analytical expectation (orange circle)2 .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 Fig. 3.
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 Bethe-Heitler 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 in 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 11 cm −3 , 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 Eq. C30.
We note that the radiative calculations we presented have some limitations.First, cold pairs are "passive" (i.e., nonemitting) particles in our code, since radiative processes of non-relativistic 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 rea-son, 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).

SUMMARY AND DISCUSSION
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 highenergy 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 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 reconnecting magnetic fields.Protons accelerated in the reconnection region will also carry a fraction ∼ 0.1 of the magnetic energy density.To avoid an "energy crisis", the hard ∝ E −1 p proton spectrum established at low energies must soften (∝ E −s p , 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 2 ν Q ν ) pk ∝ M −1 L 2 X , which comes from the direct connection between X-ray, magnetic, and non-thermal 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 by the attenuation of Bethe-Heitler-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 non-thermal proton population may play an important role also 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 Eq. 15) TeV .
(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 by the physical density of the initialized electron-proton torus)3 .To the best of our knowledge, 3D PIC simulations in this regime (1 ≪ σ ≪ σ p ) have yet to properly quantify how the post-break 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 electronpositron 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 selfconsistently 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 similarmost notably, by assuming that Comptonization happens in a magnetically-dominated corona, one can infer the coronal density and field strength directly from X-ray observationsbut 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.

APPENDIX A. BETHE-HEITLER ENERGY LOSSES
In the main text, we have considered as a dominant loss channel the photohadronic interaction of reconnectionaccelerated 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 Bethe-Heitler (BH) interactions with the optical-ultraviolet (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 , and 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 γ (E γ )/dE γ reads A1) where and we take the BH cross section σ ϕe (ϵ) from the numerical fits of Stepney & Guilbert (1983).For the X-ray photon target density given by Eq. 8, which we report here for convenience ) we obtain the energy loss rate as where ξ X (γ p ) is the dimensionless function defined by For comparison, the photohadronic timescale for The ratio between these two timescales does not depend on the X-ray luminosity or the size of the region; noting that the function ξ X (γ p )/γ 3 p has a maximum value of about 2 × 10 −4 , we easily verify that the ratio between the two timescales is bounded by 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 ≤ E * p .Nevertheless, the energy lost via pγ interactions is divided into 3 particle species (neutrinos, gamma-rays, and pairs), while the proton energy lost via Bethe-Heitler is channeled to pairs alone.Therefore, the emission of Bethe-Heitler pairs and neutrinos may be comparable, as shown in Fig. 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 multi-temperature 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 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 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 which is valid in the limit r OUV ≫ L. We can now estimate the BH energy loss timescale off the OUV field as We can again determine the ratio between this energy loss timescale and the photohadronic timescale, noting that the function ξ OUV (γ p )/γ 3 p has a maximum value of about 9 × 10 −25 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 of the low-energy photons injected in 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 which is again much longer that the t pγ loss timescale.

B. IMPACT OF FAST COOLING ON THE PROTON DISTRIBUTION
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 η * p the fraction of magnetic energy density that would go into non-thermal protons in the absence of cooling.Thus, in the absence of cooling effects on spectral slopes, we would have (B14) where s ≃ 3 as in the main text.Notice that, as compared to Eqs. 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 −3 p rather than E −4 p .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 η * p defined here is η p = 3η * p /4. Accounting for cooling, the steady-state spectrum can be approximated as 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 (B16) Therefore, the peak neutrino spectrum, which lies in the range 0.05E p,cool < E p < 0.05E p,br , is namely, the neutrino luminosity is a fixed fraction of the Xray luminosity, independently 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 n γγ e = 1.9 × 10 6 L 4 X,43 η * 2 p,−1 In the fast-cooling regime E p,cool < E p,br , we obtain These results allow us to extend the regions in parameter space shown by Fig. 2 to the regime with E p,cool ≲ E p,br ≃ 25 TeV.We show the corresponding regions in Fig. 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 amount of pairs can be attained in this regime compared to the marginally-fast cooling.
C. SPHERICAL GEOMETRY The numerical code we use for the radiative calculations (see App. 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 as R the size of the spherical region.In relating the planar to spherical case, we will choose an effective size for the spherical region by equating the volumes for the sphere 4πR 3 /3 and the reconnection region β rec L 3 .Therefore, R is related to the full-length of the current sheet as R = L(3β rec /(4π)) 1/3 ≃ 0.28 Lβ The energy density of X-rays is now connected to the luminosity by From the equality we obtain the magnetic energy density The corresponding magnetic field is The acceleration timescale is also changed as The escape timescale is left unchanged and is For the photohadronic timescale, the change in the X-rays energy density leads to the revised estimate Equating this with the escape timescale we find the critical energy We similarly find the maximum cooling-limited proton energy as PeV.
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 /3 and the volume of the layer L 3 β rec , leading to the estimate 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 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 Eq. ( 19)), due to the spherical region being significantly smaller in all directions.Nonetheless, differences found in n γγ e between the planar and spherical geometry do not affect proton cooling, and in turn our neutrino estimates.
For completeness we show in Fig. 5 the regions of the 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 Fig. 2, for a compact enough region (R ∼ 1 − 2 r g ) one can simultaneously produce a sizable amount of pairs and saturate the observed neutrino flux, without entering the fast-cooling regime (see also Figs. 3 and 6).This is due to the more compact spherical geometry, which allows for a larger production even at lower luminosities.
D. NUMERICAL APPROACH 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 (Bethe-Heitler) pair production, 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 Bethe-Heitler 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 non-relativistic 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 non-linear 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 App.C, we consider a spherical region with a radius R = L(3β rec /(4π)) 1/3 ≃ 0.28 Lβ 1/3 −1 .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 Sec. 3.3), are then "injected" into the spherical region at a rate given by where q 0 = (s − 2)/(s − 1)η * p cu B R −1 E −2 p,br , assuming E p,rad ≫ E p,br ≫ E p,min .Here, we introduce η * p to refer to the fraction of energy carried by the proton population in the absence of cooling, as explained in App.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 6 M ⊙ , E X,min = 100 eV, E X,max = 100 keV, X-ray photon index -2, L X = 5 × 10 43 erg s −1 , β rec = 0.1, η X = 0.5, η * 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 12 cm.
Figure 6 shows the steady-state proton distribution we obtain numerically (thick solid line) and analytically as explained in Sec.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: • 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).
• About 50% of the energy carried initially 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.• We find that the density of pairs with γ e ∼ 1 is n e,cold ∼ 3 × 10 11 cm −3 .This value is consistent with the analytical estimate in Eq.C30 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 Eq. 1).

E. NEUTRINO EMISSION OUTSIDE OF THE RECONNECTION LAYER
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 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 s, where σ pp = 3×10 −26 cm 2 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 timescale is of the same order of magnitude as the escape timescale from the layer, and 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 an 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 −4 (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.
Figure2.Regions of the L − L X,[2,10]  keV parameter space leading to a neutrino signal from NGC 1068 consistent with IceCube observations.The axis at the top shows the bolometric intrinsic X-ray luminosity LX , which we use in our equations.We show in red the requirements on the cooling break (Eq.12) and in blue the constraints on the peak neutrino luminosity (Eq.17).We show in green the region where the pairs injected by the attenuation of hadronically-produced MeV photons have a number density from 5% to 100% of the pair density inferred from the Compton opacity.We use benchmark values M7 = 0.67, ηp = 0.05, E p,br = 25 TeV.

Figure 3 .
Figure3.Spectral energy distribution of photons (thick solid yellow line) and all-flavor neutrinos (thick solid orange line) from the reconnection region.Single-flavor neutrino spectra, not including mixing (i.e., flavor oscillations), are also shown (dashed and dotted orange lines).For comparison, we show the estimated peak neutrino luminosity in Eq.C29 (orange circle).The non-thermal cascade spectrum is overplotted with a dotted black line.Other lines indicate emission from different processes (see inset legend) before accounting for γγ pair production.Gamma-ray data by Fermi-LAT(Abdollahi et al. 2020) and upper limits by MAGIC(Acciari et al. 2019) are also included (light and dark green markers respectively).Solid and dashed black lines show the best-fit all-flavor IceCube neutrino spectrum with the 68% confidence interval, adopted fromAbbasi et al. (2022).

Figure 4 .
Figure 4. Same as Fig. 2, accounting for the corrections arising in the fast-cooling regime.We choose η * p = 0.067 to match the value of ηp = 0.05 used in the marginally fast-cooling regime of Fig. 2. Below the thin black dashed line, protons are fast-cooling (E p,cool < E p,br ).

Figure 5 .
Figure5.Same as Fig.2, but for spherical geometry.We show the regions in the parameter space R − L X,[2−10] keV , where R is now the size of the spherical region.

Figure 6 .
Figure 6.Proton (thick blue line) and all-flavor neutrino (thick orange line) energy density spectra normalized to the magnetic energy density.The steady state proton distribution has ηp = 0.05.The spectrum of protons "injected" into the radiative region with η * p = 0.1 (i.e., the reconnection downstream) is overplotted with a transparent solid line.Vertical lines indicate characteristic proton energies discussed in the text.Protons with energies Ep ≳ E p,cool cool due to photomeson interactions with the X-ray coronal photons.Cyan and yellow markers indicate the analytical expression for the steady-state proton and neutrino energy distribution, which agrees well with the numerical solution.