© The Authors 2025

1. Introduction

Supermassive black holes in active galactic nuclei (AGNs) have long been regarded as potential accelerators to extreme energies (e.g., Eichler 1979a; Berezinsky & Ginzburg 1981; Protheroe & Szabo 1992), whether in the accretion disk itself (e.g., Dermer et al. 1996; Aharonian & Neronov 2005; Lynn et al. 2014; Kimura et al. 2019a,b), in the magnetosphere (e.g., Levinson & Boldt 2002; Aharonian et al. 2002), around the base of the jet (e.g., Rieger & Mannheim 2000, 2002; Rieger & Aharonian 2009), or more recently, through tidal disruption events (e.g., Farrar & Piran 2014; Guépin et al. 2018; Murase et al. 2020a; Winter & Lunardini 2023).

The recent announcement the arrival directions of (a subset of) high-energy (ϵ_ν ≳ 1 − 10 TeV) neutrinos are correlated with local AGNs, in particular, from the nearby (d ∼ 10 Mpc) Seyfert galaxy NGC 1068 (IceCube Collaboration 2022; Halzen 2023; Abbasi et al. 2025), has placed these theoretical considerations on a solid footing, triggered a surge of interest in these systems, and more generally, paved the way for a novel view of the origin of high-energy neutrinos (Neronov et al. 2024; Padovani et al. 2024). If it is confirmed, the detection of these neutrinos, with a flux higher by an order of magnitude than the gamma-ray flux that would be expected in association, suggests that they are produced through hadronic interactions of protons with energies ϵ_p ≳ 30 TeV in compact environments, where gamma-gamma absorption through pair production becomes high enough to limit the outgoing gamma-ray flux (Murase et al. 2020b). Particle acceleration can be envisaged at accretion shocks, in the magnetized turbulence, through magnetic reconnection, or through a combination of these processes (Murase et al. 2020b, 2024; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fang et al. 2023; Fiorillo et al. 2024a,b; Mbarek et al. 2024; Das et al. 2024; Ambrosone 2024). In principle, the corona is strongly magnetized and turbulent, and stochastic acceleration therefore emerges as a prime candidate. The main uncertainty on the expected ion and secondary neutrino spectra in this context stems from the modeling of particle acceleration and from the description of the corona itself.

Several authors have recently pointed out differences between the magnitude of the neutrino flux inferred from NGC 1068 and that expected from the overall population (i.e., NGC 1068 should be significantly more luminous than other Seyferts, see Waxman 2024) and with the available energy budget of the NGC 1068 corona itself (Inoue et al. 2024). Pending confirmation of the IceCube results, which reach the 4.2σ confidence level so far, we take these observations at face value and provide a fresh and detailed discussion of stochastic acceleration in a turbulent magnetized corona here in order to discuss the conditions under which the overall flux and observed spectral shape can be recovered. To this end, we followed recent insights into the physics of particle acceleration in turbulent and sheared velocity flows to investigate the acceleration rate in the turbulent corona and at the base of the black hole outflow. We payed particular attention to the transport equation in stochastic acceleration, and we addressed the possible backreaction of accelerated particles on the turbulence and its consequences for phenomenology. As we discuss below, this backreaction, which generically leads to the damping of the turbulence and therefore to the modification of the accelerated particle spectra (e.g., Lemoine et al. 2024, and references therein), appears likely based on the high ion energy flux inferred from IceCube data. We also examined the possible reacceleration of particles by shear acceleration at the base of the jet, close to the black hole. We have a hierarchical picture in mind in which particles can be accelerated by differing processes depending on their energy. This picture is motivated by the vast gap that separates the kinetic scales (e.g., the gyroradius of thermal ions r_L = 3 × 10³ cm for a 1 GeV proton in a magnetic field with a strength of B = 10³ G) from the macroscopic scale of the source (e.g., the black hole gravitational radius r_g ≡ GM/c² ≃ 1.5 × 10¹² M₇ cm)¹. In these environments, preacceleration in reconnecting current sheets is generic, and it can energize ions up to a Lorentz factor of about σ, where the magnetization parameter is defined as the ratio of the magnetic and plasma energy density (Werner et al. 2018). Unless σ takes very high values (e.g., Fiorillo et al. 2024a; Mbarek et al. 2024), particles must draw most of their energy from their interaction with the turbulence. However, as they gain energy, their mean free path increases so much that they may eventually become sensitive to the coherent shear structure of the velocity flow around the black hole and are therefore further energized through shear acceleration (e.g., Rieger 2019).

We thus adopt this point of view and focus on the following key questions: The conditions under which these sources can push protons to ∼30 − 300 TeV energies, and the extent to which can they account for the spectrum inferred from IceCube neutrino observations. Correspondingly, we do not seek to reproduce the multiwavelength spectrum observed from the prototype NGC 1068, which has been discussed and modeled to fine levels of detail (Murase et al. 2020b, 2024; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fang et al. 2023; Fiorillo et al. 2024a,b; Mbarek et al. 2024; Das et al. 2024). However, we properly implement the associated rates of energy loss that limit the acceleration of multi-TeV to PeV protons.

This paper is organized as follows. In Sect. 2 we discuss the implementation of particle acceleration, first in a turbulent hot corona, and then, in a sheared velocity flow. We show that various parameters control the ion energy spectra. We discuss the results in Sect. 3 and provide three scenarios that can reproduce the inferred high-energy proton spectra. One scenario considers stochastic acceleration in turbulence, another scenario combines stochastic acceleration with shear acceleration, and the third scenario involves stochastic acceleration, properly accounting for the feedback of accelerated particles on the turbulence. We emphasize that in the first two scenarios, different sources are expected to exhibit (possibly largely) different spectra because of the strong dependence of the derived spectra on the physical conditions (e.g., the Alfvénic velocity v_A). We also argue that self-regulation of the stochastic acceleration by turbulence damping provides a well-motivated satisfactory fit to the inferred spectrum without requiring ad hoc normalization of the flux of high-energy ions. This opens new avenues of research into the physics of these cosmic accelerators. We summarize our findings in Sect. 4.

2. Modeling particle acceleration

2.1. Radiation fields and general setup

We considered a generic setup in which a supermassive black hole with a mass of M_BH and size r_g is embedded in a luminous accretion disk (L_disk) corona environment. For simplicity, the corona that encompasses the main disk was assumed to be quasi-spherical and compact, with a characteristic length of r_co ≤ 100 r_g (e.g., Fabian et al. 2015). In contrast to the disk, the corona was taken to be hot, with proton temperatures up to virial, that is, $T_{\mathrm{p}}\simeq 1.2\times {10}^{11}(30/\widehat{r})$ K, where $\widehat{r}\equiv r/r_{\mathrm{g}}$. While still high, radiative cooling results in electron temperatures that are significantly lower (T_e ≪ T_p), as evidenced by the observed X-ray spectral cutoffs (ϵ_cX) at tens to hundreds of keV (e.g., Kamraj et al. 2018, 2022; Kammoun et al. 2024). In principle, energy dissipation through reconnection in magnetically dominated regions could represent an important means to maintain these electron temperatures (e.g., Di Matteo et al. 1997; Merloni & Fabian 2001; Liu et al. 2002; Sironi & Beloborodov 2020). Characteristic coronal plasma β parameters (β_p = p_gas/p_mag) are expected to be about unity, but they likely decrease with increasing distance from the midplane (e.g., De Villiers et al. 2003, 2005). While AGN corona have commonly been thought to be optically thin (τ_T < 1; Zdziarski 1985; Stern et al. 1995), recent modeling suggested that in some sources, τ_T may reach values up to a few (e.g., Kara et al. 2017; Kamraj et al. 2022). Characteristic Alfvén speeds are about $v_{\mathrm{A}}\simeq 0.15{(30/\widehat{r})}^{1/2}{\beta }_{\mathrm{p}}^{-1/2}c$, which makes efficient stochastic acceleration feasible². In terms of β_p, the thermal plasma β parameter, the corona magnetic field may reach strengths of $B=\sqrt{8\pi n_{\mathrm{p}}{m}_{\mathrm{p}}{c}^2/(3\widehat{r}{\beta }_{\mathrm{p}})}\simeq 4\times {10}^3(30/\widehat{r})M_7^{-1/2}{\beta }_{\mathrm{p}}^{-1/2}{\tau }_{\mathrm{T}}^{1/2}$ G, where τ_T = n_pσ_Tr_c (cf., Padovani et al. 2024). In general, however, there are uncertainties because the physical origin of the corona and its emission are not well understood. In particular, magnetic field estimates based on the detection of coronal radio synchrotron emission (hybrid corona), for example, rather suggest B ∼ O(10 G) on scales of ∼80 r_g (Inoue & Doi 2018; Inoue et al. 2020). It seems conceivable that the magnetic field decays away from the midplane. In the following, we take B = 10³ G as the reference value for the corona in application for NGC 1068. We discuss implications for varying fields, but we already note that most of our results are insensitive to this choice. Unless the magnetic field takes values well in excess of 10⁴ G, proton synchrotron losses can indeed be safely neglected. Furthermore, the rate of acceleration in stochastic acceleration is quantified by v_A and ℓ_c (the outer scale of the turbulence), not by B per se. The strength of the magnetic field enters the expression for the particle mean free path, which controls escape losses and the efficiency of shear acceleration, but with an exponent of about one-third (see Appendix A).

As mentioned above, we focused on the suprathermal proton spectra produced by turbulent and/or shear acceleration and discuss the extent to which this can reproduce a spectrum as inferred from the IceCube observations of NGC 1068 (IceCube Collaboration 2022). This proton spectrum has a generic spectral index s ∼ −3 in the energy interval ∼30 − 300 TeV, and its normalization is such that the suprathermal particle pressure at ∼30 TeV is expected to carry a fraction from ∼0.01 to ∼0.3 of the total plasma pressure (Murase et al. 2020b; Kheirandish et al. 2021; Eichmann et al. 2022; Das et al. 2024). In the figures that follow, we express the energy spectra in units of p_gas and accordingly normalize the butterfly diagram characterizing the inferred proton energy spectrum to 0.3, corresponding to a pressure ratio of 0.1. Because of the experimental and modeling uncertainties, this normalization is itself uncertain to within a factor of a few.

2.2. Particle acceleration in turbulence

Stochastic particle acceleration has received increased attention in recent years, with notable inputs from numerical simulations (Zhdankin et al. 2017, 2019; Isliker et al. 2017; Pecora et al. 2018; Comisso & Sironi 2018, 2019; Kimura et al. 2019a; Wong et al. 2020; Trotta et al. 2020; Bresci et al. 2022; Pezzi et al. 2022; Comisso & Sironi 2022; Pugliese et al. 2023) and theoretical developments (Lynn et al. 2014; Brunetti & Lazarian 2016; Lemoine 2019, 2021, 2022; Lemoine & Malkov 2020; Demidem et al. 2020; Sioulas et al. 2020a,b; Xu & Lazarian 2023).

In phenomenological applications, stochastic acceleration is commonly modeled through a purely diffusive Fokker-Planck equation, which is characterized by a (energy-dependent) diffusion coefficient D_ϵϵ, an approach mostly motivated by simplicity and by arguments following a quasilinear description of weak turbulence theory (e.g., Schlickeiser 1989). In this approach, all information regarding the acceleration rate is encoded in the diffusion coefficient D_ϵϵ, in particular, the acceleration rate ν_acc ≡ D_ϵϵ/ϵ². The general transport equation that describes the evolution of the particle distribution n_ϵ ≡ dn/dϵ then takes the form

$$\begin{array}{c}\hfill {\partial }_t{n}_{ϵ}={\mathcal{L}}_{\mathrm{stoch}}{n}_{ϵ}-{\partial }_{ϵ}\left({\dot{ϵ}}_{\mathrm{loss}}{n}_{ϵ}\right)-\frac{n_{ϵ}}{{\tau }_{\mathrm{esc}}},\end{array}$$(1)

where ℒ_stoch denotes the diffusive operator describing stochastic acceleration, that is,

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{stoch}}{n}_{ϵ}\equiv {\partial }_{ϵ}\left(D_{ϵϵ}{\partial }_{ϵ}{n}_{ϵ}\right)-2{\partial }_{ϵ}(D_{ϵϵ}{n}_{ϵ}/ϵ)\end{array}$$(2)

for the purely diffusive Fokker-Planck scheme. In Eq. (1), ${\dot{ϵ}}_{\mathrm{loss}}$ denotes the proton energy loss rate (here, a negative quantity by convention) associated with Bethe-Heitler pair production pγ → e⁻e⁺, and hadronic p − p and p − γ interactions that lead to neutrino production. These energy losses become significant at energies ϵ ≳ 10 TeV (see, e.g., Murase 2022), and their implementation is discussed in Appendix A.

The timescale τ_esc characterizes the time over which particles escape from the acceleration region through spatial diffusion. In a highly dynamic turbulence, it is important to recall that particle transport can occur both through scattering on magnetic inhomogeneities (with mean free path λ_scatt) and through turbulent transport, with a characteristic diffusion coefficient κ_turb ∼ ℓ_cv_A/3. As the scattering mean free path is about ${\lambda }_{\mathrm{scatt}}\sim r_{\mathrm{L}}^{1/3}{\ell }_{\mathrm{c}}^{2/3}$ (see Appendix A), diffusive escape is in practice governed by turbulent transport. This point seems to have gone mostly unnoticed, even though it brings an important constraint on the geometry of the corona. A large ℓ_c diminishes the acceleration rate, but increases the escape rate, while a high v_A increases both acceleration and escape rates. For clarity, we adopted a diffusion coefficient κ = κ_turb + λ_scattc/3 and discuss the consequences in Sect. 3. The diffusive escape timescale is then written τ_esc ≡ r_co²/(2κ).

Additionally, we also took the influence of particle escape from the accelerating region through advective transport into account. To do this, we bound the integration time to τ_adv, which is defined in terms of the advection velocity v_adv as τ_adv ≡ r_co/v_adv. Thus, we did not seek from the outset a stationary solution to the transport equation, Eq. (1). In our formulation, the proton spectrum achieves stationarity because the injection is itself stationary in time. Specifically, at any given time t, the proton spectrum contained in the turbulent corona consists of particles that have been injected at times t′∈[t − τ_adv,  t] and accelerated for a duration τ ≡ t − t′. This spectrum can be expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{cor}}}{\mathrm{d}ϵ}={\displaystyle {\int }_0^{{\tau }_{\mathrm{adv}}}}\mathrm{d}\tau {\displaystyle {\int }_0^{+\infty }}\mathrm{d}ϵ\prime G_{\mathrm{turb}}(ϵ;ϵ\prime ,\tau )\frac{\mathrm{d}{\dot{n}}_{\mathrm{inj}}}{\mathrm{d}ϵ\prime }\end{array}$$(3)

in terms of the Green function G_turb(ϵ; ϵ′, τ), which is associated with Eq. (1). This determines the probability of a transition from ϵ′ to ϵ over a time interval τ, and of the injection rate spectrum dṅ_inj/dϵ. We assumed that the latter is constant in time, and as discussed below, that it is monoenergetic without loss of generality.

The above implementation accommodates injection of particles at one side of the corona, then advection, or uniform stationary injection throughout the corona equally well. In both cases, dn_cor/dϵ can be regarded as the stationary average spectrum. In the former case, the time-dependent spectrum at location x can be written as dn_turb(x)/dϵ = ∫dϵ′ G_turb(ϵ;  ϵ′, x/v_adv)dn_inj/dϵ′ in terms of the injection density dn_inj/dϵ′ at the initial boundary, and the substitutions x → v_advτ, dn_inj/dϵ′ → τ_adv dṅ_inj/dϵ′ lead to Eq. (3) above, where the time integral represents an average over the corona. In the latter case, particles are injected throughout the corona at a rate dṅ_inj/dϵ′ (independent of x), so that

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{turb}}(x)}{\mathrm{d}ϵ}={\displaystyle {\int }_0^x}\frac{\mathrm{d}x\prime }{v_{\mathrm{adv}}}{\displaystyle {\int }_0^{+\infty }}\mathrm{d}ϵ\prime G_{\mathrm{turb}}(ϵ;ϵ\prime ,(x-x\prime )/v_{\mathrm{adv}})\frac{\mathrm{d}{\dot{n}}_{\mathrm{inj}}}{\mathrm{d}ϵ\prime }·\end{array}$$(4)

Then, dn_cor/dϵ is again recovered by taking the average of dn_turb(x)/dϵ over x throughout the corona.

The advection term is sometimes included in the transport equation as an escape term of the form −n_ϵ/τ_adv, which implicitly implies that particles can spend a time > τ_adv in the acceleration region and thus become accelerated to much higher energies. However, advection is a systematic and not a random loss, and it must therefore be treated differently. This impacts the particle spectra in a significant way. In the present description, in particular, advection effectively cuts off the spectrum at a maximum energy set by ν_accτ_adv ∼ 1 (or lower, if other losses dominate).

Recent numerical experiments of stochastic particle acceleration, which either tracked test particles in a magnetohydrodynamic (MHD) simulation or relied on a full kinetic description using the particle-in-cell (PIC) technique, have revealed a richer and more complex landscape than the simple Fokker-Planck formulation above. These simulations indicate a diffusion coefficient of the form D_ϵϵ ≃ 0.2γ²(v_A/c)²c/ℓ_c for particles with a Lorentz factor γ ≡ ϵ/mc², where v_A denotes the Alfvén velocity in the total magnetic field, and ℓ_c is the coherence length of the turbulence. However, they also demonstrated that the purely diffusive Fokker-Planck approach does not successfully account for the particle spectra (Isliker et al. 2017; Pecora et al. 2018; Lemoine & Malkov 2020), unless an energy-dependent (as well as v_A-dependent) advection coefficient is added, which must be extracted from numerical simulations (Wong et al. 2020). This is not a trivial issue because the results of these kinetic simulations must be extrapolated far in time in order to make the connection with phenomenology. For instance, for our fiducial values v_adv = 0.02 c, r_co = 30 r_g and ℓ_c = 10 r_g, the advection timescale reads τ_adv ≃ 150ℓ_c/c, while PIC simulations typically run over ≃10ℓ_c/c.

In Lemoine (2021, 2022), we argued that particles are accelerated through a generalized Fermi process, in which particles gain energy as they cross intermittent regions of dynamic, curved, and/or compressed magnetic field lines. This model was benchmarked on MHD simulations in a regime relevant to black hole coronae, namely a sub- to mildly relativistic Alfvénic velocity, and the time-dependent Green functions for the spectra obtained in this approach also match those seen in PIC simulations in similar conditions (Comisso & Sironi 2022). To account for the intermittent nature of the accelerating structures, it proves necessary to go beyond a purely diffusive Fokker-Planck treatment and to consider the full probability distribution function of the random forces acting on the particles. The transport can then be modeled using Eq. (1) above, with a more general diffusion operator characterized by

$$\begin{array}{c}\hfill {\mathcal{L}}_{\mathrm{stoch}}{n}_{ϵ}\equiv {\displaystyle {\int }_0^{+\infty }}\mathrm{d}ϵ\prime [\frac{\phi (ϵ|ϵ\prime )}{t_{ϵ\prime }}{n}_{ϵ\prime }(t)-\frac{\phi (ϵ\prime |ϵ)}{t_{ϵ}}{n}_{ϵ}(t)],\end{array}$$(5)

where φ(ϵ|ϵ′) represents the kernel describing the probability of jumping from ϵ′ to ϵ over a time interval t_ϵ′, as described in Lemoine (2022) (see also Isliker et al. 2017 for related considerations). This kernel was determined through dedicated measurements carried out in an MHD simulation with an Alfvénic velocity v_A = 0.4 c. It can nevertheless be extrapolated to other Alfvénic velocities by rescaling its width in proportion to (v_A/0.4 c)², and we do this below.

In order to leverage these recent findings and to gauge the influence of the underlying theoretical uncertainty on the predicted particle spectra, we adopt both pictures in the following, that is, we integrate the transport equation, Eq. (1), using either the Fokker-Planck kernel (Eq. (2), hereafter referred to as model (1)) or that describing the generalized Fermi process (Eq. (5), hereafter model (2)). Our main motivation here is to provide a more accurate description of the proton and secondary neutrino spectra and of the corresponding acceleration rate.

For practical applications, the generalized Fermi process produces time-dependent Green functions that take the form of broken power laws in the absence of escape and losses. These solutions obviously differ from the lognormal shape (e.g., Kardashev 1962) that follows from a purely diffusive Fokker-Planck kernel with D_ϵϵ ∝ ϵ² (in the absence of escape and losses). When a continuous injection of particles is included, both formulations nevertheless lead to hard spectra beyond the injection Lorentz factor γ₀ (Schlickeiser 1984), up to a maximum energy that is either dictated by the finite time (τ_adv) over which particle acceleration takes place, or by energy losses.

We do not discuss the physics of injection here and simply assumed that a given number of particles was injected into turbulent acceleration with a characteristic Lorentz factor γ₀ ∼ 1. The value of this Lorentz factor, or the exact shape of the injected spectrum, is not relevant because turbulent acceleration effectively loses the dependence on the initial conditions when an acceleration to high energies has taken place. This injection could either take place through magnetic reconnection in localized patches of sufficient magnetization (e.g., Mbarek et al. 2024) through fast, small-scale turbulent acceleration in localized regions, as observed in some simulations (e.g., Pecora et al. 2018; Trotta et al. 2020), or simply because a fraction of the energy density is contained in supra-thermal particles accelerated elsewhere and accreted into the corona. The exact energy fraction contained in these particles is relevant, however, as we discuss in the following.

The spectral shape is controlled by the following parameters: D_ϵϵ or φ(ϵ|ϵ′) for the acceleration rate. They both depend on v_A, ϵ and ℓ_c; τ_adv ≡ r_co/v_adv, which sets the integration time before advective escape at velocity v_adv, removes the particles τ_esc and τ_loss, which control the timescales associated with diffusive escape and radiative energy losses. We followed previous treatments for the latter and recall their value in Appendix A. Figure 1 shows the results obtained for the Fokker-Planck model (solid lines) and the generalized Fermi process (dashed lines) for various choices of v_A (upper panel), ℓ_c (middle panel), and v_adv (lower panel). The butterfly diagram shows the range of values inferred from the IceCube detection of neutrinos from NGC 1068, as discussed above. The energy spectra shown here is expressed in units of the background plasma pressure (Sect. 2.1). However (see also below), the spectra shown here were normalized in an ad hoc way to a total energy density of unity in these units. In general, turbulent acceleration allows for a favorable situation by yielding sufficiently hard spectra, which effectively minimizes the required energy budget if the peak occurs close to the observed energy range.

Fig. 1. Proton energy spectra ϵ2 dncor/dϵ (per log-interval of energy) vs. proton energy ϵ predicted by stochastic acceleration in a turbulent corona, starting from mono-energetic protons with a Lorentz factor of γ0 ∼ 1. In each panel, the solid line corresponds to a solution obtained by integrating the Fokker-Planck equation up to time τadv = rco/vadv (model (1) including energy losses, advection, and escape as described in the text), and the dashed line shows the corresponding solution for the same parameters for the generalized Fermi model described by Eq. (5) (model (2)). In all panels, rco is set to 30rg, and in the top and middle panels, vadv ≃ 0.02 c. The upper panel shows the dependence on vA, assuming otherwise a coherence length scale for the turbulence ℓc = 10 rg and an advection velocity vadv = 0.02 c. From left to right, as indicated, vA/c = 0.05, vA/c = 0.10, vA/c = 0.20, and vA/c = 0.40. The middle panel examines the dependence on ℓc, assuming otherwise vA/c = 0.20 and vadv/c = 0.02. From right to left, as indicated, ℓc = 3 rg, ℓc = 10 rg, and ℓc = 30 rg. The bottom panel examines the influence of advection, characterized here by the advection velocity vadv (which sets the advection timescale τadv ≡ rco/vadv), assuming otherwise vA/c = 0.20 and ℓc = 3 rg. From right to left, as indicated, vadv/c = 0.02, vadv/c = 0.05, vadv/c = 0.16 and a (more extreme) case vadv/c = 0.49. The butterfly indicates the range of values needed to reproduce a neutrino spectrum as observed by IceCube for NGC 1068. The energy density spectrum is expressed in units of the background plasma pressure (see Sect. 2.1). In these units, each individual spectrum has been normalized to an integrated fraction of unity for the sake of illustration.

The trivial observation that hardly any of the spectra fit this butterfly diagram can be misleading. The important point is that the predicted spectra are extremely sensitive to the choice of parameters. This is easily understood, noting that for D_ϵϵ ∝ ϵ², the mean energy increases exponentially fast in time, that is, ⟨ϵ⟩∝exp(4ν_acct), with ν_acc ≡ D_ϵϵ/ϵ² the acceleration rate defined earlier. Consequently, in a fixed time interval t = τ_adv over which ν_acct ≫ 1 (needed to achieve acceleration to high energies), changing ν_acc by a factor of a few changes the maximum energy by orders of magnitude. The dependence of ν_acc on the parameters is recalled above, ν_acc ∝ v_A²/ℓ_c. We assumed δB/B ≳ 1, otherwise ν_acc should be further suppressed by a factor (δB/B)². The dependence on the advection velocity is also easily understood because the number of e-folds of the increase in the maximum energy scales as ν_accτ_adv ∝ 1/v_adv. From this, we can also infer the dependence on the scale of the corona. The spectra that reach PeV energies are cut off because of catastrophic hadronic losses. Other spectra are cut off because of the limited time they spend in the corona.

Clearly, parameter values that fit the inferred butterfly diagram might be found, and examples are shown below. This comes at the price of some fine-tuning, however. One important corollary is that under the above assumptions, we do not expect all AGNs, even those with a similar X-ray luminosity, to produce neutrinos with a similar spectrum. If the physical conditions that dictate the rate of acceleration (e.g., v_A) show an extended distribution over the neutrino-emitting AGN subpopulation, the population as a whole is instead expected to produce an extended neutrino energy spectrum, from ≲GeV to ∼0.1 PeV. This has consequences for our prediction of the diffuse flux from this population of AGNs and for the modeling of the spectrum of a given source, as we discuss further in Sect. 3.

Another important lesson drawn from Fig. 1 is that fitting the inferred proton spectrum raises an energy budget issue, as noted earlier (Sect. 2.1), which also represents an intriguing coincidence. Namely, in order to match the flux, we must inject a tiny fraction of the thermal plasma into stochastic acceleration, such that after energization by a factor ∼10⁵, its pressure is comparable to that of its parent population. If more particles had been injected, the energy requirement would have become excessive, while in the opposite case, the proton flux falls short of the required flux. Moreover, if suprathermal particles carry an energy density comparable to that in the thermal component, they must backreact on the flow, and in particular, damp the turbulence that feeds them in energy. In this situation, stochastic acceleration becomes nonlinear and self-regulated by damping (Eilek 1979; Eichler 1979b; Litvinenko 2012; Kakuwa 2016; Lemoine et al. 2024). We discuss this possibility further in Sect. 3.

Alternatively, the turbulent corona might not shape the entire proton spectrum, but only part of it, and protons escaping the corona might be further accelerated in the sheared velocity flows at the base of the jet. We first discuss how this shear acceleration can be implemented, and we then examine the general combination of the two scenarios in Sect. 3.

2.3. Particle acceleration in sheared flows

Accreting black hole systems are known to drive relativistic jets and winds, and thereby provide an environment that is conducive to shear-type particle acceleration (e.g., Kimura et al. 2018; Rieger 2019; Lemoine 2019; Rieger & Duffy 2022; Webb et al. 2023; Wang et al. 2024). This might provide a complementary mechanism for proton energization if turbulent acceleration in the corona were to saturate at multi-TeV energies. To be efficient, shear acceleration commonly requires seed injection of energetic particles, which is in line with a hierarchical acceleration scenario. In our context, we envisaged that the black hole vicinity supports a fast, sheared outflow; for example, that the hot X-ray emitting corona forms the base of an outward-moving jet, which may be supported by radiation pressure or be driven magnetically (e.g., Beloborodov 1999; Malzac et al. 2001; Merloni & Fabian 2002; King et al. 2017; Liska et al. 2022; Poutanen et al. 2023; Dexter & Begelman 2024; Sridhar et al. 2025). For the prototype Seyfert galaxy NGC 1068, jet-like features are indeed apparent on larger scales, although jet speeds are most likely nonrelativistic (≲0.1 c) on scales of several dozen parsec (∼10⁷ r_g) from the nucleus (Gallimore et al. 1996; Roy et al. 2000; May & Steiner 2017). Part of its observed gamma-ray emission may in fact be related to the jet (e.g., Lenain et al. 2010; Salvatore et al. 2024). The radio-inferred (large-scale and time-averaged) jet power is about L_jet ∼ 10⁴³ erg/s (Padovani et al. 2024). On the other hand, the minimum (isotropic-equivalent) cosmic-ray power required to account for the observed neutrino emission is about several 10⁴² erg/s (e.g., Murase 2022, see also Appendix A). If part of this emission is indeed related to the outflow, even if it were relativistically boosted, it seems conceivable that cosmic-ray acceleration may contribute to the deceleration of the flow (viscous drag) as it propagates outward. We did not include this possible backreaction and self-regulation of the acceleration process in our modeling, but we note that it may impact on the bulk flow profile (e.g., shear broadening) and on the turbulence generation.

For reference, we explore a model below in which turbulent acceleration in the corona saturates at γ_p ≤ 10⁴ (10 TeV) either due to v_A ≲ 0.1 c or ℓ_c ≫ r_g in the corona. The corona thus sustains a reservoir of energetic seed particles that can diffuse into the outflow and experience further jet-shear acceleration. In this process, acceleration is related to a particle that effectively samples the velocity difference while it is scattered across a shearing flow (Rieger 2019). For simplicity, we considered a collimated mildly relativistic outflow or jet, with a Lorentz factor of γ_j ≡ 1/(1 − u_z²/c²)^1/2 ≲ 3, characterized by a lateral velocity shear u_z(r) with a width of $\mathrm{\Delta }r={\widehat{r}}_{\mathrm{sh}}{r}_{\mathrm{g}}$. Energetic protons undergoing shear acceleration are subjected to radiative losses and diffusive escape (see Appendix A). For the phenomenological application, the acceleration process can be modeled through a diffusive Fokker-Planck equation for n(ϵ) of the type Eq. (1), with mean momentum diffusion coefficient $D_{ϵϵ}\propto p^2{\overline{D}}_{\mathrm{p}}$, where for energetic protons, ϵ = pc = γm_pc². The coefficient ${\overline{D}}_{\mathrm{p}}\equiv a_{\mathrm{g}}\lambda /c$ can be obtained from the local (r-dependent) shear coefficient D_p(r) = (1/15)γ_j(r)⁴(∂u_z/∂r)²λ/c by suitable spatial averaging, for example, a_g = (2/3)(c/Δr)²/w with w = 116 (ln[(1 + β₀)/(1 − β₀)])⁻² in the case of a linearly decreasing shear flow profile with a maximum (spine) speed β₀ = (1 − 1/γ_j, max²)^1/2 (Rieger & Duffy 2019, 2022). In the former expression, λ = cτ_c = η (Δr)(r_L/Δr)^1/3 denotes the scattering mean free path (r_L the gyroradius), for which we took a scaling with η = 2, following Rieger & Duffy (2019) and Wang et al. (2023) (see also Appendix A). Particles can escape from the acceleration region either due to cross-field diffusion (laterally), that is, τ_esc = (Δr)²/(2κ) with κ = λc/3, or through advective transport (along z), τ_adv, z. The accelerated particle distribution on scales of the corona is of interest here. While shear acceleration likely continues to operate along the jet, we considered the former to dominate hadronic interactions and thus neutrino production. Accordingly, we approximated advective escape by ${\tau }_{\mathrm{adv},z}\simeq z/⟨u_z⟩=\stackrel{\sim }{z}{r}_{\mathrm{co}}/⟨u_z⟩$, with ⟨u_z⟩ the mean velocity (corresponding to β₀c/2), and $\stackrel{\sim }{z}\equiv z/r_{\mathrm{co}}$ lower than a few here to ensure that the target photon density is high enough. As before (Sect. 2.2), this was implemented as a constraint on the time available for particle acceleration. Figure 2 presents exemplary proton energy spectra as a function of time (in units of τ_c(γ₀)) assuming that seed protons are continuously injected with a Lorentz factor of γ₀ = 10⁴ into a shearing flow with γ_j = 3, B = 50 G, ${\widehat{r}}_{\mathrm{sh}}=2$ and $\stackrel{\sim }{z}=2$, assuming a corona size of r_co = 30 r_g.

Fig. 2. Exemplary time-dependent proton energy distribution obtained from jet-shear acceleration with a Fokker-Planck model assuming continuous injection of seed protons with γ0 = 104. As a result of advective escape, the predicted particle spectrum will appear somewhat softer (red, t = τadv, z) compared to the otherwise quasi steady-state distribution (yellow). The gray shaded region illustrates the spectral characteristics required to reproduce the neutrino emission for NGC 1068 as measured by IceCube.

In general, the predicted spectra are sensitive to the choice of parameters, for example, the time available for particle acceleration (set by τ_adv, z). One important corollary of the latter is that shear modeling imposes constraints on the size and geometry of the corona. In addition, the required maximum outflow speed to model the CR parent distribution behind the IceCube neutrino spectrum is sensitive to the flow profile. In comparison to power-law or Gaussian-type velocity profiles, a linearly decreasing shear flow profile requires higher speeds to reproduce hard particle spectra. For the former, γ_j ≤ 2 are sufficient (Rieger & Duffy 2022). This suggests that mildly relativistic outflow speeds are sufficient to account for an inferred proton spectrum such as in NGC 1068. Jet deceleration due to cosmic-ray-induced viscous drag, interaction with a central magnetic or cosmic-ray-driven wind, along with the efficient generation of turbulence and entrainment of ambient material, may explain the relatively low velocity of the jet in NGC 1068 on larger scales. Our current model at this stage remains exploratory in nature since the innermost outflow characteristics in Seyfert-2 AGNs remain uncertain, and since we neglected a possible shear contribution due to flow rotation. The latter may facilitate transport of angular momentum from the underlying accretion flow, and it night in effect contribute to further spectral hardening (Rieger & Mannheim 2002). In general, while particle acceleration in sheared flows will not lead to a canonical power-law shape, faster flows will (ceteris paribus) result in harder particle spectra. This agrees with the above corollary, according to which not all X-ray bright Seyfert AGNs are expected to have similar neutrino spectra.

3. Discussion

In the preceding section, we have examined the general features of particle acceleration in a turbulent corona, with possible reacceleration in the sheared velocity flow outside of this corona. We interpret these results here and combine them to suggest ways in which the generic proton spectrum inferred from IceCube observations can be reproduced, and we draw conclusions regarding the phenomenology. More specifically, we investigated the following three scenarios.

First, we considered pure stochastic acceleration in the turbulent corona, as modeled in Sect. 2.2, and we demonstrate that a reasonable set of parameters can be found to reproduce the high-energy general slope. This was known from previous studies (e.g., Murase et al. 2020b; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fiorillo et al. 2024b); the main novelty here lies in the transport equation we used, in particular, for the implementation of generalized Fermi acceleration, and in the treatment of escape and advection. As noted above (Sect. 2.2), matching the spectrum requires some tuning in terms of the physical parameters and in terms of flux normalization. This suggests that other sources might display vastly different neutrino spectra even for otherwise similar coronal X-ray luminosities. Interestingly, the IceCube experiment reports differing spectra from other Seyfert galaxies (e.g., Abbasi et al. 2025, 2024). In the case of NGC 4151 (d ∼ 16 Mpc), for instance, the characteristic neutrino energy seems higher by a factor ∼5, and the high-energy slope is slightly harder than in NGC 1068.

The spectra obtained in this first scenario are shown in the top panel of Fig. 3. While the parameters are the same for the two spectra shown, namely r_co = 30 r_g, ℓ_c = 10 r_g, v_adv = 0.02 c and v_A ≃ 0.2 c, their flux normalizations are different. The ratio of integrated energy density of suprathermal particles to background plasma pressure is ≃4.3 for model (1) and 1.0 for model (2). Because the ratio of the pressures is lower by a factor 3, the dynamical influence of the suprathermal particles should not be ignored in model (1). This further motivates the need to study the impact of backreaction on the turbulence, which is discussed further below. Finally, we note that other sets of parameters might be found, for example, a lower Alfvén velocity with a smaller coherence length scale. However, the general spectral shape should be preserved.

Fig. 3. Proton energy spectra dup/dlnϵ ≡ ϵ2 dn/dϵ vs. proton energy ϵ for three characteristic scenarios described in the text. For each panel, the energy density distribution is shown in units of the background plasma pressure. Top panel: Stochastic acceleration in the turbulent corona, with a characteristic Alfvénic velocity vA ≃ 0.2 c (see text for details), for model (1) and (2), corresponding to the solid and dashed lines, respectively. Middle panel: Combination of stochastic acceleration in the corona with shear acceleration at the base of the jet (see text for the details of the physical parameters). The stochastic contribution is indicated by dotted lines, the shear contribution by dashed lines, and the sum of the two by solid lines. The different colors correspond to models (1) and (2). For the top and middle panels, the normalization of the energy density in accelerated particles has been fixed in order to reproduce the inferred proton flux. Bottom panel: Modeling of stochastic acceleration, accounting for self-regulation by particle feedback on the turbulence (damping), for a characteristic Alfvén velocity vA ≃ 0.2 c. The plasma β parameter is 0.5 (see Appendix A).

In a second scenario, we investigated the interplay of turbulent and shear acceleration. Shear acceleration requires high-energy seed particles because the rate at which acceleration proceeds increases with the particle mean free path (e.g. Rieger 2019; Rieger & Duffy 2022; Webb et al. 2023; Wang et al. 2024), so that in the absence of turbulent preacceleration, this mechanism would be effectively inefficient. In our setting, particles that have been preaccelerated in a turbulent corona are assumed to enter a region of sheared velocity flow, either through diffusive escape or through advection from the corona. For clarity, we concentrated on the latter case, in which particles first cross a turbulent corona and then enter a region of sheared velocity flow. There is a sizable uncertainty here, not only on the overall geometry of the problem, but also on the exact fraction of particles that are able to transit from one region to the next, as well as on the parameters characterizing the shear flow. We thus used the same parameters as in Sect. 2.3 to describe the base of a mildly relativistic jet. To model the particle spectra, we convolved the Green function, which characterizes shear acceleration G_shear(ϵ;  ϵ′, τ_shear), that is, the probability of jumping from ϵ′ to ϵ over the timescale τ_shear spent in the shear region, with an injection distribution that describes the time-dependent particle spectrum of the turbulent region. Denoting the latter as dn_turb/dϵ, we write the stationary spectrum in the shear region as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{shear}}}{\mathrm{d}ϵ}=\frac{1}{{\tau }_{\mathrm{adv},z}}{\displaystyle {\int }_0^{{\tau }_{\mathrm{adv},z}}}\mathrm{d}\tau {\displaystyle \int \mathrm{d}}ϵ\prime G_{\mathrm{shear}}(ϵ;ϵ\prime ,\tau )\frac{v_{\mathrm{adv}}}{v_{\mathrm{shear}}}\frac{\mathrm{d}{n}_{\mathrm{turb}}}{\mathrm{d}ϵ\prime },\end{array}$$(6)

with τ_adv, z the advection time in the sheared region, as before, and

$$\begin{array}{c}\hfill \frac{\mathrm{d}{n}_{\mathrm{turb}}}{\mathrm{d}ϵ\prime }={\displaystyle \int \mathrm{d}}ϵ″G_{\mathrm{turb}}(ϵ\prime ;ϵ″,{\tau }_{\mathrm{adv}})\frac{\mathrm{d}{n}_{\mathrm{inj}}}{\mathrm{d}ϵ″},\end{array}$$(7)

where dn_inj/dϵ″ characterizes (Sect. 2.2) the injection spectrum at the entry of the turbulent region, as before. We introduced a ratio v_adv/v_shear as a prefactor in Eq. (6) to account for the difference in advection velocity of each region (conservation of particle current). We assumed that the shear region extends over r_co and took ⟨u_z⟩=β₀ c/2 ≃ 0.5 c.

Finally, the total stationary spectrum was obtained as the sum of the spectra in each region, dn_cor/dϵ (Eq. (3)) and dn_shear/dϵ (Eq. (6)). Although the overall normalization was left here as a free parameter, the relative normalization of these two was fixed by the volume of each region. The result is shown in the middle panel of Fig. 3. This figure illustrates that the shear contribution acts as re-energization of the seed population produced by turbulent acceleration. This shapes an effective power-law tail close to ∝ϵ⁻³ (per number), which is close to the tail inferred from IceCube observations. The parameters we adopted are the same as in Sect. 2.3 for the shear layer, and v_A = 0.20 c, v_adv = 0.02 c, ℓ_c = 10 r_g, and r_co = 30 r_g as in Fig. 1 for the stochastic acceleration. For simplicity, we ignored possible boosting effects for the shear contribution. In general, if the shear flow were subrelativistic, shear acceleration would become too slow to have any effect, as quantified by the dependence of the diffusion coefficient on the Lorentz factor of the flow (Sect. 2.3). We also note that we considered a scenario in which particles were injected at one boundary of the corona, accelerated by turbulence, and then advected into the shear layer. If we had considered that a fraction of particles could be advected into the black hole instead of the jet, the shear contribution would have been reduced. Similarly, if we had assumed that particles were injected into stochastic acceleration uniformly throughout the corona (see Sect. 2.2), then the spectrum injected into the shear process would have been the steady-state spectrum dn_cor/dϵ instead of the spectrum discussed above. The overall shear contribution would have been effectively reduced in this case as well, and, a normalization of the overall spectrum to the inferred values without invoking boosting would have implied stronger energy requirements. In our case, the ratios of the integrated energy density in suprathermal particles to background plasma pressure are 4.1 and 2.7 for models (1) and (2), respectively.

Finally, we investigated the consequences of self-regulation of the acceleration process by turbulent damping. Because in the two scenarios discussed previously, the pressure carried by the accelerated particles is comparable to that in the background plasma, this last scenario appears to be quite likely. We also note that in the above cases, the flux normalization was set in such a way as to reproduce the normalization inferred from IceCube observations, and as mentioned previously, this amounts to selecting a precise fraction of particles that are extracted from the thermal pool are and then accelerated to high energies. We thus accounted for the backreaction of the accelerated particles on the turbulence following the prescriptions of Lemoine et al. (2024). As detailed therein, backreaction becomes effective when the rate at which particles draw energy from the turbulent cascade exceeds the energy at which the cascade is replenished by energy injection. The former rate can be written ∼ν_accu_p (u_p is the suprathermal proton energy density), with ν_acc ∼ v_A²/(cℓ_c), while the latter reads ∼(v_A/ℓ_c) u_B (u_B is the magnetic energy density). Turbulent damping therefore becomes effective when u_p/u_B ∼ c/v_A. When we assume a plasma β_p unity (see Appendix A), this also means that damping takes place when $u_{\mathrm{p}}/p_{\mathrm{gas}}\sim c/v_{\mathrm{A}}\sim \mathcal{O}(1)$ here. When the particles start to damp the turbulence, the stochastic acceleration rate stalls, and acceleration becomes self-regulated.

To fix the parameters, we assumed that the cascade rate on the outer scale was 0.5 v_A/ℓ_c, as expected on general grounds. We used r_co = 30 r_g, v_adv = 0.02 c, and ℓ_c = 10 r_g as before, and we injected particles at γ₀ ∼ 1 as before, with an energy density lower by a factor of four than contained in the magnetized turbulence. We then tracked their acceleration using the two effective models (a) and (b) described by Lemoine et al. (2024), which mimic models (1) and (2) above. In detail, model (a) corresponds to the Fokker-Planck model and assumed v_A = 0.21 c, while model (b) describes advection in momentum space of a broken power-law spectrum assuming v_A = 0.28 c. Finally, we computed the global stationary spectrum over the corona by integrating the time-dependent spectra as in Eq. (3). The resulting spectrum is shown in the lower panel of Fig. 3. We stress that in this scenario, no ad hoc normalization of the energy content was adopted. This energy content is dictated by the backreaction of particles on the turbulence, which, as discussed above, implies that the total energy density in suprathermal particles must be about p_gasc/v_A. In units of p_gas, the energy density associated with the particle spectra shown in Fig. 3 is 3.4 in model (a) and 4.1 in model (b). The general shape of the spectrum does not depend much on the method used to carry out the integration. The key parameter that determines whether the overall flux can match the detailed spectrum at energies 30 − 300 TeV is the Alfvén velocity (more generally, ν_accτ_adv), which determines the location of the high-energy cutoff. The fact that the inferred spectrum can be reproduced at the price of reasonable parameters without an arbitrary normalization provides additional support and motivation for this scenario.

The global spectral shape of self-regulated acceleration can be understood as follows. In a first stage, particles are accelerated as in the test-particle picture of the first scenario that we investigated previously in this section. This occurs as long as the energy carried by the suprathermal particles remains below the threshold for backreaction, and it explains the similarity of the spectral slopes at low energies with that shown in the upper panel of Fig. 3. When feedback becomes efficient, acceleration stalls at the momentum determining the peak of the energy distribution, while particles with a higher momentum can continue to be accelerated because they interact with wavelength modes of larger scales that have not yet been damped. This tends to shape time-dependent energy spectra that are approximately flat at high energies. The average of these time-dependent spectra over the crossing time of the corona, accounting for energy and escape losses, then shapes the observed spectrum, similar to or slightly steeper than −2 in slope (per number). One phenomenological implication of this modeling is that in order to match a spectrum steeper than −2 in slope, the spectrum should be close to the cutoff region, that is, it should reveal increasing curvature toward higher energies. For higher acceleration rates or higher ν_accτ_adv, the maximum energy is pushed to higher values, and the spectrum in a given region therefore approaches −2 in slope.

We therefore recall the dependence of the cutoff location on ν_accτ_adv. However, if this combination of parameters were inhomogeneous throughout the corona, for example, because the Alfvénic velocity itself varied in space (as proposed in the context of X-ray binaries; see Nättilä 2024), the sharp high-energy cutoff would be replaced by a softer dependence associated with the tail of values of ν_accτ_adv. A more general interpretation of the proton spectrum is therefore called for, in which the flux level is determined by backreaction as modeled above, while the slope of the spectrum in the region probed by IceCube is dictated by the distribution of the acceleration rates (and advection times) experienced by the particle population. We defer the detailed investigation of this scenario to future work.

4. Conclusions

We have investigated scenarios of stochastic particle acceleration in a turbulent magnetized black hole corona and in the sheared velocity flow at the base of the jet to examine to which degree they can account for the proton energies inferred from IceCube neutrino observations of the Seyfert prototype NGC 1068, namely ϵ ∼ 30 − 300 TeV, and for the overall level of the energy density. To this end, we used recent results for stochastic particle acceleration in turbulence and paid particular attention to the transport equation, that is, we considered the influence of turbulent transport on escape losses and incorporated advection losses as a time limit on particle acceleration (Sect. 2.2). In Sect. 3 we compared three characteristic scenarios to the proton spectrum inferred in the case of NGC 1068 to conclude the following.

In a first scenario, which portrayed particle acceleration in the test-particle limit in a turbulent corona, we confirmed previous results that suggested that a characteristic Alfvénic velocity v_A ≃ 0.1 − 0.2 c for a coherence length ℓ_c ∼ 3 − 10 r_g can reproduce the general shape of the observed spectrum. However, this comes at the cost of an ad hoc normalization of the proton spectrum up to an appreciable fraction of the background plasma pressure, which in some cases exceeds unity. In terms of phenomenology, the spectrum is expected to reveal an increasingly strong curvature toward higher energies.

In a second scenario, we investigated the possibility that stochastic acceleration consumes less energy, but is completed by a phase of reacceleration in the sheared velocity flow at the base of a mildly relativistic jet. This flow, whether interpreted as an inner disk wind or jet wall, was envisioned to be launched from small radii well inside the corona; beyond this, however, it does not favor any specific model for the corona geometry. While this scenario admittedly introduces new parameters into the problem, it shows that shear reacceleration can play a role and lead to a satisfactory fit to the observed spectrum for a reasonable choice of parameters. In particular, depending on the conditions for escape, the particle distribution might extend to even higher energies, as shown in Fig. 3. In our context, the relative importance of shear acceleration depends on the magnitude of the shear, which determines the rate at which particles are re-energized in the jet boundary layer, just as it depends on the advection velocity in the corona. Similarly to the first scenario, the overall flux normalization was chosen freely and required that the energy density in high-energy protons lies close to the overall pressure of the background plasma.

This apparent coincidence is intriguing with respect to particle injection into the stochastic process, because it means that a specific fraction of particles must be extracted from the thermal pool that is smaller by some orders of magnitude than unity, yet such that after acceleration, it carries an amount of energy that is comparable to the pressure of their parent population. Based on the observation that the gas pressure is β_p times the magnetic energy density u_B, and that β_p ∼ 1, this observation may indicate that the acceleration process is as efficient as can be, that is, that the high-energy particles extract as much energy as they can from the turbulent cascade.

We therefore examined a third scenario in which turbulent acceleration was self-regulated by this feedback, which takes the form of damping of the turbulent cascade in the inertial range. In this scenario, the problem was fixed by the rate at which turbulent energy is injected on the outer scale ℓ_c, with a generic value of ∼u_Bv_A/ℓ_c. Adopting this rate and characteristic Alfvénic velocities as above, solving for the concomitant dynamical evolution of the turbulence and of the accelerated particles, we were able to reproduce the inferred proton spectra without having to specify the flux normalization. We did not explore the possibility that the accelerated particles might provide a backreaction at a more global (hydrodynamical) level on the structure of the flow itself, such as the liftoff of the hydrostatic equilibrium in the corona or jet deceleration in the case of shear acceleration. These issues clearly deserve further investigation.

This self-regulated acceleration process becomes operative provided the fraction χ₀ of particles that are extracted from the thermal pool and injected into stochastic acceleration verifies ${\chi }_0\gtrsim {(v_{\mathrm{A}}/c)}^{-1}{\beta }_{\mathrm{p}}^{-1}{({ϵ}_{\max }/{\overline{ϵ}}_{\mathrm{th}})}^{-1}$ in terms of the plasma β parameter (β_p), the Alfvénic velocity v_A (here playing the role of the characteristic eddy velocity on the outer scale of the turbulence), ϵ_max the maximum energy of the nonthermal energy spectrum, and ${\overline{ϵ}}_{\mathrm{th}}\sim k_{\mathrm{B}}T$ the mean thermal energy. If the injection fraction χ₀ is smaller than the above, feedback has not set in by the time the particles are accelerated to ϵ_max, and hence, particle acceleration proceeds as in a test-particle picture. However, if ${ϵ}_{\max }/{\overline{ϵ}}_{\mathrm{th}}\gg 1$, that is, if particles are able to reach high energies, feedback becomes natural, and the high-energy part of the accelerated spectrum then depends weakly on the details of the injection.

In this scenario of self-regulated stochastic acceleration, the source outputs as much energy as available in high-energy protons and achieves approximate equipartition between the supra-thermal proton and turbulent magnetic energy densities. This feature is all the more appealing in the context of high apparent neutrino luminosities such as were reported for NGC 1068. One direct consequence of this is that the overall neutrino luminosity of the source is not only controlled by the X-ray luminosity, which governs p − γ interactions, but also by the magnetic energy content. Our study indicates that other parameters are likely to influence the neutrino spectrum. We notably remarked that the maximum energy is strongly sensitive to the acceleration rate, for example, to the Alfvén velocity, the coherence length, and to the advection timescale through the corona. This dependence arises because the mean particle energy increases exponentially fast in turbulent acceleration up to a number of e-folds that is set by a combination of these parameters. An important corollary of this is that different sources might exhibit different spectra. An alternative way of viewing this strong dependence is that if the corona were inhomogeneous in terms of v_A (e.g.), the final spectrum is a sum of spectra obtained for differing v_A, or in other words, the final spectrum is shaped by the distribution of v_A (and other parameters) in the corona. Overall, it is therefore tempting to interpret the flux level as the consequence of self-regulated particle acceleration by turbulence damping (or possibly, viscous damping of the sheared velocity flow), and the general slope of the spectrum as given by the distribution of acceleration rates within the corona.

Finally, as pointed out to us by A. Beloborodov and A. Levinson, the possibility of turbulence damping introduces new issues. In particular, it may jeopardize the energization of electrons, which is expected to power the X-ray luminosity of the corona. We note, however, that in our description, the advection of particles through the corona effectively leads to stratification in terms of mean energy because time-dependent spectra also depend on the position. Close to the boundary where particles are injected into the corona, turbulence is therefore not affected, and electron energization is efficient. Furthermore, the damping of turbulence on scales of ≲l, for instance, does not prevent the formation of reconnection layers on scales ∼l, which offer ancillary energization channels. In an electron-ion plasma, reconnection can energize particles to a mean energy ∼(1 + σ)m_pc², where the magnetization parameter σ ≃ (v_A/c)² (at σ ≲ 1) denotes the ratio of the magnetic to plasma energy density (Werner et al. 2018). While turbulent reconnection is therefore not expected to affect the proton spectrum at the highest energies for moderate values of σ, it may contribute to electron energization. In a conservative model, it suffices to energize electrons to Lorentz factors γ_e < 10² to power the X-ray luminosity, and these electrons behave as mildly or subrelativistic protons. Overall, it appears reasonable to assume that the electrons can indeed be energized as needed in spite of turbulent damping. A detailed investigation of this question will form the basis of future work.

Acknowledgments

It is a pleasure to thank A. Beloborodov and A. Levinson, F. Oikonomou and E. Peretti as well as the members of the CN-6 connector of the Munich Excellence Cluster ORIGINS (in particular, P. Padovani, X. Rodrigues, and E. Resconi) for stimulating discussions. FMR acknowledges support by a DFG grant (RI 1187/8-1) and the kind hospitality of the IAP Paris.

References

Appendix A: Modelling the loss and escape rates

We incorporate energy and escape losses as follows. While the advection timescale τ_adv is momentum-independent, τ_esc ≃ r_co²/2κ varies with momentum through the spatial diffusion coefficient κ. We model transport as a combination of turbulent advection with an effective diffusion coefficient κ_turb = ℓ_cv_A/3, and scattering on magnetic inhomogeneities with a mean free path ${\lambda }_{\mathrm{scatt}}=r_{\mathrm{L}}^{1/3}{\ell }_{\mathrm{c}}^{2/3}$, see for example Berezinskii et al. (1990), or Lemoine (2023) and Kempski et al. (2023) for recent discussions. Accordingly, we set κ ≡ κ_turb + λ_scattc/3. As discussed in the main text, this indicates that particle escape is primarily governed by turbulent transport in the energy interval of phenomenological relevance, see Fig. A.1 for an illustration with our fiducial parameters ℓ_c = 10 r_g, v_A = 0.2 c, r_co = 30r_g and B = 10³ G. Note that, when computing shear acceleration, we disregard the effect of turbulent transport, since we assume that the turbulence in the sheared region is not highly dynamic.

Fig. A.1. Energy losses for: Ldisk = 5 × 1044 erg/s, LX = 0.8 × 1044 erg/s, ΓX = 1.95, ϵX,max = 128 keV, rco = 30 rg, rg = 2.35 × 1012M7.2 cm, B = 103 G, np = 2 × 109 cm−3. The diffusive escape time is shown for vA = 0.2 c, while the advection time assumes vadv = 0.02 c (as expected for the radial inflow velocity at rco with viscosity parameter α = 0.1).

To model radiative losses, we consider a disk emission represented by a multicolor black-body L_ϵ ∝ ϵ^4/3exp(−ϵ/ϵ_cd) for ϵ > 0.1ϵ_cd and L_ϵ ∝ ϵ³ for ϵ < 0.1ϵ_cd. The disk cutoff is ϵ_cd = 31.5 eV. The integrated disk luminosity is L_disk = 5 × 10⁴⁴ erg/s. The X-ray emission is a powerlaw L_ϵ ∝ ϵ^{1 − Γ_X}exp(−ϵ/ϵ_cX) for ϵ > ϵ_cd (i.e., above the disk cutoff energy), and L_ϵ ∝ ϵ² below. The X-ray cutoff is set at ϵ_cX = 128 keV and the integrated luminosity L_X = 0.8 × 10⁴⁴ erg/s.

We assumed a black hole of mass M_BH = 10^7.2 M_⊙, and corresponding gravitational radius r_g ≃ 2.35 × 10¹²M_7.2 cm (Schwarzschild radius r_s = 2r_g). The Eddington reference luminosity is L_Edd = 2 × 10⁴⁵M_7.2 erg/s. We note that for NGC 1068, the all-flavour, integrated (1.5-15 TeV) isotropic-equivalent neutrino luminosity is of order L_ν ≃ 10⁴²(d/10 Mpc)² erg/s ∼L_Edd/10³, and the estimated bolometric luminosity L_bol ∼ 5 × 10⁴⁴ erg/s (e.g., IceCube Collaboration 2022; Padovani et al. 2024). We take r_co = 30r_g for the equatorial extension of the corona. The thermal background proton number density is set to n_th = 2 × 10⁹ cm⁻³, corresponding to a plasma β parameter β_p ∼ 1 (we use 1 everywhere). Advective escape (radial inflow) in the corona is parameterized by v_adv, where we take v_add = α_dv_k(r_co) as the reference value, with v_k(r) the Keplerian velocity field and α_d ≃ 0.1 the viscosity parameter.

Figure A.1 shows the energy losses timescales for the given background photon spectral energy distribution. For comparison, for the chosen parameters, the required cosmic-ray power for a NGC 1068 type source is of order L_CR ∼ L_ν(ct_cool/r_co)∼(5 − 10)×10⁴³ erg/s, equivalent to an energy density of u_CR ≃ L_CR/(4πr_co²c)∼(3 − 6)×10⁴ erg/cm³ and comparable to the energy density u_B ≃ 4 × 10⁴ erg/cm³ of the magnetic field.

All Figures

Fig. 1. Proton energy spectra ϵ2 dncor/dϵ (per log-interval of energy) vs. proton energy ϵ predicted by stochastic acceleration in a turbulent corona, starting from mono-energetic protons with a Lorentz factor of γ0 ∼ 1. In each panel, the solid line corresponds to a solution obtained by integrating the Fokker-Planck equation up to time τadv = rco/vadv (model (1) including energy losses, advection, and escape as described in the text), and the dashed line shows the corresponding solution for the same parameters for the generalized Fermi model described by Eq. (5) (model (2)). In all panels, rco is set to 30rg, and in the top and middle panels, vadv ≃ 0.02 c. The upper panel shows the dependence on vA, assuming otherwise a coherence length scale for the turbulence ℓc = 10 rg and an advection velocity vadv = 0.02 c. From left to right, as indicated, vA/c = 0.05, vA/c = 0.10, vA/c = 0.20, and vA/c = 0.40. The middle panel examines the dependence on ℓc, assuming otherwise vA/c = 0.20 and vadv/c = 0.02. From right to left, as indicated, ℓc = 3 rg, ℓc = 10 rg, and ℓc = 30 rg. The bottom panel examines the influence of advection, characterized here by the advection velocity vadv (which sets the advection timescale τadv ≡ rco/vadv), assuming otherwise vA/c = 0.20 and ℓc = 3 rg. From right to left, as indicated, vadv/c = 0.02, vadv/c = 0.05, vadv/c = 0.16 and a (more extreme) case vadv/c = 0.49. The butterfly indicates the range of values needed to reproduce a neutrino spectrum as observed by IceCube for NGC 1068. The energy density spectrum is expressed in units of the background plasma pressure (see Sect. 2.1). In these units, each individual spectrum has been normalized to an integrated fraction of unity for the sake of illustration.

In the text

Fig. 2. Exemplary time-dependent proton energy distribution obtained from jet-shear acceleration with a Fokker-Planck model assuming continuous injection of seed protons with γ0 = 104. As a result of advective escape, the predicted particle spectrum will appear somewhat softer (red, t = τadv, z) compared to the otherwise quasi steady-state distribution (yellow). The gray shaded region illustrates the spectral characteristics required to reproduce the neutrino emission for NGC 1068 as measured by IceCube.

In the text

Fig. 3. Proton energy spectra dup/dlnϵ ≡ ϵ2 dn/dϵ vs. proton energy ϵ for three characteristic scenarios described in the text. For each panel, the energy density distribution is shown in units of the background plasma pressure. Top panel: Stochastic acceleration in the turbulent corona, with a characteristic Alfvénic velocity vA ≃ 0.2 c (see text for details), for model (1) and (2), corresponding to the solid and dashed lines, respectively. Middle panel: Combination of stochastic acceleration in the corona with shear acceleration at the base of the jet (see text for the details of the physical parameters). The stochastic contribution is indicated by dotted lines, the shear contribution by dashed lines, and the sum of the two by solid lines. The different colors correspond to models (1) and (2). For the top and middle panels, the normalization of the energy density in accelerated particles has been fixed in order to reproduce the inferred proton flux. Bottom panel: Modeling of stochastic acceleration, accounting for self-regulation by particle feedback on the turbulence (damping), for a characteristic Alfvén velocity vA ≃ 0.2 c. The plasma β parameter is 0.5 (see Appendix A).

In the text

	Fig. A.1. Energy losses for: Ldisk = 5 × 1044 erg/s, LX = 0.8 × 1044 erg/s, ΓX = 1.95, ϵX,max = 128 keV, rco = 30 rg, rg = 2.35 × 1012M7.2 cm, B = 103 G, np = 2 × 109 cm−3. The diffusive escape time is shown for vA = 0.2 c, while the advection time assumes vadv = 0.02 c (as expected for the radial inflow velocity at rco with viscosity parameter α = 0.1).
In the text