INTRODUCTION

Active galactic nuclei (AGN), namely the compact regions at the center of active galaxies, are powerful emitters of electromagnetic radiation across a broad range of wavelengths. The source of their power is associated with the accretion of matter onto a central supermassive black hole (SMBH). AGN are observed throughout the whole electromagnetic spectrum, from radio up to hard X-rays and often also in gamma rays. The spectral components common to all non-jetted AGN (for jetted AGN, the jet emission might hide these components) are (see, e.g., Padovani et al. (2017)): (1) an optical thermal peak emitted by the accretion disk, (2) a second thermal peak in the infrared associated to a parsec-scale dusty torus, and (3) an X-ray power-law spectrum with an photon index close to -2 associated with a region somehow close to the SMBH, often referred to as the AGN corona. The size and precise location of the corona relative to the SMBH and the accretion disk are matter of debate in the scientific community (see e.g. Cackett et al. 2021, and references therein). Nevertheless, it Corresponding author: damiano.fiorillo@nbi.ku.dk is generally believed that the corona could extend from a few up to hundreds of gravitational radii.

The large amount of power emitted by AGN in the form of electromagnetic radiation has led to their widespread consideration as candidates for the acceleration of high-energy cosmic rays (Berezinsky 1977;Berezinsky & Ginzburg 1981;Silberberg & Shapiro 1979;Eichler 1979). This possibility has recently come to the forefront after the IceCube collaboration reported the measurement, at a 4.2 Ã level of significance, of a high-energy neutrino flux in the 1.5-15 TeV energy range from the active galaxy NGC 1068 (Aartsen et al. 2020;Abbasi et al. 2022). The detection of high-energy neutrinos is usually interpreted as the smoking-gun signature of the interaction of hadronic cosmic rays with matter or radiation; indeed, recent work by Das et al. (2024) has ruled out, in a nearly modelindependent way, the possibility that the neutrinos could have originated from leptonic or beta decay processes. Therefore, the IceCube measurement requires hadronic cosmic rays accelerated up to tens of TeV. On the other hand, the absence of a comparable TeV gamma-ray flux (Acciari et al. 2019) implies that these neutrinos must be produced in a region optically thick to µµ pair production. Such a region could be located close to the central black hole, where the optical and X-ray fields would absorb the gamma-rays from neutral pion decays (associated with neutrinos) and reprocess them to lower energies (Berezinsky & Ginzburg 1981). Thus, the GeV gamma-rays observed by Fermi-LAT (Abdo et al. 2010) would be reasonably produced in a different, outer region, e.g., a weak jet (Lenain et al. 2010), the circumnuclear starburst region (Yoast-Hull et al. 2014;Ambrosone et al. 2021;Eichmann et al. 2022), a large scale AGN-driven outflow (Lamastra et al. 2016), failed line-driven winds (Inoue et al. 2022), or an ultra-fast outflow (Peretti et al. 2023). On the other hand, neutrino emission is strongly inferred to be associated with the corona (see Padovani et al. 2024, for a comprehensive review).

The main question that remains open is what mechanism would accelerate protons in the corona up to the required energies of tens of TeV. Various possibilities have been explored in the literature. In the case of diffusive shock acceleration (Inoue et al. 2020), or gyroresonant stochastic acceleration (Murase et al. 2020), proton acceleration to high energies is limited by photohadronic cooling, which sets the energy scale for the observed IceCube signal. Mbarek et al. (2024) considered proton re-acceleration in turbulence after a pre-acceleration phase in intermittent reconnection layers in the vicinity of the black hole as a possible explanation for the neutrino signal. Magnetic reconnection was considered in Kheirandish et al. (2021), modeled as a power-law injection up to the maximum energies allowed by the acceleration, leading to the conclusion that it was disfavored as an explanation for the neutrinos because of the large number of events expected above tens of TeV. On the other hand, in our previous work (Fiorillo et al. 2024), we showed that a reconnection layer might be entirely responsible for the proton acceleration. This scenario is markedly distinct from the previous ones, in that individual protons could be accelerated to energies much higher than hundreds of TeV without efficient cooling, but the maximum global energy density that the proton population can reach is limited by the available magnetic energy density, which determines the energy scale of the peak in the IceCube neutrino spectrum.

In this work, we revisit the scenario of stationary turbulent acceleration proposed in Murase et al. (2020) incorporating recent advancements from particle-in-cell (PIC) simulations in describing the particle acceleration process. Specifically, particle acceleration is governed by non-resonant interactions with turbulent fluctuations (Comisso & Sironi 2019;Wong et al. 2020;Zhdankin et al. 2020;Bresci et al. 2022), resulting in an energy-independent timescale for the particle acceleration process. Furthermore, we explicitly account for the fact that the volumetric rate of proton energization is bounded by the rate of turbulent magnetic energy dissipation. Adopting these constraints, we find that previously considered scenarios are difficult to reconcile with the high neutrino flux measured by IceCube. In this work, on the other hand, we propose a physically motivated model involving non-resonant proton acceleration in strong turbulence, showing that it is consistent with both the normalization and spectral shape of the IceCube signal.

2. TURBULENT MODEL: CORONAL PROPERTIES Before providing a detailed description of the proton and neutrino spectral features inside the corona, in this section we outline the physical parameters of the turbulent coronal model using order-of-magnitude arguments.

In this model, the corona is a spherical region with radius R, with a magnetic field of strength B and stationary turbulence, surrounding a black hole with mass M = M 7 10 7 M » , consistent with the order of magnitude in Padovani et al. (2024). The gravitational radius of the black hole is r g = GM/c 2 ≃ 1.5 × 10 12 M 7 cm, where G is the gravitational constant and c is the speed of light. The corona contains a population of leptons, with number density n e , and protons, with number density n p . The number density of leptons can be directly inferred by the Compton opacity of the corona, whose characteristic values are in the ballpark of Ä T ≃ 0.5 (e.g. Ricci et al. 2018). This implies a lepton number density

where à T is the Thomson cross section.

In our previously proposed reconnection model (Fiorillo et al. 2024) we assumed that the coronal plasma was predominantly composed of leptons, with n e k n p . The rationale is that, if reconnecting current sheets form in the magnetospheric region, the bulk of the protons are expected to be accelerated at a few r g from the event horizon. In order to explain the neutrino signal, a proton density much smaller than the lepton density estimated from Eq. ( 1) is required, thus leading to the conclusion that the corona must be pairdominated within a reconnection-driven model. On the other hand, in a turbulence-driven model, which is not generally associated with the black hole magnetosphere but rather with the accretion flow, one does not necessarily expect the bulk of the protons to be accelerated; in other words, the proton number density estimated from the neutrino signal refers to the number density of non-thermal protons, n nt p . The latter might be a small fraction of the total amount of protons in the corona, n nt p j n p , but there is no fixed prediction as to what fraction. Thus, the neutrino signal in this case does not allow us to uniquely infer the proton and pair composition of the corona, and one could still have a system with n p ≃ n e . We maintain this assumption throughout the subsequent discussion.

The proton number density impacts the properties of the corona, since it affects the magnetization of the plasma, defined as

where B is the total magnetic field strength. With our assumption that n p ≃ n e , protons dominate the mass density and directly determine the plasma magnetization Ã. The strength of the turbulent magnetic field fluctuations, characterized by the root-mean-squared value of the turbulent magnetic field, ¶B, is generally smaller than the total magnetic field strength. Therefore, we denote

as the magnetization associated with the turbulent field alone. With this definition, we have (v A /c) 2 = à tur /(1 + à tur ) ≃ min(à tur , 1), where v A indicates the Alfvén velocity associated with the turbulent component of the magnetic field. We can also relate the total plasma magnetization à to the commonly used proton plasma beta, ´= 8Ãn p k B T p /B 2 . As in Murase et al. (2020), assuming the virial temperature k B T p = GM m p /3R, allows us to express the proton plasma beta as

In our model, a fraction of the coronal protons are accelerated by scattering off turbulent fluctuations. In the historical resonant picture, particles with a gyroradius Ä g are resonantly accelerated by turbulent fluctuations of a comparable size ¼ ∼ Ä g . The rate of energy gain in this process is proportional to the amount of turbulent power at wavenumbers k ∼ Ä -foot_0 g , resulting in an energy-dependent acceleration rate. However, even within the framework of quasi-linear theory, the interaction between charged particles and turbulent fluctuations, interpreted as waves, is not strongly resonant due to significant broadening by wave damping (Chandran 2000;Demidem et al. 2020). In recent years, it has emerged that non-resonant scattering on turbulent fluctuations with size ¼ ≳ Ä g generally constitutes the dominant driver of particle acceleration. This argument is supported by particle-in-cell (PIC) simulations (Comisso & Sironi 2019;Wong et al. 2020;Zhdankin et al. 2020;Bresci et al. 2022), which showed that the characteristic timescale associated with stochastic particle acceleration, t acc , is essentially independent of particle energy and dominated by the largest scale structures that carry most of the magnetic turbulent energy.

As a typical size for the largest scale structures, we should take the turbulence coherence length (i.e., the turbulencedriving scale), ℓ, which is taken to be a fraction of the coronal size, ℓ = ¸R, with ¸f 1. According to PIC simulations, the diffusion coefficient in momentum space resulting from stochastic particle acceleration is (Comisso & Sironi 2019) D p ≃ 0.1Ã tur p 2 c/ℓ. 1 Therefore, as an order of magnitude for the acceleration timescale, we take

The factor à tur in the denominator accounts for the typical energization in every scattering event.

Turbulent magnetic field fluctuations play an important role also in proton confinement. For protons with gyro-radii Ä g j ℓ, scattering on intermittent small-scale field reversals becomes important in strong turbulence, which constitutes the main focus of our study. On general grounds, the proton mean free path depends on the gyro-radius as

with ¶ > 0, namely longer proton residence times for lowerenergy protons. The exact value of ¶, affected by turbulence intermittency, is a subject of active investigation, though theoretical arguments and simulations suggest 0.3 f ¶ f 0.5 (Lemoine 2023;Kempski et al. 2023). Here, we use ¶ = 1/3, but varying ¶ within the range written above does not alter our conclusions. In an effective random-walk fashion, the residence time within the corona is t esc ≃ R 2 /¼ mfp c. Considering that the minimum escape timescale is R/c, when angular diffusion is so slow that the particles freely stream out of the system, we write the escape timescale as

where E p indicates the proton energy.

The dominant factor limiting the acceleration of protons up to the Hillas limit of the turbulent cascade, E p,max = eBℓ, is their cooling. The bulk of the energy injected in nonthermal protons must peak at a characteristic energy around E p,peak ≃ 20 TeV in order to explain a neutrino peak energy E ν,peak ≃ 1 TeV. Therefore, one must assume that protons above E p,peak are cooled faster than they are accelerated by the magnetized turbulence.

At E p = 20 TeV, there are three competing processes that cool protons, namely pp scattering on the dense coronal matter, photopion scattering on the X-ray coronal field, and Bethe-Heitler (BH) scattering on the optical-ultraviolet (OUV) field produced in the accretion disk. The pp energy loss timescale can be estimated as

Here » p ≃ 0.5 is the typical inelasticity in pp collisions, and à pp (E p ) is the total pp cross section; for our numerical calculations, we use the expression valid above 10 GeV from Aharonian & Atoyan (1996).

The cooling rates due to photohadronic interactions depend sensitively on the electromagnetic spectrum inside the corona. For this work, we use the electromagnetic spectrum shown in Fig. 1 in the form of the photon luminosity dL γ /dE γ , defined as the amount of energy emitted by the AGN per unit time, differential in the photon energy E γ . The multiwavelength AGN spectrum from optical to X-ray is inferred following Marconi et al. (2004) under the assumption of an X-ray luminosity L X = 10 44 erg s -1 , where the X-ray luminosity is integrated from 2 keV to 10 keV. The far-infrared torus component is computed according to the prescription of Mullaney et al. (2011). From the photon luminosity, we extract the photon number density per unit energy BH target pγ target Figure 1. Electromagnetic luminosity emitted from the AGN as a function of the photon energy Eγ. We identify by dashed lines the energies corresponding to the dominant target for Bethe-Heitler (BH) and pγ interactions for a proton with energy Ep = 20 TeV, as identified in the text.

From the photon spectrum, we can directly extract the timescales for photohadronic interactions. The first process we consider is Bethe-Heitler (BH) or pair photoproduction pµ → pe + e -. The timescale for BH energy loss can be obtained by integrating the photon spectrum over the cross section and the inelasticity for the scattering. For this calculation, we use the fits provided in Chodorowski et al. (1992). The dominant scatterers for a proton with a typical energy E p = 20 TeV are photons with an approximate energy E γ ≃ m e m p c 4 /E p ≃ 20 eV, which therefore lie in the high-energy tail of the big blue bump often referred to as the soft-excess.

The second process we consider is inelastic photopion production, often called simply pµ interaction, which usually leads to the production of pions in the final state. At relatively low center-of-mass energies, the interaction is dominated by intermediate resonances, the most prominent of which is the ∆ + resonance, while at higher center-of-mass energies it is dominated by deep inelastic scattering leading to multi-pion production. We model the energy loss rate for pµ scattering following Atoyan & Dermer (2003) and Dermer & Atoyan (2003). In this case, the center-of-mass energy for the scattering is typically higher than in BH interactions, due to the higher threshold for pion production, so for a proton with energy E p = 20 TeV the dominant target photons have an energy E γ ≃ ϵm e m p c 4 /Efoot_1 p ≃ 10 keV, corresponding to the X-ray spectrum originating in the AGN corona, with ϵ = 390.

For completeness, we also account for synchrotron energy losses 2

although these are strongly subdominant compared to all the other cooling channels. We collect all these results together, where the energydependent timescales are evaluated at a characteristic energy E p = 20 TeV:

Note that the acceleration timescale depends on the magnetization associated with the turbulent field alone, Ã tur = ¸BÃ, while the cooling and escape timescales either do not depend on the magnetization or depend on the total magnetization Ã.

To estimate the relevant ranges for these parameters, we first notice that t acc cannot be much larger than t pp , otherwise pp losses, which are mostly flat in energy above a few GeV, would lead to a very soft proton spectrum, impeding acceleration to the required energies. By comparing t acc and t pp in Eq. ( 10), we obtain the constraint à tur ¸≳ 0.18 .

It is reasonable to assume that the magnetic field is controlled by the dynamics of the accretion flow. Therefore, in the vicinity of the black hole, we assume that the coherence length can be as small as the gravitational radius, ¸≳ r g /R. For a corona with radius R = 20 -40 r g , then we obtain ¸≳ 0.03 -0.05. This constraint alone suggests values of à tur ≳ 0.01. We will later find that this range is suggested independently also by the IceCube signal, whose normalization requires comparable values of à tur .

A further consideration is that BH losses should become competitive with particle acceleration at around 20 TeV in order to explain the position of the peak of the IceCube spectrum. From this argument, by equating t acc and t BH in Eq. ( 10), it follows that Using the constraint (11) in Eq. ( 12) results in the size of the corona being bounded by the condition

For the remainder of this study, we adopt the benchmark values M 7 = 1 and R = 20r g . Another essential constraint arises from the global energetics of the corona. Assuming the corona is magnetically powered, the rate of magnetic energy dissipation will set the upper limit on its energetics. This should be regarded as an upper limit on the energy available to energize the non-thermal protons producing the IceCube neutrino signal. We assume that the dissipation of the turbulent magnetic field ¶B occurs on the magnetic reconnection timescale (Comisso & Sironi 2019), which is given by (Comisso 2024) 3

with ¸rec ≃ 0.1 (Comisso & Bhattacharjee 2016;Cassak et al. 2017). Then, the rate of turbulent magnetic energy dissipation is

3 Here, we consider that the mean magnetic field acts as the non-reconnecting (guide) magnetic field component within magnetic reconnecting layers forming in the turbulent environment.

Thus we find

Therefore, for à f 1, the rate of dissipated magnetic energy depends on the magnetization associated with the turbulent field alone, à tur = ¸BÃ, rather than on à and ¸B individually.

Note that the rate of magnetic energy dissipation is comparable with the X-ray luminosity, which is consistent with the scenario of a magnetically-powered corona. However, this is not a required assumption of our subsequent analysis.

In Fig. 2 we show the energy-dependent timescales for the benchmark scenario we will consider throughout this work, with the parameters chosen so that the neutrino spectrum peaks in the correct energy range (see below). The main competition to the acceleration mechanism is BH energy losses, which dominate over acceleration above 20 TeV. pµ and synchrotron losses are subdominant in this energy range. pp is about one order of magnitude slower than the acceleration timescale, and at 20 TeV is comparable with the pµ energy loss. The timescale for proton escape is much longer than the other timescales, and it does not strongly affect the proton spectrum around 20 TeV. Nevertheless, as we will see, it plays a crucial role in shaping the proton spectrum at lower energies.

At this stage, we observe two crucial opposing tendencies: the magnetization associated with the turbulent field of the corona cannot be too low, because otherwise the rate of dissipated magnetic energy cannot account for the observed neutrino luminosity, but neither it can be too large, otherwise BH losses are not efficient enough to counteract proton acceleration (see Eq. ( 10)) and fail to produce a peak in the neutrino spectrum at approximately 1 TeV. These observations concur to point to the range

as the relevant one for the corona of NGC 1068.

COSMIC-RAY TRANSPORT AND NEUTRINO PRODUCTION

In this section, we discuss the dynamics of proton acceleration and cooling beyond the order-of-magnitude estimates presented in the previous section, and compute the neutrino spectrum produced in the corona. Our discussion is phrased in terms of a population of protons, which we describe by their phase-space distribution f p (p), where p is their relativistic momentum, so that the number density of protons per unit energy is f p (p)4Ãp 2 dp. Thus, the usual proton distribution per unit volume and energy is

for ultra-relativistic protons.

In our model, protons are accelerated over the timescale t acc by the magnetized turbulence. This mechanism entails all non-thermal particles being accelerated at the same rate through scattering off large-scale turbulent fluctuations. A natural question is whether this process is indeed the most relevant in establishing a population of non-thermal particles. There are at least two aspects that are not captured by this model.

Firstly, there are rare scattering events which result in large energy gains, much larger than the mean squared momentum change ∆p 2 ≃ Ã tur p 2 . These events lead to faster acceleration for a minority of particles, but their rate is determined by the amount of power in turbulent structures with very large velocity gradients. In Gaussian turbulence models, this power is exponentially suppressed. In realistic settings, the suppression is only of the power-law type, due to the intermittent formation of structures with large gradients. Therefore, lowenergy particles are initially accelerated at a faster rate due to rare scattering events occurring on these structures with large gradients, with a typical spectrum f p (p) ∝ p -4 (Lemoine 2022).

Secondly, particles can also be accelerated in reconnecting current sheets (Comisso & Sironi 2018, 2019;Nättilä & Beloborodov 2021;Comisso & Sironi 2022). If the guide field in the reconnection layers is not much larger than the reconnecting magnetic field, these particles undergo acceleration over a much shorter timescale, t rec ≃ (Ä g /c)/(¸1 /2 B ¸rec ) (Comisso & Sironi 2018, 2019). However, for the magnetization values relevant to this work (à tur ≲ 1), reconnection typically results in a soft spectrum f p (p) ∝ p -s with s g 4 -5 (depending on the guide field). Once these particles escape the reconnection region, they enter the turbulent region; if their energy exceeds the threshold at which cooling dominates over turbulent acceleration, they would quickly be cooled down to lower energy; instead, if their energy is below this threshold, they enter with a spectrum steeper than the ones accelerated in regular turbulence, which as we will see typically predicts f p (p) ∝ p -s with s ∼ 3. Hence, we will neglect the contribution of reconnection to the acceleration of high-energy non-thermal protons. We stress that this is true only in environments with relatively low à p = B 2 /4Ãn p m p c 2 , whereas in the scenario proposed in Fiorillo et al. (2024) we considered a very large value of à p , making reconnection the dominant acceleration mechanism.

Scattering on typical magnetic inhomogeneities, on the other hand, is a slower acceleration process, but it acts equally on all non-thermal particles. In a steady state, this mechanism becomes the most significant, as it tends to generate much harder spectra by energizing all non-thermal particles with the same efficiency. For this reason, here we focus on the simplest scenario in which only the mean energy gain resulting from turbulent stochastic acceleration is accounted for. The intermittency associated with the turbulent fluctuations, as well as the role of reconnection, is neglected. This assumption allows us to treat the turbulent acceleration using a Fokker-Planck approach, because rare events with large energization are not considered. Thus, the evolution of f p (p) can be written as

(19) On the right-hand side, the first term accounts for the diffusive acceleration over the timescale t acc . This form of the acceleration term is analogous to the one conventionally used for gyroresonant mechanisms (see, e.g., Berezinsky et al. (1990) for a derivation), and was indeed already employed in Murase et al. (2020) within this context. However, in the gyroresonant framework the acceleration timescale is energy-dependent, as particle acceleration depends on the turbulent power at a scale corresponding to the energy-dependent gyroradius. In contrast, first-principles PIC simulations have shown that the energization of particles in strong turbulence is dominated by nonresonant interactions with large-scale turbulent structures (Comisso & Sironi 2019;Wong et al. 2020;Bresci et al. 2022), so the acceleration timescale is energy-independent and given by Eq. ( 5).

The second term on the right-hand side of Eq. ( 19) accounts for particle cooling, characterized by the timescale

The third term describes the escape of particles, and the last term, q p (p), describes the low-energy injection of protons in the regions of turbulent acceleration, which comes, e.g., from the pre-acceleration by reconnection layers within the turbulent cascade. The specific form of q p (p) is inessential, since it is only relevant at very low energies, much smaller than the ones to which protons are accelerated. We are interested in the stationary solution of Eq. ( 19). Thus, we solve it without the time derivative. The normalization of the solution would be determined by the low-energy source term q p (p), which however is not directly accessible by observations. However, we can relate the normalization to the rate with which the magnetic turbulence energizes the non-thermal particle spectrum. From Eq. ( 19), we have that the amount of energy transferred to the non-thermal protons per unit time and volume is

obtained by integrating the acceleration term in Eq. ( 19), multiplied by p 3 dp, over all energies. Thus, we can normalize the proton spectrum such that the energy flowing to the protons is a fixed fraction of the rate of magnetic energy dissipation,

Below, we will show results for the extreme case F = 1, to show the maximum proton and neutrino production rates that could be achieved by turbulence. Fig. 3 shows the stationary proton spectrum obtained in this scenario, for varying values of à tur and ¸. At high energies, as expected, the spectrum is rapidly suppressed by the fast cooling. Lower values of à tur and higher values of ¸move this transition at lower energies, since they lead to slower acceleration. At low energies, the dominant process shaping the spectrum is the acceleration. One easily verifies that, in this regime, Eq. ( 19) in steady-state admits two independent solutions, f p ∝ p -3 and f p ∝ p 0 . Below the peak, the solution transitions from the former (at low energies) to the latter (close to the peak), which shows as a "pile-up" region. The position of the transition is determined by the rate of escape, despite the latter being slower than any other process. Indeed, in the absence of escape, Eq. ( 19) does not admit any steady solution, since the number of particles would be perpetually increased by the injection term with no process countering it. The steady solution we find arises from the balance between injection and the slow escape. Note that in the gyroresonant framework, the spectral index in the lowenergy region would be highly sensitive to the turbulence spectrum. Instead, since the interaction between protons and turbulence is non-resonant, the behavior of the distribution in this region turns out to be universal. For neutrinos, the production is dominated by the peak region, and therefore the timescale for escape (which determines the break in the proton spectrum) does not play a substantial role in shaping the neutrino spectrum.

With the proton distribution in hand, we can now determine the neutrino spectrum steadily produced in the corona. To achieve this, we compute the neutrinos produced from the decay of pions originating in pp and pµ collisions. To calculate the spectra of the produced neutrinos, we adopt the fit functions provided in Kelner et al. (2006) and Kelner & Aharonian (2008). For pµ interactions, we consider as a target the photon spectrum from Fig. 1. We assume a distance between the source and the Earth of d = 10.1 Mpc (Tully et al. 2008;Padovani et al. 2024).

The bottom panels of Fig. 3 show the all-flavor neutrino flux at Earth, compared with the neutrino signal observed by IceCube assuming flavor equipartition at Earth (thus we multiply by 3 the IceCube muon neutrino flux). The neutrino signal is dominated by a competition of pp and pµ interactions. pµ interactions produce a harder neutrino spectrum, due to the photohadronic efficiency increasing with energy, and dominate close to the peak of the spectrum, while pp production dominates at lower energies. pp interactions lead to a neutrino flux dΦ ν /dE ν ∝ E -1 ν -the pile-up features in the proton spectrum are mostly washed out in the neutrino production and the neutrino spectrum tracks the proton spectrum f p ∝ p -3 below the break. The pµ neutrino spectrum is hardened as dΦ ν /dE ν ∝ E -2+sγ ν (see, e.g., Winter (2012); Fiorillo et al. (2021)), where s γ is the local spectral index of the photons n γ (E γ ) ∝ E -sγ γ , which in the X-ray coronal spectrum is close to s γ ≃ 2. As ¸increases, the acceleration rate decreases, causing the neutrino spectrum to peak at lower energies. On the contrary, increasing à tur tends to increase the acceleration rate. Notice that even with ¸= 0.1, à tur = 0.03 is unable to match the IceCube signal, so the range à tur ≳ 0.01 again emerges as a natural requisite to fit the normalization of the IceCube spectrum. In turn, for à tur = 0.01, the coherence length, and therefore ¸, cannot be too large, otherwise acceleration becomes too slow and the neutrino signal is underpredicted compared to the IceCube measurement.

As a benchmark case, we now focus on the choice à tur = ¸= 0.3 (we assume again ¸B = 0.3, therefore à = 1). For the adopted benchmark coronal radius, the proton plasma beta is ´≃ 0.03. For this case, it is now helpful to quantify the main global properties of the proton population. We consider now a fraction F = 0.3, which agrees with the IceCube observations, as shown in Fig. 4. The total number density of non-thermal protons, after integrating over energy, is n nt p ≃ 1.9 × 10 4 cm -3 , and as we can see from Fig. 3, is mostly dominated by protons at the peak energy, close to the pile-up region. This should be compared with the density of thermal protons in the corona n p ≃ 2.5 × 10 10 cm -3 , which is about six orders of magnitude higher.

While a comprehensive assessment of the non-thermal to thermal proton ratio is beyond the scope of the present work, a few observations can be made. In the large-amplitude turbulence scenario under consideration, one can estimate the fraction of injected protons by observing that within one ℓ/c, the volume of plasma going through the reconnecting current sheets is given by (Comisso & Sironi 2019) V rec ∼ ¸rec 4ÃR 3 /3. Furthermore, for current sheets with a guide

10 0 10 2 10 4 10 6 Neutrino energy, E ν [GeV] 10 -11 10 -10 10 -9 10 -8 10 -7 10 -6

All-flavor flux,

pp pγ magnetic field comparable to the reconnecting component, as expected when ¶B/B ∼ 1/2, the typical power-law slope for protons with p g p 0 is dN p /dpdV ∝ p - [3,4] (Comisso 2024), where p 0 is the characteristic momentum of the protons processed by reconnection. If protons start to be accelerated by scattering off large-scale turbulent fluctuations at an injection momentum of order p inj ∼ 10 2 p 0 , then current sheets inject a fraction ∼ 10 -4 -10 -6 of the processed protons. Combining this with the fraction of protons processed by reconnection, the overall fraction of thermal protons injected into turbulent acceleration is ∼ 10 -4 ¸rec -10 -6 ¸rec ∼ 10 -5 -10 -7 , in accordance with the normalization adopted for the IceCube neutrino signal. The energy density of the non-thermal proton population is u p ≃ 8 × 10 5 erg/cm 3 . For comparison, the magnetic field energy density of the turbulent field for our benchmark scenario is u B = Ã tur n p m p c 2 /2 ≃ 1.9 × 10 7 erg/cm 3 . Hence, protons have an energy density more than an order of magnitude lower than the magnetic field energy density, suggesting a negligible impact of proton acceleration on the turbulent cascade. On the other hand, if the non-thermal protons were to draw an energy comparable to that of the turbulence fluctuations, damping effects on the turbulent cascade might become significant (Lemoine et al. 2024), potentially resulting in a steepening of the proton energy spectrum.

Finally, let us comment on the possible electromagnetic signatures. pp and pµ interactions produce à 0 , which decay to photons subsequently interacting via µµ scattering with the photon field and producing an electromagnetic cascade. In addition, BH interactions produce e + e -which emit photons via synchrotron and Compton radiation, contributing to the cascade. For this work, we do not delve into a full computation of the output of this cascade as we did in Fiorillo et al. (2024), since the neutrino signal by itself is already informative on the coronal properties. However, we verify explicitly that the corona is sufficiently optically thick to µµ interactions to absorb photons in the GeV range, impeding them from exceeding the upper limits from MAGIC (Acciari et al. 2019).

To show this, we determine the optical thickness Ä γγ for a photon with energy E γ traveling through the corona of radius R. For the interaction cross section of µµ pair production, we use the approximation presented in Eq. (4.7) of Coppi & Blandford (1990). For the reconnection-based model, we use the benchmark size of the corona in the spherical model of Fiorillo et al. (2024), namely R = 1.4 × 10 12 cm. Fig. 5 shows the optical thickness Ä γγ in both scenarios. In the reconnection-based model, the corona is much more compact, leading to a much denser target of photons. However, in the turbulence-based model the assumed photon luminosity is higher by about an order of magnitude, which seems to be somewhat required for BH losses to compete with acceleration at the relevant energies. The two factors cancel out, so that overall the optical thickness for the two models is pretty similar. In both cases, photons above about 10 MeV are strongly attenuated.

4. DISCUSSION Magnetized turbulence is one of the most compelling candidates as a mechanism for accelerating the protons responsible for the neutrino signal originating from NGC 1068. We have shown that this scenario is feasible if strong turbulence ( ¶B ∼ B) with magnetization in the range of 0.01 ≲ Ã tur ≲ 1 characterize the corona, since otherwise the rate of dissipated magnetic energy is insufficient to explain the normalization of the IceCube observed flux. In this regime, the often-used assumption that the interaction between cosmic rays and magnetized turbulence is gyroresonant breaks down, since gyroresonances are strongly broadened, and non-resonant interactions dominate. One can still define a diffusion coefficient for the proton energy, but the corresponding acceleration timescale t acc , namely the typical time for the acceleration of non-thermal protons, is essentially energy-independent, as inferred from PIC simulations (Comisso & Sironi 2019;Wong et al. 2020;Zhdankin et al. 2020;Bresci et al. 2022). This contrasts sharply with the gyroresonant model, where t acc depends on the amount of power in the turbulent cascade at scales resonant with the gyroradius of the particle. On the other hand, in the more appropriate non-resonant model, the largest scales of the turbulent cascade control the acceleration rate of all particles.

From these two main changes -relatively large à tur and acceleration timescale independent of energy -follow significant consequences. The rate of acceleration is much faster than in the scenarios with gyroresonant stochastic acceleration considered previously (Murase et al. 2020). Therefore, in order to explain the soft IceCube neutrino spectrum, BH cooling must effectively compete with the stochastic acceleration of non-thermal protons, which in turn requires more compact coronae compared to prior studies. The typical timescale at which the proton spectrum evolves is of the order of 10 4 s. Note that the turbulent scenario envisioned here is truly stationary, with turbulence sustained steadily by energy injection at large scales. In this sense, our model differs from the turbulent acceleration considered in Mbarek et al. (2024), which concluded that turbulence would produce a spectrum f p ∝ p -5 . This conclusion is appropriate only for the case of transient turbulence, in which particles are rapidly accelerated by intermittent fluctuations. In contrast, if turbulence is maintained in a steady state, as in our case, protons are accelerated into a much harder spectrum with f p ∝ p -3 .

Due to the compactness of the corona, our scenario of turbulence-driven acceleration does not lead to a sizable gamma-ray flux at GeV energies. Thus, neutrinos are the only messenger that at present can pinpoint some of the properties of the non-thermal protons. We leave a more detailed analysis of the gamma-ray emission that would arise from the cascade of gamma-rays injected by photohadronic interactions for future work. Our findings have shown that the neutrino signal by itself is already sufficient to constrain the allowed range for à tur , ¸, and the corona size.

Our deduction of a more compact corona seems required to ensure that the timescale of BH losses is comparable to the short particle acceleration timescale. This is found in agreement with model-independent multi-messenger constraints discussed in Murase (2022). One might consider reducing the acceleration rate by decreasing à tur , which could potentially lead to scenarios with a larger corona (e.g. Murase et al. (2020)). However, crucially, lowering à tur slows down pro-ton acceleration, making it easier to accommodate the peak energy of the neutrino spectrum, but also lowers the rate of dissipated magnetic energy, reducing the energy budget available for neutrino production. For example, for a magnetic field B ≃ 1 kG and a turbulence strength ¸B ≃ 0.01, the Alfvèn velocity associated with the turbulent field is v 2 A /c 2 = à tur = 2¸Bu B /n p m p c 2 ≃ 5 × 10 -5 with n p = 10 10 cm -3 . Hence, for a coronal volume V ≃ 10 43 cm 3 , similarly to what was determined in some prior studies, the rate of dissipated magnetic energy, assuming a turbulence scale of ℓ = 10 14 cm, is L B = ¸Bu B V (v A /ℓ) ≃ 9 × 10 39 erg/s, two orders of magnitude lower than the observed neutrino luminosity, and therefore not sufficient to power the required non-thermal proton energization.

The scenario we consider here is also entirely distinct from our previous proposal of reconnection-driven proton acceleration outlined in Fiorillo et al. (2024). In that study, we considered a much more compact corona in the black hole magnetosphere with very high à p = B 2 /4Ãn p m p c 2 (lowdensity region). The acceleration mechanism in that case is entirely due to magnetic reconnection. Additionally, the reconnection-based scenario is likely to result in multiple flaring events, each associated with the formation of a transient reconnection layer, whose typical lifetime is 100 s. Instead, in the present scenario, turbulence is for all practical purposes stationary. It is an intriguing question whether future measurements can truly disentangle the two scenarios of particle acceleration driven by reconnection or turbulence. The transient nature of the reconnection-based picture would certainly be a useful observational test; while it is unlikely that neutrino measurements could exhibit these transient features, future MeV telescopes would be sensitive to the photons from the electromagnetic cascade which would exhibit a similar temporal variability. Another element of difference, which is perhaps the most interesting for future measurements of the neutrino signal at IceCube-Gen2, is the neutrino spectrum; in the reconnection-based scenario, the signal has a power-law spectrum with f p ∝ p -5 (for strong guide fields, as we hypothesized in Fiorillo et al. (2024)), whereas in the turbulence-based scenario considered here cooling eventually inhibits particle acceleration and therefore the spectrum cuts off exponentially. Overall, despite operating on fundamentally different physical principles compared to the reconnectionbased model, the turbulence-based scenario convincingly accounts for the main features of the neutrino signal observed at IceCube.