# Neutrinos from stochastic acceleration in black hole environments

Martin Lemoine, Frank Rieger

<figure><img src="image_1.png" /></figure>

## To cite this version:

Martin Lemoine, Frank Rieger. Neutrinos from stochastic acceleration in black hole environments. Astronomy & Astrophysics - A&A, 2025, 697, pp.A124. 10.1051/0004-6361/202453296. hal-04833890

HAL Id: hal-04833890

https://hal.science/hal-04833890v1

Submitted on 4 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Neutrinos from stochastic acceleration in black hole environments

Martin Lemoine¹ and Frank Rieger²,³

¹ Astroparticule et Cosmologie (APC), CNRS - Université Paris Cité, 75013 Paris, France
e-mail: mlemoine@apc.in2p3.fr

² Max Planck Institute for Plasma Physics (IPP), Boltzmannstraße 2, 85748 Garching, Germany

³ Institute for Theoretical Physics, Heidelberg University, Philosophenweg 12, 69120 Heidelberg, Germany
e-mail: frank.rieger@ipp.mpg.de

May 14, 2025

## ABSTRACT

Recent experimental results from the IceCube detector and their phenomenological interpretation suggest that the magnetized turbulent corona of nearby X-ray luminous Seyfert galaxies can produce ~ 1 – 10 TeV neutrinos via photo-hadronic interactions. We investigate the physics of stochastic acceleration in these environments in detail and examine the conditions under which the inferred proton spectrum can be explained. To this end, we used recent findings on particle acceleration in turbulence and paid particular attention to the transport equation, notably for transport in momentum space, turbulent transport outside of the corona, and advection through the corona. We first remark that the spectra we obtained are highly sensitive to the value of the acceleration rate, for instance, to the Alfvénic velocity. Then, we examined three prototype scenarios, one scenario of turbulent acceleration in the test-particle picture, another scenario in which particles were preaccelerated by turbulence and further energized by shear acceleration, and a final scenario in which we considered the effect of particle backreaction on the turbulence (damping), which self-regulates the acceleration process. We show that it is possible to obtain satisfactory fits to the inferred proton spectrum in all three cases, but we stress that in the first two scenarios, the energy content in suprathelial protons has to be fixed in an ad hoc manner to match the inferred spectrum at an energy density close to that contained in the turbulence. Interestingly, self-regulated acceleration by turbulence damping naturally brings the suprathelial particle energy content close to that of the turbulence and allowed us to reproduce the inferred flux level without additional fine-tuning. We also suggest that based on the strong sensitivity of the highest proton energy to the Alfvénic velocity (or acceleration rate), any variation in this quantity in the corona might affect (and in fact, set) the slope of the high-energy proton spectrum.

**Key words.** Acceleration of particles – Turbulence – Black hole physics – Magnetohydrodynamics (MHD)

## 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. 2019b,a), 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. 2020b; Winter & Lunardini 2023).

The recent announcement the arrival directions of (a subset of) high-energy ($\epsilon_v \ge 1 - 10$ TeV) neutrinos are correlated with local AGNs, in particular, from the nearby ($d \sim 10$ Mpc) Seyfert galaxy NGC 1068 (IceCube Collaboration et al. 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 $\epsilon_p \ge 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. 2020a). 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. 2020a; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fang et al. 2023; Fiorillo et al. 2024b; Murase et al. 2024; Mbarek et al. 2024; Fiorillo et al. 2024a; Das et al. 2024; Ambrósone 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\sigma$ 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 \times 10^3$ cm for a 1 GeV proton in a magnetic field with a strength of $B = 10^3$ G) from the macroscopic scale of the source (e.g., the black hole gravitational radius $r_g \equiv GM/c^2 \approx 1.5 \times 10^{12}$ M$_7$ cm)$^1$. In these environments, preacceleration in reconnecting current sheets is generic, and it can energize ions up to a Lorentz factor of about $\sigma$, where the magnetization parameter is defined as the ratio of the magnetic and plasma energy density (Werner et al. 2018). Unless $\sigma$ takes very high values (e.g., Fiorillo et al. 2024b; 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 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. 2020a; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fang et al. 2023; Fiorillo et al. 2024b; Murase et al. 2024; Mbarek et al. 2024; Fiorillo et al. 2024a; 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 quasispherical and compact, with a characteristic length of $r_{co} \le 100r_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_p \approx 1.2 \times 10^{11}$ (30/[^3] K, where[^3] $\hat{r} \equiv r/r_g$. While still high, radiative cooling results in electron temperatures that are significantly lower ($T_e \ll T_p$), as evidenced by the observed X-ray spectral cutoffs ($\epsilon_{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 $\beta$ parameters ($\beta_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 ($\tau_T < 1$; Zdziarski 1985; Stern et al. 1995, recent modeling suggested that in some sources, $\tau_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_A \approx 0.15$ (30/[^3] $)^{1/2} \beta_p^{-1/2} c$, which makes efficient stochastic acceleration feasible$^2$. In terms of $\beta_p$, the thermal plasma $\beta$ parameter, the corona magnetic field may reach strengths of $B = \sqrt{8\pi n_p m_p c^2 / (3[^3]\hat{r}\beta_p)} \approx 4 \times 10^3$ (30/[^3] $) M_7^{-1/2} \beta_p^{-1/2} \tau_T^{1/2}$ G, where $\tau_T = n_p \sigma_T r_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 \sim O(10$ G) on scales of $\sim 80r_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^3$ 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^4$ G, proton synchrotron losses can indeed be safely neglected. Furthermore, the rate of acceleration in stochastic acceleration is quantified by $v_A$ and $\ell_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 et al. 2022). This proton spectrum has a generic spectral index $s \sim -3$ in the energy interval $\sim 30-300$ TeV, and

<footnote>$^1$ We use the usual shorthand notation $Q_x = Q/10^x$, where $Q$ is expressed in cgs units, with the exception of solar masses for $M$.</footnote>

<footnote>$^2$ $v_A$ sets the characteristic velocity of turbulent eddies on the outer scale because we assumed a strongly turbulent corona with $\delta B \sim B$ and $\beta_p \sim 1$.</footnote>
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. 2020a; 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; Isliker et al. 2017; Pecora et al. 2018; Comisso & Sironi 2018, 2019; Kimura et al. 2019b; Zhdankin et al. 2019; 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; Lemoine & Malkov 2020; Demidem et al. 2020; Sioulas et al. 2020b,a; Lemoine 2021, 2022; 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_{\epsilon\epsilon}$, 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_{\epsilon\epsilon}$, in particular, the acceleration rate $v_{acc} \equiv D_{\epsilon\epsilon}/\epsilon^2$. The general transport equation that describes the evolution of the particle distribution $n_{\epsilon} \equiv dn/d\epsilon$ then takes the form

$$ \partial_t n_{\epsilon} = \mathcal{L}_{\text{stoch}} n_{\epsilon} - \partial_{\epsilon} (\dot{\epsilon}_{\text{loss}} n_{\epsilon}) - \frac{n_{\epsilon}}{\tau_{\text{esc}}}, \quad (1) $$

where $\mathcal{L}_{\text{stoch}}$ denotes the diffusive operator describing stochastic acceleration, that is,

$$ \mathcal{L}_{\text{stoch}} n_{\epsilon} \equiv \partial_{\epsilon} (D_{\epsilon\epsilon} \partial_{\epsilon} n_{\epsilon}) - 2 \partial_{\epsilon} (D_{\epsilon\epsilon} n_{\epsilon} / \epsilon) \quad (2) $$

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

The timescale $\tau_{\text{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 $\lambda_{\text{scatt}}$) and through turbulent transport, with a characteristic diffusion coefficient $\kappa_{\text{turb}} \sim \ell_c v_A / 3$. As the scattering mean free path is about $\lambda_{\text{scatt}} \sim r_L^{1/3} \ell_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 $\ell_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 $\kappa = \kappa_{\text{turb}} + \lambda_{\text{scatt}} c / 3$ and discuss the consequences in Sect. 3. The diffusive escape timescale is then written $\tau_{\text{esc}} \equiv r_{\text{co}}^2 / (2\kappa)$.

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 $\tau_{\text{adv}}$, which is defined in terms of the advection velocity $v_{\text{adv}}$ as $\tau_{\text{adv}} \equiv r_{\text{co}}/v_{\text{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' \in [t - \tau_{\text{adv}}, t]$ and accelerated for a duration $\tau \equiv t - t'$. This spectrum can be expressed as

$$ \frac{dn_{\text{cor}}}{d\epsilon} = \int_{0}^{\tau_{\text{adv}}} d\tau \int_{0}^{+\infty} d\epsilon' G_{\text{turb}}(\epsilon; \epsilon', \tau) \frac{d\dot{n}_{\text{inj}}}{d\epsilon'} \quad (3) $$

in terms of the Green function $G_{\text{turb}}(\epsilon; \epsilon', \tau)$, which is associated with Eq. (1). This determines the probability of a transition from $\epsilon'$ to $\epsilon$ over a time interval $\tau$, and of the injection rate spectrum $d\dot{n}_{\text{inj}}/d\epsilon$. 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_{\text{cor}}/d\epsilon$ can be regarded as the stationary average spectrum. In the former case, the time-dependent spectrum at location $x$ can be written as $dn_{\text{turb}}(x)/d\epsilon = \int d\epsilon' G_{\text{turb}}(\epsilon; \epsilon', x/v_{\text{adv}})dn_{\text{inj}}/d\epsilon'$ in terms of the injection density $dn_{\text{inj}}/d\epsilon'$ at the initial boundary, and the substitutions $x \rightarrow v_{\text{adv}}\tau$, $dn_{\text{inj}}/d\epsilon' \rightarrow \tau_{\text{adv}}d\dot{n}_{\text{inj}}/d\epsilon'$ 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\dot{n}_{\text{inj}}/d\epsilon'$ (independent of $x$), so that

$$ \frac{dn_{\text{turb}}(x)}{d\epsilon} = \int_{0}^{x} \frac{dx'}{v_{\text{adv}}} \int_{0}^{+\infty} d\epsilon' G_{\text{turb}}(\epsilon; \epsilon', (x-x')/v_{\text{adv}}) \frac{d\dot{n}_{\text{inj}}}{d\epsilon'} \quad (4) $$

Then, $dn_{\text{cor}}/d\epsilon$ is again recovered by taking the average of $dn_{\text{turb}}(x)/d\epsilon$ over $x$ throughout the corona.

The advection term is sometimes included in the transport equation as an escape term of the form $-n_{\epsilon}/\tau_{\text{adv}}$, which implicitly implies that particles can spend a time $> \tau_{\text{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 $v_{\text{acc}}\tau_{\text{adv}} \sim 1$ (or lower, if other losses dominate).

Recent numerical experiments of stochastic particle acceleration, which either tracked test particles in a magnetohydrodynamical (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_{\epsilon\epsilon} \approx 0.2\gamma^2(v_A/c)^2c/\ell_c$ for particles with a Lorentz factor $\gamma \equiv \epsilon/mc^2$, where $v_A$ denotes the Alfvén velocity in the total magnetic field, and $\ell_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 \text{ c}$, $r_{co} = 30 \text{ r}_g$ and $\ell_c = 10 \text{ r}_g$, the advection timescale reads $\tau_{adv} \approx 150 \ell_c / c$, while PIC simulations typically run over $\approx 10 \ell_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

$$ \mathcal{L}_{\text{stoch}} n_{\epsilon} \equiv \int_{0}^{+\infty} d\epsilon' \left[ \frac{\varphi(\epsilon|\epsilon')}{t_{\epsilon'}} n_{\epsilon'}(t) - \frac{\varphi(\epsilon'|\epsilon)}{t_{\epsilon}} n_{\epsilon}(t) \right], \quad (5) $$

where $\varphi(\epsilon|\epsilon')$ represents the kernel describing the probability of jumping from $\epsilon'$ to $\epsilon$ over a time interval $t_{\epsilon'}$, 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 \text{ c}$. It can nevertheless be extrapolated to other Alfvénic velocities by rescaling its width in proportion to $(v_A/0.4 \text{ c})^2$, 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., Kardashiev 1962) that follows from a purely diffusive Fokker-Planck kernel with $D_{\epsilon\epsilon} \propto \epsilon^2$ (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 $\gamma_0$ (Schlickeiser 1984), up to a maximum energy that is either dictated by the finite time ($\tau_{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 $\gamma_0 \sim 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_{\epsilon\epsilon}$ or $\varphi(\epsilon|\epsilon')$ for the acceleration rate. They both depend on $v_A$, $\epsilon$ and $\ell_c$; $\tau_{adv} \equiv r_{co}/v_{adv}$, which sets the integration time before advective escape at velocity $v_{adv}$, removes the particles $\tau_{esc}$ and $\tau_{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), $\ell_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.

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_{\epsilon\epsilon} \propto \epsilon^2$, the mean energy increases exponentially fast in time, that is, $\langle\epsilon\rangle \propto \exp(4v_{acc}t)$, with $v_{acc} \equiv D_{\epsilon\epsilon}/\epsilon^2$ the acceleration rate defined earlier. Consequently, in a fixed time interval $t = \tau_{adv}$ over which $v_{acc}t \gg 1$ (needed to achieve acceleration to high energies), changing $v_{acc}$ by a factor of a few changes the maximum energy by orders of magnitude. The dependence of $v_{acc}$ on the parameters is recalled above, $v_{acc} \propto v_A^2/\ell_c$. We assumed $\delta B/B \gtrsim 1$, otherwise $v_{acc}$ should be further suppressed by a factor $(\delta B/B)^2$. 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 $v_{acc}\tau_{adv} \propto 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 $\lesssim \text{GeV}$ to $\sim 0.1 \text{ 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 $\sim 10^5$, its pressure is comparable to that of its parent population. If more particles had been
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Proton energy spectra $\epsilon^2 dn_{\text{cor}}/d\epsilon$ (per log-interval of energy) vs. proton energy $\epsilon$ predicted by stochastic acceleration in a turbulent corona, starting from mono-energetic protons with a Lorentz factor of $\gamma_0 \sim 1$. In each panel, the solid line corresponds to a solution obtained by integrating the Fokker-Planck equation up to time $\tau_{\text{adv}} = r_{\text{co}}/v_{\text{adv}}$ [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, $r_{\text{co}}$ is set to $30r_g$, and in the top and middle panels, $v_{\text{adv}} \simeq 0.02c$. The upper panel shows the dependence on $v_A$, assuming otherwise a coherence length scale for the turbulence $\ell_c = 10r_g$ and an advection velocity $v_{\text{adv}} = 0.02c$. From left to right, as indicated, $v_A/c = 0.05$, $v_A/c = 0.10$, $v_A/c = 0.20$, and $v_A/c = 0.40$. The middle panel examines the dependence on $\ell_c$, assuming otherwise $v_A/c = 0.20$ and $v_{\text{adv}}/c = 0.02$. From right to left, as indicated, $\ell_c = 3r_g$, $\ell_c = 10r_g$, and $\ell_c = 30r_g$. The bottom panel examines the influence of advection, characterized here by the advection velocity $v_{\text{adv}}$ (which sets the advection timescale $\tau_{\text{adv}} \equiv r_{\text{co}}/v_{\text{adv}}$), assuming otherwise $v_A/c = 0.20$ and $\ell_c = 3r_g$. From right to left, as indicated, $v_{\text{adv}}/c = 0.02$, $v_{\text{adv}}/c = 0.05$, $v_{\text{adv}}/c = 0.16$ and a (more extreme) case $v_{\text{adv}}/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.</figcaption></figure>

injected, the energy requirement would have become excessive, while in the opposite case, the proton flux falls short of the re-

quired 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 ($\le 0.1c$) on scales of several dozen parsec ($\sim 10^7 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_{\text{jet}} \sim 10^{43}$ 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^{42}$ 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 $\gamma_p \le 10^4$ (10 TeV) either due to $v_A \le 0.1c$ or $\ell_c \gg 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 $\gamma_j \equiv 1/(1-u_z^2/c^2)^{1/2} \le 3$, characterized by a lateral velocity shear $u_z(r)$ with a width of $\Delta r = \hat{r}_{\text{sh}}r_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
<figure><img src="image_4.png" /><figcaption>**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 $\gamma_0 = 10^4$. As a result of advective escape, the predicted particle spectrum will appear somewhat softer (red, $t = \tau_{\text{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.</figcaption></figure>

diffusive Fokker-Planck equation for $n(\epsilon)$ of the type Eq. (1), with mean momentum diffusion coefficient $D_{\epsilon\epsilon} \propto p^2 \bar{D}_p$, where for energetic protons, $\epsilon = pc = \gamma m_p c^2$. The coefficient $\bar{D}_p \equiv a_g \lambda/c$ can be obtained from the local (r-dependent) shear coefficient $D_p(r) = (1/15)\gamma_j(r)^4(\partial u_z/\partial r)^2\lambda/c$ by suitable spatial averaging, for example, $a_g = (2/3)(c/\Delta r)^2/w$ with $w = 116 (\ln[(1+\beta_0)/(1-\beta_0)])^{-2}$ in the case of a linearly decreasing shear flow profile with a maximum (spine) speed $\beta_0 = (1-1/\gamma_{j,max}^2)^{1/2}$ (Rieger & Duffy 2019, 2022). In the former expression, $\lambda = c\tau_c = \eta(\Delta r)(r_L/\Delta r)^{1/3}$ denotes the scattering mean free path ($r_L$ the gyroradius), for which we took a scaling with $\eta = 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, $\tau_{\text{esc}} = (\Delta r)^2/(2\kappa)$ with $\kappa = \lambda c/3$, or through advective transport (along $z$), $\tau_{\text{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_{\text{adv,z}} \approx z/\langle u_z \rangle = \bar{z} r_{\text{co}}/\langle u_z \rangle$, with $\langle u_z \rangle$ the mean velocity (corresponding to $\beta_0 c/2$), and $\bar{z} \equiv z/r_{\text{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 $\tau_c(\gamma_0)$) assuming that seed protons are continuously injected with a Lorentz factor of $\gamma_0 = 10^4$ into a shearing flow with $\gamma_j = 3$, $B = 50$ G, $\hat{r}_{\text{sh}} = 2$ and $\bar{z} = 2$, assuming a corona size of $r_{\text{co}} = 30 r_g$.

In general, the predicted spectra are sensitive to the choice of parameters, for example, the time available for particle acceleration (set by $\tau_{\text{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, $\gamma_j \le 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 might 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. 2020a; Inoue et al. 2020; Kheirandish et al. 2021; Murase 2022; Eichmann et al. 2022; Fiorillo et al. 2024a); 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 \sim 16$ Mpc), for instance, the characteristic neutrino energy seems higher by a factor $\sim 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_{\text{co}} = 30 r_g$, $\ell_c = 10 r_g$, $v_{\text{adv}} = 0.02 c$ and $v_A \approx 0.2 c$, their flux normalizations are different. The ratio of integrated energy density of suprathermal particles to background plasma pressure is $\approx 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 co-
herence length scale. However, the general spectral shape should be preserved.

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_{\text{shear}}(\epsilon; \epsilon', \tau_{\text{shear}})$, that is, the probability of jumping from $\epsilon'$ to $\epsilon$ over the timescale $\tau_{\text{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_{\text{turb}}/d\epsilon$, we write the stationary spectrum in the shear region as

$$ \frac{dn_{\text{shear}}}{d\epsilon} = \frac{1}{\tau_{\text{adv,z}}} \int_{0}^{\tau_{\text{adv,z}}} d\tau \int d\epsilon' G_{\text{shear}}(\epsilon; \epsilon', \tau) \frac{v_{\text{adv}}}{v_{\text{shear}}} \frac{dn_{\text{turb}}}{d\epsilon'}, \quad (6) $$

with $\tau_{\text{adv,z}}$ the advection time in the sheared region, as before, and

$$ \frac{dn_{\text{turb}}}{d\epsilon'} = \int d\epsilon'' G_{\text{turb}}(\epsilon'; \epsilon'', \tau_{\text{adv}}) \frac{dn_{\text{inj}}}{d\epsilon''}, \quad (7) $$

where $dn_{\text{inj}}/d\epsilon''$ characterizes (Sect. 2.2) the injection spectrum at the entry of the turbulent region, as before. We introduced a ratio $v_{\text{adv}}/v_{\text{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_{\text{co}}$ and took $\langle u_z \rangle = \beta_0 c/2 \approx 0.5 c$.

Finally, the total stationary spectrum was obtained as the sum of the spectra in each region, $dn_{\text{cor}}/d\epsilon$ [Eq. (3)] and $dn_{\text{shear}}/d\epsilon$ [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 $\propto \epsilon^{-3}$ (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_{\text{adv}} = 0.02 c$, $\ell_c = 10 r_g$, and $r_{\text{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,

<figure><img src="image_5.png" /><figcaption>Fig. 3. Proton energy spectra $du_p/d\ln\epsilon \equiv \epsilon^2 dn/d\epsilon$ vs. proton energy $\epsilon$ 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 $v_A \approx 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 $v_A \approx 0.2 c$. The plasma $\beta$ parameter is 0.5 (see Appendix A).</figcaption></figure>

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_{\text{cor}}/d\epsilon$ 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 ~ $v_{\text{acc}}u_p$ ($u_p$ is the suprathermal proton energy density), with $v_{\text{acc}} \sim v_A^2/(c\ell_c)$, while the latter reads ~ ($v_A/\ell_c$) $u_B$ ($u_B$ is the magnetic energy density). Turbulent damping therefore becomes effective when $u_p/u_B \sim c/v_A$. When we assume a plasma $\beta_p$ unity (see App. A), this also means that damping takes place when $u_p/p_{\text{gas}} \sim c/v_A \sim 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/\ell_c$, as expected on general grounds. We used $r_{\text{co}} = 30 r_g$, $v_{\text{adv}} = 0.02 c$, and $\ell_c = 10 r_g$ as before, and we injected particles at $\gamma_0 \sim 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_{\text{gas}}c/v_A$. In units of $p_{\text{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, $v_{\text{acc}}\tau_{\text{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 $v_{\text{acc}}\tau_{\text{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 $v_{\text{acc}}\tau_{\text{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ätilä 2024), the sharp high-energy cutoff would be replaced by a softer dependence associated with the tail of values of $v_{\text{acc}}\tau_{\text{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 $\epsilon \sim 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 \approx 0.1-0.2 c$ for a coherence length $\ell_c \sim 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 \sim 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 $l_c$, with a generic value of $\sim u_B v_A / l_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 $\chi_0$ of particles that are extracted from the thermal pool and injected into stochastic acceleration verifies $\chi_0 \gtrsim (v_A/c)^{-1} \beta_p^{-1} (\epsilon_{max}/\epsilon_{th})^{-1}$ in terms of the plasma $\beta$ parameter ($\beta_p$), the Alfvénic velocity $v_A$ (here playing the role of the characteristic eddy velocity on the outer scale of the turbulence), $\epsilon_{max}$ the maximum energy of the nonthermal energy spectrum, and $\bar{\epsilon}_{th} \sim k_B T$ the mean thermal energy. If the injection fraction $\chi_0$ is smaller than the above, feedback has not set in by the time the particles are accelerated to $\epsilon_{max}$, and hence, particle acceleration proceeds as in a test-particle picture. However, if $\epsilon_{max}/\bar{\epsilon}_{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 suprathelial 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-\gamma$ 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 $\lesssim l$, for instance, does not prevent the formation of reconnection layers on scales $\sim l$, which offer ancillary energization channels. In an electron-ion plasma, reconnection can energize particles to a mean energy $\sim (1+\sigma)m_p c^2$, where the magnetization parameter $\sigma \approx (v_A/c)^2$ (at $\sigma \lesssim 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 $\sigma$, it may contribute to electron energization. In a conservative model, it suffices to energize electrons to Lorentz factors $\gamma_e < 10^2$ 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.

*Acknowledgements.* 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

Abbasi, R., Ackermann, M., Adams, J., et al. 2025, ApJ, 981, 131

Abbasi, R., Ackermann, M., Adams, J., et al. 2024, arXiv e-prints, arXiv:2406.07601

Aharonian, F. & Neronov, A. 2005, ApJ, 619, 306

Aharonian, F. A., Belyanin, A. A., Derishev, E. V., Kocharovsky, V. V., & Kocharovsky, V. V. 2002, Phys. Rev. D, 66, 023005

Ambrosone, A. 2024, J. Cosmology Astropart. Phys., 2024, 075

Beloborodov, A. M. 1999, MNRAS, 305, 181

Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S. 1990, Astrophysics of cosmic rays (Amsterdam: North-Holland, 1990, edited by Ginzburg, V.L.)

Berezinsky, V. S. & Ginzburg, V. L. 1981, in International Cosmic Ray Conference, Vol. 1, International Cosmic Ray Conference, 238

Bresci, V., Lemoine, M., Gremillet, L., et al. 2022, Phys. Rev. D, 106, 023028

Brunetti, G. & Lazarian, A. 2016, Mon. Not. R. Astron. Soc., 458, 2584

Comisso, L. & Sironi, L. 2018, Phys. Rev. Lett., 121, 255101

Comisso, L. & Sironi, L. 2019, Astrophys. J., 886, 122

Comisso, L. & Sironi, L. 2022, ApJ, 936, L27

Das, A., Zhang, B. T., & Murase, K. 2024, ApJ, 972, 44

De Villiers, J.-P., Hawley, J. F., & Krolik, J. H. 2003, ApJ, 599, 1238

De Villiers, J.-P., Hawley, J. F., Krolik, J. H., & Hirose, S. 2005, ApJ, 620, 878
Demidem, C., Lemoine, M., & Casse, F. 2020, Phys. Rev. D, 102, 023003

Dermer, C. D., Miller, J. A., & Li, H. 1996, ApJ, 456, 106

Dexter, J. & Begelman, M. C. 2024, MNRAS, 528, L157

Di Matteo, T., Blackman, E. G., & Fabian, A. C. 1997, MNRAS, 291, L23

Eichler, D. 1979a, ApJ, 232, 106

Eichler, D. 1979b, Astrophys. J., 229, 413

Eichmann, B., Oikonomou, F., Salvatore, S., Dettmar, R.-J., & Tjus, J. B. 2022, ApJ, 939, 43

Eilek, J. A. 1979, Astrophys. J., 230, 373

Fabian, A. C., Lohfink, A., Kara, E., et al. 2015, MNRAS, 451, 4375

Fang, K., Lopez Rodriguez, E., Halzen, F., & Gallagher, J. S. 2023, ApJ, 956, 8

Farrar, G. R. & Piran, T. 2014, arXiv e-prints, arXiv:1411.0704

Fiorillo, D. F. G., Comisso, L., Peretti, E., Petropoulou, M., & Sironi, L. 2024a, ApJ, 974, 75

Fiorillo, D. F. G., Petropoulou, M., Comisso, L., Peretti, E., & Sironi, L. 2024b, ApJ, 961, L14

Gallimore, J. F., Baum, S. A., O'Dea, C. P., & Pedlar, A. 1996, ApJ, 458, 136

Guépin, C., Kotera, K., Barausse, E., Fang, K., & Murase, K. 2018, A&A, 616, A179

Halzen, F. 2023, arXiv e-prints, arXiv:2305.07086

IceCube Collaboration, Abbasi, R., Ackermann, M., et al. 2022, Science, 378, 538

Inoue, Y. & Doi, A. 2018, ApJ, 869, 114

Inoue, Y., Khangulyan, D., & Doi, A. 2020, ApJ, 891, L33

Inoue, Y., Takasao, S., & Khangulyan, D. 2024, PASJ, 76, 996

Isliker, H., Vlahos, L., & Constantinescu, D. 2017, Phys. Rev. Lett., 119, 045101

Kakuwa, J. 2016, Astrophys. J., 816, 24

Kammoun, E., Lohfink, A. M., Masterson, M., et al. 2024, Frontiers in Astronomy and Space Sciences, 10, 1308056

Kamraj, N., Brightman, M., Harrison, F. A., et al. 2022, ApJ, 927, 42

Kamraj, N., Harrison, F. A., Baloković, M., Lohfink, A., & Brightman, M. 2018, ApJ, 866, 124

Kara, E., García, J. A., Lohfink, A., et al. 2017, MNRAS, 468, 3489

Kardashev, N. S. 1962, Soviet Ast., 6, 317

Kempski, P., Fielding, D. B., Quataert, E., et al. 2023, MNRAS, 525, 4985

Kheirandish, A., Murase, K., & Kimura, S. S. 2021, ApJ, 922, 45

Kimura, S. S., Murase, K., & Mészáros, P. 2019a, Phys. Rev. D, 100, 083014

Kimura, S. S., Murase, K., & Zhang, B. T. 2018, Phys. Rev. D, 97, 023026

Kimura, S. S., Tomida, K., & Murase, K. 2019b, MNRAS, 485, 163

King, A. L., Lohfink, A., & Kara, E. 2017, ApJ, 835, 226

Lemoine, M. 2019, Phys. Rev. D, 99, 083006

Lemoine, M. 2021, Phys. Rev. D, 104, 063020

Lemoine, M. 2022, Phys. Rev. Lett., 129, 215101

Lemoine, M. 2023, Journal of Plasma Physics, 89, 175890501

Lemoine, M. & Malkov, M. A. 2020, Mon. Not. R. Astron. Soc., 499, 4972

Lemoine, M., Murase, K., & Rieger, F. 2024, Phys. Rev. D, 109, 063006

Lenain, J. P., Ricci, C., Türler, M., Dorner, D., & Walter, R. 2010, A&A, 524, A72

Levinson, A. & Boldt, E. 2002, Astroparticle Physics, 16, 265

Liska, M. T. P., Musoke, G., Tchekhovskoy, A., Porth, O., & Beloborodov, A. M. 2022, ApJ, 935, L1

Litvinenko, Y. E. 2012, Astron. Astrophys., 544, A94

Liu, B. F., Mineshige, S., & Shibata, K. 2002, ApJ, 572, L173

Lynn, J. W., Quataert, E., Chandran, B. D. G., & Parrish, I. J. 2014, Astrophys. J., 791, 71

Malzac, J., Beloborodov, A. M., & Poutanen, J. 2001, MNRAS, 326, 417

May, D. & Steiner, J. E. 2017, MNRAS, 469, 994

Mbarek, R., Philippov, A., Chernoglazov, A., Levinson, A., & Mushotzky, R. 2024, Phys. Rev. D, 109, L101306

Merloni, A. & Fabian, A. C. 2001, MNRAS, 321, 549

Merloni, A. & Fabian, A. C. 2002, MNRAS, 332, 165

Murase, K. 2022, ApJ, 941, L17

Murase, K., Karwin, C. M., Kimura, S. S., Ajello, M., & Buson, S. 2024, ApJ, 961, L34

Murase, K., Kimura, S. S., & Mészáros, P. 2020a, Phys. Rev. Lett., 125, 011101

Murase, K., Kimura, S. S., Zhang, B. T., Oikonomou, F., & Petropoulou, M. 2020b, ApJ, 902, 108

Nätilä, J. 2024, Nature Communications, 15, 7026

Neronov, A., Savchenko, D., & Semikoz, D. V. 2024, Phys. Rev. Lett., 132, 101002

Padovani, P., Resconi, E., Ajello, M., et al. 2024, Nature Astronomy, 8, 1077

Pecora, F., Servidio, S., Greco, A., et al. 2018, J. Plasma Phys., 84, 725840601

Pezzi, O., Blasi, P., & Matthaeus, W. H. 2022, ApJ, 928, 25

Poutanen, J., Veledina, A., & Beloborodov, A. M. 2023, ApJ, 949, L10

Protheroe, R. J. & Szabo, A. P. 1992, Phys. Rev. Lett., 69, 2885

Pugliese, F., Brodiano, M., Andrés, N., & Dmitruk, P. 2023, ApJ, 959, 28

Rieger, F. M. 2019, Galaxies, 7, 78

Rieger, F. M. & Aharonian, F. A. 2009, A&A, 506, L41

Rieger, F. M. & Duffy, P. 2019, ApJ, 886, L26

Rieger, F. M. & Duffy, P. 2022, ApJ, 933, 149

Rieger, F. M. & Mannheim, K. 2000, A&A, 353, 473

Rieger, F. M. & Mannheim, K. 2002, A&A, 396, 833

Roy, A. L., Wilson, A. S., Ulvestad, J. S., & Colbert, J. M. 2000, in EVN Symposium 2000, Proceedings of the 5th european VLBI Network Symposium, ed. J. E. Conway, A. G. Polatidis, R. S. Booth, & Y. M. Pihlström, 7

Salvatore, S., Eichmann, B., Rodrigues, X., Dettmar, R. J., & Becker Tjus, J. 2024, A&A, 687, A139

Schlickeiser, R. 1984, Astron. Astrophys., 136, 227

Schlickeiser, R. 1989, ApJ, 336, 243

Sioulas, N., Isliker, H., & Vlahos, L. 2020a, Astrophys. J., 895, L14

Sioulas, N., Isliker, H., Vlahos, L., Koumtzis, A., & Pisokas, T. 2020b, Mon. Not. R. Astron. Soc., 491, 3860

Sironi, L. & Beloborodov, A. M. 2020, ApJ, 899, 52

Sridhar, N., Ripperda, B., Sironi, L., Davelaar, J., & Beloborodov, A. M. 2025, ApJ, 979, 199

Stern, B. E., Poutanen, J., Svensson, R., Sikora, M., & Begelman, M. C. 1995, ApJ, 449, L13

Trotta, D., Franci, L., Burgess, D., & Hellinger, P. 2020, ApJ, 894, 136

Wang, J., Rieger, F. M., & Mizuno, Y. 2024, ApJ, 967, 36

Wang, J.-S., Reville, B., Mizuno, Y., Rieger, F. M., & Aharonian, F. A. 2023, MNRAS, 519, 1872

Waxman, E. 2024, Talk given at RICAP 24

Webb, G. M., Xu, Y., Biermann, P. L., et al. 2023, ApJ, 958, 169

Werner, G. R., Uzdensky, D. A., Begelman, M. C., Cerutti, B., & Nalewajko, K. 2018, MNRAS, 473, 4840

Winter, W. & Lunardini, C. 2023, ApJ, 948, 42

Wong, K., Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2020, Astrophys. J., 893, L7

Xu, S. & Lazarian, A. 2023, ApJ, 942, 21

Zdziarski, A. A. 1985, ApJ, 289, 514

Zhdankin, V., Uzdensky, D. A., Werner, G. R., & Begelman, M. C. 2019, Phys. Rev. Lett., 122, 055101

Zhdankin, V., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2017, Phys. Rev. Lett., 118, 055103
## Appendix A: Modelling the loss and escape rates

<figure><img src="image_6.png" /><figcaption>**Fig. A.1.** Energy losses for: $L_{\text{disk}} = 5 \times 10^{44}$ erg/s, $L_X = 0.8 \times 10^{44}$ erg/s, $\Gamma_X = 1.95$, $\epsilon_{X,max} = 128$ keV, $r_{\text{co}} = 30 r_g$, $r_g = 2.35 \times 10^{12} M_{7.2}$ cm, $B = 10^3$ G, $n_p = 2 \times 10^9$ cm$^{-3}$. The diffusive escape time is shown for $v_A = 0.2$ c, while the advection time assumes $v_{\text{adv}} = 0.02$ c (as expected for the radial inflow velocity at $r_{\text{co}}$ with viscosity parameter $\alpha = 0.1$).</figcaption></figure>

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_{\text{CR}} \sim L_V(ct_{\text{cool}}/r_{\text{co}}) \sim (5-10) \times 10^{43}$ erg/s, equivalent to an energy density of $u_{\text{CR}} \approx L_{\text{CR}}/(4\pi r_{\text{co}}^2 c) \sim (3-6) \times 10^4$ erg/cm$^3$ and comparable to the energy density $u_B \approx 4 \times 10^4$ erg/cm$^3$ of the magnetic field.

We incorporate energy and escape losses as follows. While the advection timescale $\tau_{\text{adv}}$ is momentum-independent, $\tau_{\text{esc}} \approx r_{\text{co}}^2/2\kappa$ varies with momentum through the spatial diffusion coefficient $\kappa$. We model transport as a combination of turbulent advection with an effective diffusion coefficient $\kappa_{\text{turb}} = \ell_c v_A/3$, and scattering on magnetic inhomogeneities with a mean free path $\lambda_{\text{scatt}} = r_L^{1/3} \ell_c^{2/3}$, see for example Berezinskii et al. (1990), or Lemoine (2023) and Kempski et al. (2023) for recent discussions. Accordingly, we set $\kappa \equiv \kappa_{\text{turb}} + \lambda_{\text{scatt}} c/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 $\ell_c = 10 r_g$, $v_A = 0.2$ c, $r_{\text{co}} = 30 r_g$ and $B = 10^3$ 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.

To model radiative losses, we consider a disk emission represented by a multicolor black-body $L_{\epsilon} \propto \epsilon^{4/3} \exp(-\epsilon/\epsilon_{cd})$ for $\epsilon > 0.1\epsilon_{cd}$ and $L_{\epsilon} \propto \epsilon^3$ for $\epsilon < 0.1\epsilon_{cd}$. The disk cutoff is $\epsilon_{cd} = 31.5$ eV. The integrated disk luminosity is $L_{\text{disk}} = 5 \times 10^{44}$ erg/s. The X-ray emission is a powerlaw $L_{\epsilon} \propto \epsilon^{1-\Gamma} \exp(-\epsilon/\epsilon_{cX})$ for $\epsilon > \epsilon_{cd}$ (i.e., above the disk cutoff energy), and $L_{\epsilon} \propto \epsilon^2$ below. The X-ray cutoff is set at $\epsilon_{cX} = 128$ keV and the integrated luminosity $L_X = 0.8 \times 10^{44}$ erg/s.

We assumed a black hole of mass $M_{\text{BH}} = 10^{7.2} M_{\odot}$, and corresponding gravitational radius $r_g \approx 2.35 \times 10^{12} M_{7.2}$ cm (Schwarzschild radius $r_s = 2r_g$). The Eddington reference luminosity is $L_{\text{Edd}} = 2 \times 10^{45} 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_V \approx 10^{42} (d/10 \text{ Mpc})^2$ erg/s $\sim L_{\text{Edd}}/10^3$, and the estimated bolometric luminosity $L_{\text{bol}} \sim 5 \times 10^{44}$ erg/s (e.g., IceCube Collaboration et al. 2022; Padovani et al. 2024). We take $r_{\text{co}} = 30 r_g$ for the equatorial extension of the corona. The thermal background proton number density is set to $n_{\text{th}} = 2 \times 10^9 \text{ cm}^{-3}$, corresponding to a plasma $\beta$ parameter $\beta_p \sim 1$ (we use 1 everywhere). Advective escape (radial inflow) in the corona is parameterized by $v_{\text{adv}}$, where we take $v_{\text{add}} = \alpha_d v_k(r_{\text{co}})$ as the reference value, with $v_k(r)$ the Keplerian velocity field and $\alpha_d \approx 0.1$ the viscosity parameter.