1 INTRODUCTION↩︎

Cosmic rays (CRs) are high-energy charged particles propagating through interstellar and intergalactic space, spanning energies from \(10^{9}\,\mathrm{eV}\) to \(10^{20}\,\mathrm{eV}\). Understanding their origin is not only a fundamental problem in high-energy astrophysics but also essential for galaxy evolution. Despite extensive studies, the origin of CRs across different energy ranges remains unresolved. Supernova remnants (SNRs) are regarded as the most promising candidates up to the knee energy (a few \(\times10^{15}\,\mathrm{eV}\)) [1]–[3], whereas the sources and acceleration mechanisms at and above the knee energy are largely unknown.

Active galactic nuclei (AGN) are among the leading candidates for the origin of high-energy CRs across and beyond the knee energy. Proposed acceleration sites in AGN include relativistic jets [4]–[8], coronae [9]–[13], accretion flows [14], [15], and outflows [16]–[20]. In addition, AGN are considered promising multi-messenger sources because the acceleration of CRs is expected to produce high-energy neutrinos and gamma rays [17]–[19], [21]–[25].

Among the AGN environments, ultrafast outflows (UFOs) have recently attracted growing attention. UFOs are mildly relativistic winds, launched at tens of percent of the speed of light, and identified through blueshifted Fe K-shell absorption features (Fe XXV He\(\alpha\) at 6.7 keV and Fe XXVI Ly\(\alpha\) at 7.0 keV) or P-Cygni–like profiles. They are observed in several tens of percent of AGN [26], [27], suggesting that they are more common than powerful jets. Observationally, the detection of GeV gamma-rays in some AGN—often difficult to explain by AGN coronae alone—provides further motivation to consider UFOs as promising acceleration sites [25], [28]. In fact, GeV gamma-ray emission has also been reported in several UFOs in AGN [22].

UFOs colliding with ambient gas are expected to form shock structures. These shocks have been proposed as promising sites for diffusive shock acceleration (DSA) [2], [3], [29]–[31], in which CRs repeatedly cross the shock front via pitch-angle scattering off magnetic turbulence, gaining energy efficiently. Recent studies indicate that protons may be accelerated up to \(10^{18}\,\mathrm{eV}\) and heavy nuclei up to \(10^{20}\,\mathrm{eV}\) at reverse shocks of the UFOs [20], [32].

The maximum energy achievable via DSA depends primarily on magnetic-field strength, turbulence amplitude, and shock velocity (Eq. 27 ). Owing to their mildly relativistic high velocities, UFO shocks are expected to accelerate CRs efficiently. However, observational constraints on the magnetic-field strength remain largely uncertain. Many theoretical studies therefore introduce a phenomenological parameter, \(\epsilon_B\), to characterize the fraction of shock energy in magnetic fields. Moreover, Bohm diffusion (\(\delta B/B_0=1\)) are often assumed without explicitly treating the physics of magnetic field amplification, and saturation. Such simplifications potentially introduce systematic uncertainties in estimates of \(E_{\max}\) and, consequently, predictions of neutrino and gamma-ray emission.

The nonresonant hybrid instability (NRH instability, also known as the Bell instability) offers a self-consistent mechanism for determining key magnetic-field parameters such as \(\varepsilon_B\) and \(\delta B/B_0\) in the shock vicinity [33]. This instability is driven by the CR current escaping upstream of the shock, which induces a return current in the background plasma and amplifies circularly polarized Alfvén waves. Magnetic-field amplification regions have been observed in SNRs [34]–[36], and these are often interpreted as possible evidence of NRH instability.

NRH instability has been investigated both analytically and numerically. Analytically, [33] derived the linear growth theory. Numerically, a wide range of methods has been applied, including full particle-in-cell (PIC) simulations [37], kinetic hybrid approaches [38], [39], MHD-PIC methods [40], and hybrid MHD–Vlasov-Fokker-Planck (VFP) models [41], [42]. While PIC simulations capture kinetic effects from shock formation to nonthermal injection, they are computationally prohibitive for tracing particle acceleration orders of magnitude larger than the injection. Hybrid MHD–VFP models reduce computational costs by expanding the CR Boltzmann equation into multipoles, but these studies in most case fixed the CR current, preventing a fully self-consistent treatment of particle acceleration and field amplification.

[43]–[45] developed a numerical code capable of simultaneously evolving the background plasma and CR current in one-dimensional setups. Their results demonstrated that finite spatial extent of the upstream CR current supresses the growth of NRH instability compared to previous theoretical models that assume infinite spatial extent of the CR current. Nevertheless, under favorable conditions, NRH instability can amplify magnetic fields sufficiently to enable \(10^{15}~\mathrm{eV}\) acceleration in very early stage of young SNRs.

In this work, we employ the numerical code of [43] to investigate the growth and saturation of NRH instability in the context of UFOs. By systematically varying control parameters such as background magnetic field \(B_0\), injection rate \(\eta\), and initial amplitude of the magnetic fluctuation \(\xi_{B,\mathrm{ini}}\) defined by Eq. 21 , we assess whether the maximum CR energy \(E_{\max}\) in UFO shocks can be determined self-consistently, without resorting to phenomenological parameters.

The structure of this paper is as follows. In Sec. 2, we describe the governing MHD–CR Boltzmann equations, diffusion coefficient, and the adopted boundary and initial conditions, including \(B_0\), \(\xi_{B,\mathrm{ini}}\), \(\eta\), and \(\mathrm{p}\gamma\) cooling. In Sec. 3, we present the main results: (i) simulations without the NRH term reproduce analytical DSA, (ii) in weak-field regimes (\(B_0=10^{-5}\)–\(10^{-4}\,\mathrm{G}\)), NRH instability amplifies the magnetic field, driving \(E_{\max}\) and \(\varepsilon_B\) to converge regardless of \(\xi_{B,\mathrm{ini}}\), (iii) in strong-field regimes (\(B_0\gtrsim10^{-3}\,\mathrm{G}\)), NRH instability is ineffective and \(E_{\max}\) depends on the initial conditions, (iv) variations in interstellar medium (ISM) density affect acceleration efficiency, and (v) \(\mathrm{p}\gamma\) cooling can limit the maximum energy.

2 Basic Equations and Simulation Setup↩︎

To investigate particle acceleration in UFOs, we adopt the following physical assumptions. More detailed initial conditions are presented in Sec. 2.2.

1. When a UFO collides with the ISM, it generates a shock structure composed of a reverse shock, contact discontinuity, and forward shock. In this study, we focus on the local region around the reverse shock and analyze particle acceleration there (see Fig. 2).

2. CRs are assumed to be pure protons.

3. The wind is assumed to be steady and uniform¹.

4. The background magnetic field is assumed to be uniform in the \(x\) direction, \[B_x \equiv B_0 = \mathrm{const}. \label{eq:background95magnetic95field95UFO}\tag{1}\] In addition, broadband circular-polarized Alfvén-wave turbulence is superposed on the \(y\) and \(z\) components of magnetic field (see, Sec. 2.2.5 for detail). Note that the NRH instability can grow if there is \(B_x\) with coherent length larger than the scale of NRH instability given below by Eq. 22 .

5. The spatial dependence of physical quantities is restricted to one dimension along the \(x\)-axis, and curvature effects are neglected. On the propagation timescale of UFOs across \(\sim\mathrm{pc}\) scales, curvature effects on both fluid dynamics and CR transport are negligible.

2.1 Basic Equations↩︎

In this study, we solve the coupled system of magnetohydrodynamic (MHD) equations for the fluid component and telegraph-type diffusion–convection equations for the CR distribution function [41], [43]–[45]. The MHD equations are expressed as follows. The continuity equation is given by \[\frac{\partial \rho}{\partial t} + \frac{\partial (\rho v_x)}{\partial x} = 0. \label{eq:continuity95equation95Bell}\tag{2}\] The momentum equations are \[\frac{\partial (\rho v_x)}{\partial t} + \frac{\partial}{\partial x} \left(\rho v_x^2 + P + \frac{B_y^2 + B_z^2}{8 \pi}\right) = 0, \label{eq:EOM95Bell95x}\tag{3}\] \[\frac{\partial (\rho v_y)}{\partial t} + \frac{\partial}{\partial x} \left(\rho v_x v_y - \frac{B_x B_y}{4 \pi}\right) = -\frac{1}{c} j_x^{(\mathrm{ret})} B_z, \label{eq:EOM95Bell95y}\tag{4}\] \[\frac{\partial (\rho v_z)}{\partial t} + \frac{\partial}{\partial x} \left(\rho v_z v_x - \frac{B_z B_x}{4 \pi}\right) = \frac{1}{c} j_x^{(\mathrm{ret})} B_y. \label{eq:EOM95Bell95z}\tag{5}\] The energy equation is represented by \[\frac{\partial \epsilon}{\partial t} + \frac{\partial}{\partial x} \left\{ v_x \left( \epsilon + P + \frac{B_y^2 + B_z^2}{8 \pi} \right) - B_x \frac{\boldsymbol{B} \cdot \boldsymbol{v}}{4 \pi} \right\} = 0. \label{eq:energy95equation95Bell}\tag{6}\] The induction equations for \(B_y\) and \(B_z\) are \[\frac{\partial B_y}{\partial t} = \frac{\partial}{\partial x} \left( B_x v_y - B_y v_x \right), \label{eq:induction95equation95Bell95y}\tag{7}\] \[\frac{\partial B_z}{\partial t} = \frac{\partial}{\partial x} \left( B_x v_z - B_z v_x \right). \label{eq:induction95equation95Bell95z}\tag{8}\] Here, the total energy density is defined as \[\epsilon \equiv \frac{P}{\gamma-1}+\frac{1}{2} \rho v^2+\frac{B_y^2+B_z^2}{8 \pi}, \label{eq:definition95energy95density}\tag{9}\] which includes internal, kinetic, and magnetic energies, where \(\gamma\) is the adiabatic index. The return current \(j_x^{(\mathrm{ret})}=-j_x^{(\mathrm{CR})}\) maintains charge neutrality against the CR current \(j_x^{(\mathrm{CR})}\). It is expressed in terms of the anisotropic part of the distribution function \(f_1\), as defined later in Eq. 11 .

The CR distribution function is expanded with respect to the pitch-angle cosine relative to the background magnetic field \(\boldsymbol{B}_0\) as follows, \[f(x, \boldsymbol{p}) = f_0(x, p) + \frac{p_x}{p} f_1(x, p), \label{eq:definition95of95CR95distribution95Bell}\tag{10}\] where \(f_0(x,p)\) denotes the isotropic part and \(f_1(x,p)\) the anisotropic part. The CR current can then be written as \[j_x^{(\mathrm{CR})}=-j_x^{(\mathrm{ret})}=c e \int_{p_{\min }}^{p_{\max }} f_1(x, p) \frac{4 \pi}{3} p^2 \mathrm{~d} p. \label{eq:return95current95Bell}\tag{11}\] By defining \[F_0 \equiv f_0 p^3, \quad F_1 \equiv f_1 p^3,\] the evolution equations for \(F_0\) and \(F_1\) can be expressed as \[\begin{align} \frac{\partial F_0}{\partial t} &+ \frac{\partial\left(v_x F_0\right)}{\partial x} - \frac{1}{3} \frac{\partial v_x}{\partial x} \frac{\partial F_0}{\partial \ln p} \\ &= -\frac{c}{3} \frac{\partial F_1}{\partial x} + Q_{\mathrm{inj}} p^3 - t_{\mathrm{p}\gamma}^{-1}(p)F_0, \end{align} \label{eq:CR95evolution95equation1}\tag{12}\] \[\frac{\partial F_1}{\partial t} + \frac{\partial\left(v_x F_1\right)}{\partial x} = -c \frac{\partial F_0}{\partial x} - \frac{c^2}{3 D_{\|}(p, B)} F_1. \label{eq:CR95evolution95equation2}\tag{13}\] Here, \(Q_{\mathrm{inj}}\) is the CR injection rate, \(t_{\mathrm{p}\gamma}(p)\) is the timescale of \(\mathrm{p}\gamma\) interactions near the AGN source (see Fig. 1), and \(D_{\|}(p,B)\) is the diffusion coefficient. In this work, we incorporate the \(\mathrm{p}\gamma\) energy-loss timescale derived by [20] into the numerical framework of [43]. Note that Eqs. 12 and 13 are reduced to the diffusion convection equation in the limit \(c\rightarrow \infty\) (see, [44] for numerical tests).

The injection rate \(Q_{\mathrm{inj}}\) is modeled as follows. A fraction \(\eta\) of upstream fluid particles with momentum \(p_{\text{inj}}\) are assumed to be injected into the acceleration process [51]: \[Q_{\mathrm{inj}}(x, p)=\frac{\eta n_{\mathrm{wind}} v_{\mathrm{sh}} }{4 \pi p_{\text{inj} }^2} \delta\left(p-p_{\text{inj} }\right) \delta\left(x-x_{\mathrm{sh}}\right),\] where \(n_{\mathrm{wind}}\) is the upstream wind density, \(v_{\mathrm{sh}}\) is the shock velocity in the upstream rest frame, and \(x_{\mathrm{sh}}\) is the shock position. The relation between \(\eta\) and \(p_{\text{inj}}\) is given by \[\eta \equiv \frac{\int_{p_{\text{inj }}}^{\infty} \exp \left(-\frac{p^2}{2 m_{\mathrm{g}} k_{\mathrm{B}} T_{\text{sh }}}\right) \mathrm{d} p}{\int_0^{\infty} \exp \left(-\frac{p^2}{2 m_{\mathrm{g}} k_{\mathrm{B}} T_{\text{sh }}}\right) \mathrm{d} p},\] with \(m_{\mathrm{g}}\) the mean gas mass and \(T_{\text{sh}}\) the downstream temperature. For strong shocks with \(\gamma=5/3\), the Rankine–Hugoniot relation gives \(T_{\text{sh}}=3m_{\text{g}}v_{\text{sh}}/(16k_{\text{B}})\). We adopt \(m_{\text{g}}=1.27m_{\text{p}}\), consistent with the solar composition inferred from Cygnus A emission-line ratios [52]. For numerical reasons, the momentum range of CR is restricted to \(p_{\min }>p_{\mathrm{inj}}\) up to \(p_{\max}\). Between \(p_{\mathrm{inj}}\) and \(p_{\min}\), we assume the standard DSA spectrum \(f_0(x_{\text{sh}})\propto p^{-4}\) [43]–[45]. Accordingly, the injection rate used in simulations is rewritten as \[Q_{\mathrm{inj}}(x, p)=\frac{\eta n_{\mathrm{wind}} v_{\mathrm{sh}} p_{\mathrm{inj}}}{4 \pi p_{\min }^3}\frac{1}{\Delta p_{\text{min}}}\frac{1}{\Delta x}, \label{eq:injection95for95simulation}\tag{14}\] where the delta functions are replaced by grid intervals.

The diffusion coefficient is expressed as \[D_{\|}(p, B) = \frac{4}{3 \pi} \frac{\max \left(B_0^2, \delta B^2\right)}{\delta B^2} \frac{c E_{\mathrm{CR}}}{e \max \left(B_0, \delta B\right)}, \label{eq:diffusion95coefficient95pitch95angle}\tag{15}\] where \(\delta B^2 = B_y^2 + B_z^2\) and \(E_{\mathrm{CR}}\) is the CR energy. For \(\delta B < B_0\), Eq. 15 reduces to the pitch-angle scattering coefficient, while for \(\delta B \geq B_0\) it corresponds to the Bohm limit.

The validity of Eq. 15 is supported by two arguments. First, kinetic hybrid simulations [53] and test-particle calculations [54] have shown that NRH-amplified fields drive diffusion close to the Bohm limit. Second, NRH instability alone amplifies only small-scale magnetic fluctuations, with wavelengths much shorter than the gyroradius at \(E_{\max}\). Filamentation instability can transfer this energy to larger scales, amplifying turbulence up to the gyroradius scale [55], [56].

2.2 Initial Conditions↩︎

2.2.1 Global Fluid Initial Conditions↩︎

The wind launched from a supermassive black hole, including UFOs, is expected to form a bubble structure similar to stellar winds, as suggested by both theory and observations (see Fig. 2 left) [18], [57]–[61]. We performed one-dimensional MHD simulations of a Riemann problem, where a stationary homogeneous ISM is placed on the left side of \(x=0\) and a leftward-propagating wind is injected on the right side. The time evolution of this system generates a shock structure consisting of a reverse shock, contact discontinuity, and forward shock. The resulting density distribution of the bubble is shown in Fig. 2 (left). The detailed quantative conditions for the left and right states of the Riemann problem are given in Sec. 2.2.3 and 2.2.4. Solving Eqs. 12 and 13 for the whole system is computationally expensive. To reduce it, we extract only the local region around the reverse shock and solve Eqs. 12 and 13 in couple with the Bell-MHD Eqs. 2 9 (see, right panel of Fig. 2).

2.2.2 Region Around the Reverse Shock↩︎

The simulation domain is set to \(L_{\text{box}} = 10~\text{pc}\) with a runtime of \(t_{\text{final}} = 100~\text{yr}\).² Typical wind velocities are \(v_{\text{wind}} \sim 0.1c\)–\(0.4c\) [26], [27], [47], [50], [62]–[70]. This setup corresponds to the propagation timescale of a UFO traveling about 1 pc from the central black hole. Observationally, UFOs are typically detected at distances \(<10^{18}~\text{cm}\) [71]–[73], placing the 1 pc scale within the observed distribution range.

The wind and ISM parameters are summarized in Tab. 1. The mass outflow rate is set to \(\dot{M} \sim 0.01\)–\(1~M_{\odot}/\mathrm{yr}\) [27], [70], [71], [74].

2.2.3 Wind Region↩︎

The wind density is expressed as follows: \[\begin{align} \rho_{\text{wind}}(r) &= \frac{\dot{M}}{4\pi r^2 v_{\text{wind}}} \\ &= 8.8\times10^{-24} ~\text{g}~\text{cm}^{-3}\\ &\times\left(\frac{r}{1~\text{pc}}\right)^{-2} \left(\frac{v_{\text{wind}}}{0.2c}\right)^{-1} \left(\frac{\dot{M}}{0.1~M_{\odot}/\text{yr}}\right). \end{align} \label{eq:rho95wind95initial}\tag{16}\] The sound speed is defined as \[c_{\text{s}}\equiv\sqrt{\gamma\frac{P}{\rho}}. \label{eq:acoustic95wave95UFO}\tag{17}\] Using this, the Mach number is given by \[\mathcal{M}_{\text{wind}}\equiv\frac{v_{\text{wind}}}{c_{\text{s, wind}}}.\] The wind pressure can then be expressed as \[\begin{align} P_{\text{wind}} &= \rho_{\text{wind}}\frac{1}{\gamma} \left(\frac{v_{\text{wind}}}{\mathcal{M}_{\text{wind}}}\right)^2 \\ &= 4.7\times10^{-7}~\text{dyn/cm}^{2} ~\left(\frac{r}{1~\text{pc}}\right)^{-2} \left(\frac{\mathcal{M}_{\mathrm{wind}}}{20}\right)^{-2}\\ &\quad \times \left(\frac{\gamma}{5/3}\right)^{-1} \left(\frac{v_{\text{wind}}}{0.2c}\right) \left(\frac{\dot{M}}{0.1~M_{\odot}/\text{yr}}\right). \end{align} \label{eq:P95wind95initial}\tag{18}\]

2.2.4 ISM Region↩︎

For the ISM region, we assume a stationary homogeneous medium, representative of the narrow line region where relatively low-density gas is observed³. The density and temperature of narrow line region, estimated from emission-line diagnostics, are typically \(n_{\text{ISM}} \sim 10^2\)–\(10^4/\text{cm}^3\) and \(T_{\text{ISM}} \sim 1.0\)–\(2.5 \times 10^4~\text{K}\) [77] (see Sec. 6 of [78]). The ISM density is given by \[\begin{align} &\rho_{\text{ISM}} = m_{\text{g}}n_{\text{ISM}} \\ &= 2.1\times10^{-22}~\text{g/cm}^3 \left(\frac{m_{\text{g}}/m_{\text{p}}}{1.27}\right) \left(\frac{n_{\text{ISM}}}{10^2/\text{cm}^3}\right). \end{align} \label{eq:rho95ISM95initial}\tag{19}\] We consider three cases with \(n_{\text{ISM}} = 1, 10^2, 10^4~\text{cm}^{-3}\) (Tab. 1). The ISM pressure is represented as \[\begin{align} &P_{\text{ISM}} = n_{\text{ISM}}k_{\text{B}}T_{\text{ISM}} \\ &= 2.2\times10^{-10}~\text{dyn/cm}^{2} \left(\frac{T_{\text{ISM}}}{1.6\times10^4~\text{K}}\right) \left(\frac{n_{\text{ISM}}}{10^2/\text{cm}^3}\right). \end{align} \label{eq:P95ISM95initial}\tag{20}\]

2.2.5 Initial Magnetic Fluctuations↩︎

We superpose broadband circular-polarized Alfvén waves on to the background field in Eq. 1 . The fluctuation wavelengths are chosen such that the minimum wavelength is smaller than the maximum growth scale of the NRH instability (Eq. 22 ), while the maximum wavelength is set to be equal to the box size: \[\lambda_{\text{min}}<\lambda_{\text{B}},\quad \lambda_{\text{max}}=L_{\text{box}}=3.0\times10^{19}~\text{cm}.\] The initial fluctuation spectrum follows a Kolmogorov form, \(P_{k}(k)\propto k^{-\frac{5}{3}}\).

The fluctuation strength is defined as \[\xi_{B,\mathrm{ini}} \equiv \left\langle \frac{B_y^2 + B_z^2}{B_0^2} \right\rangle. \label{eq:definition95magnetic95fluctuation95strength}\tag{21}\] We adopt two cases: strong initial turbulence with \(\xi_{B,\mathrm{ini}}=0.1\) and weak turbulence with \(\xi_{B,\mathrm{ini}}=0.01\).

As shown in Tab. 1, we vary \(B_0\), \(\xi_{B,\mathrm{ini}}\), \(n_{\text{ISM}}\), and the presence of \(\mathrm{p}\gamma\) cooling to evaluate their impact on particle acceleration. Since observational constraints on \(B_0\) in UFOs are limited, we explore a wide range of field strengths.

2.2.6 Injection Rate↩︎

The injection rate \(\eta\) appearing in Eq. 14 is set to \(\eta=10^{-4}\) and \(10^{-3}\), as summarized in Tab. 1. Varying \(\eta\) changes the number density of CRs, thereby altering the growth efficiency of the NRH instability. Since UFO injection rates are poorly constrained observationally, we adopt \(\eta=10^{-4}\)–\(10^{-3}\) as a representative range, guided by studies of DSA in SNRs [79], [80]. This corresponds to cases where a few to tens of percent of the upstream kinetic energy is transferred into CR energy.

2.3 Boundary Conditions↩︎

For the MHD component, both the left and right boundaries are fixed. At the left boundary, we apply the wind parameters given by Eqs. 16 and 18 , with a fixed wind velocity of \(0.2c\). At the right boundary, we impose the stationary ISM parameters described by Eqs. 19 and 20 .

For the CR component, in physical space, the isotropic part \(F_0\) and anisotropic part \(F_1\) of the CR distribution function are set to zero at both boundaries⁴. In momentum space, the computational domain is defined as \[p_{\text{min}}=10^{12}~\text{eV}~c^{-1},\quad p_{\text{max}}=10^{19}~\text{eV}~c^{-1},\] with the exception that for the \(B_0=10^{-2}~\mathrm{G}\) model the minimum momentum is set to \(p_{\text{min}}=10^{13}~\text{eV}~c^{-1}\).⁵ For \(p_{\text{inj}}<p<p_{\text{min}}\), the isotropic component \(F_0\) follows the DSA spectrum \(f_0(x_{\text{sh}})\propto p^{-4}\) as expressed in Eq. 14 . Outside the numerical domain of the momentum, the anisotropic component \(F_1\) is set to zero.

2.4 Numerical Resolution↩︎

2.4.1 Spatial Resolution↩︎

The spatial resolution must satisfy three criteria simultaneously:

1. The maximum growth wavelength \(\lambda_{\text{NRH}}\) of the NRH instability must be resolved with at least 32 cells (see Fig. 21 of [44]).

2. The diffusion length \(l_{\text{diff}}\) of the lowest-energy CRs (\(p_{\text{min}}\)) must be resolved with at least 10 cells (see Fig. 20 of [44]).

3. The shortest wavelength of the initial magnetic fluctuations, \(\lambda_{\text{min}}\), must be sufficiently resolved to avoid numerical dissipation.

The maximum growth wavelength \(\lambda_{\text{NRH}}\) is given by [33] as \[\begin{align} \lambda_{\text{NRH}} &= \frac{B_0c}{j_{\text{CR}}} \\ &= 4.9\times10^{14}~\text{cm}\\ &\times\left(\frac{B_0}{10^{-5}~\text{G}}\right) \left(\frac{j_{\text{CR}}}{6.2\times10^{-10}~\text{esu s}^{-1}\text{cm}^{-2}}\right)^{-1}, \end{align} \label{eq:maximum95growth95wavelength95Bell}\tag{22}\] where \(j_{\text{CR}}\) is the CR current at the location where the NRH e-folding number \(t_{\text{adv}}/t_{\text{NRH}}\) reaches its maximum. Test calculations estimate this value as \(6.2\times10^{-10}~\text{esu s}^{-1}\text{cm}^{-2}\) for the fiducial model with \((B_0,~\eta,~\xi_{B,\mathrm{ini}})=(10^{-5}~\mathrm{G},~10^{-4},~0.1)\). The details of the NRH e-folding number are discussed in Sec. 3.1.3.

The diffusion length of the lowest-energy CRs is expressed as follows⁶: \[\begin{align} l_{\text{diff}} &\equiv \frac{D_\|}{v_{\text{sh}}} = \frac{4}{3\pi}\xi_{B}^{-1}\frac{cE_{\text{CR}}}{eB_0v_{\text{sh}}} \\ &= 8.5\times10^{14}~\text{cm} \left(\frac{\xi_{B}}{0.1}\right)^{-1} \left(\frac{E_{\text{CR}}}{10^{12}~\text{eV}}\right) \\ &\quad \times \left(\frac{B_0}{10^{-4}~\text{G}}\right)^{-1} \left(\frac{v_{\text{sh}}}{5.0\times10^9~\text{cm}/\text{s}}\right)^{-1}. \end{align} \label{eq:diffusion95length95UFO95bunkai}\tag{23}\] Here, \(v_{\text{sh}}\) is the shock velocity in the upstream rest frame.

To prevent numerical dissipation of the shortest Alfvén wavelengths, the spatial resolution must be chosen carefully. Previous studies analyzed the numerical dissipation of circularly polarized Alfvén waves using the second-order Roe flux method [81]. Based on their results, the required number of numerical cells per wavelength \(n_x\) for the shortest mode can be estimated as⁷ \[\begin{align} n_x>224 &\left(\frac{\rho_{\text{wind}}}{8.78 \times 10^{-24}~\mathrm{g/cm}^3}\right)^{\frac{1}{4}} \left(\frac{B_0}{10^{-2}~\mathrm{G}}\right)^{\frac{1}{2}} \\ &\times \left(\frac{t_{\text{final}}}{100~\mathrm{yr}}\right)^{\frac{1}{2}} \left(\frac{\lambda}{10^{-3}~\mathrm{pc}}\right)^{-\frac{1}{2}}. \end{align} \label{eq:condition95magnetic95fluctuation95divide}\tag{24}\] Accordingly, the grid resolution is set to satisfy Eq. 24 for the shortest initial wavelength \(\lambda_{\text{min}}\). Combining all conditions, we adopt the following resolution: \[\Delta x=1.7\times10^{13}~\text{cm},\quad N_{\text{cell},x}=1,760,000. \label{eq:kuukan95bunkainou95UFO}\tag{25}\]

2.4.2 Momentum-Space Resolution↩︎

The momentum space is divided uniformly in logarithmic intervals from \(p_{\text{min}}\) to \(p_{\text{max}}\). We use \(N_{\text{cell},p} = 64\) cells, consistent with the convergence tests of [44].

2.4.3 Time Step Requirement↩︎

Since the evolution equations for the CR distribution functions, Eqs. 12 and 13 , are hyperbolic, the time step must satisfy the Courant–Friedrichs–Lewy (CFL) condition, expressed as \[\Delta t\leq C_{\text{CFL}}\frac{\Delta x}{v_{\text{CR}}}\sim 3\times10^3~\text{s},\] where Eq. 25 is used. The CR velocity is assumed to be \(v_{\text{CR}} = c/\sqrt{3}\), corresponding to free streaming at relativistic speeds. We adopt \(C_{\text{CFL}}=0.8\).

3 Simulation Results↩︎

3.1 Variation of Background Magnetic Field Strength↩︎

In this section, we discuss the results for the Fiducial models in Tab. 1, which assumes a mass outflow rate of \(0.1~M_{\odot}/\mathrm{yr}\) and an ISM density of \(10^2~\mathrm{cm}^{-3}\), with \(B_0\) varied from \(10^{-5}\) to \(10^{-2}~\mathrm{G}\). Tab. 2 summarizes the maximum CR energy for each \(B_0\) with \(\xi_{B,\mathrm{ini}}=0.1\), comparing cases with and without NRH instability and for different injection rates. Tab. 3 focuses on \(\eta=10^{-4}\), contrasting cases with initial magnetic fluctuations \(\xi_{B,\mathrm{ini}}=0.1\) and \(\xi_{B,\mathrm{ini}}=0.01\). These results show that for \(B_0=10^{-5}\) and \(10^{-4}~\mathrm{G}\), the maximum energy converges to nearly the same value regardless of \(\xi_{B,\mathrm{ini}}\) owing to NRH instability (see Sec. 3.1.3). In contrast, for \(B_0>10^{-3}~\mathrm{G}\) the NRH instability becomes inefficient, and the maximum energy is determined by \(\xi_{B,\mathrm{ini}}\) (see Sec. 3.1.4).

Tab. 4 further presents the magnetic energy fraction \(\varepsilon_B\), defined as the ratio of magnetic energy to upstream kinetic energy, with \(\varepsilon_B^{\mathrm{up}}\) and \(\varepsilon_B^{\mathrm{down}}\) referring to the upstream and downstream regions, respectively: \[\begin{align} \varepsilon_B^{\mathrm{up}} &\equiv \frac{\left\langle \dfrac{\delta B^2}{8\pi} \right\rangle_{[0,\,l_d]}}{\tfrac{1}{2}\rho_{\mathrm{wind}}\,v_{\mathrm{wind}}^2}, \\[1ex] \varepsilon_B^{\mathrm{down}} &\equiv \frac{\left\langle \dfrac{\delta B^2}{8\pi} \right\rangle_{[-l_d,\,0]}}{\tfrac{1}{2}\rho_{\mathrm{wind}}\,v_{\mathrm{wind}}^2}. \end{align}\] Here, \(l_{\mathrm{d}}\equiv\frac{4}{3\pi}\frac{cE_\mathrm{max}(\mathrm{No~NRH})}{eB_0v_{\mathrm{sh}}}\) represents the diffusion length estimated from the maximum energy in the absence of NRH instability (Tab. 2), assuming \(\xi_{B,\mathrm{ini}}=1\).⁸ In all cases, \(\varepsilon_B\) is smaller upstream than downstream, indicating that magnetic amplification in the upstream region controls the overall efficiency of particle acceleration as expected.

3.1.1 Analytical prediction from diffusive shock acceleration↩︎

In the framework of the DSA model, the energy spectrum of CRs can be analytically derived. In the limit of infinite Mach number, the isotropic distribution function is expressed as [2], [3], [29], [30] \[f_0(p) \propto p^{-4}.\]

The maximum acceleration energy can also be estimated within the DSA framework. The acceleration timescale can be estimated by using Eq. 21 , \[t_{\mathrm{acc}} \equiv \frac{D_{\|}}{v_{\mathrm{sh}}^2} = \frac{4}{3 \pi} \frac{c E_{\mathrm{CR}}}{e B_0 v_{\mathrm{sh}}^2}\xi_{B}^{-1}. \label{eq:acceleration95time95UFO}\tag{26}\] By equating \(t_{\mathrm{acc}}\) with the time \(t\), the maximum energy is obtained as \[\begin{align} E_{\text{max}} &= \frac{3\pi}{4}\frac{e B_0 v_{\mathrm{sh}}^2}{c}t\xi_{B} \\ &\sim 1.8 \times 10^{15}~\mathrm{eV} \left(\frac{B_0}{10^{-5}~ \mathrm{G}}\right)\\ &\times\left(\frac{v_{\mathrm{sh}}}{5.0\times10^9~\text{cm/s}}\right)^2 \left(\frac{\xi_{B}}{0.1}\right) \left(\frac{t}{100~ \mathrm{yr}}\right), \end{align} \label{eq:order95estimation95maximum95energy95CR95Bell}\tag{27}\] which will be used in the following discussions.

3.1.2 Case without NRH instability↩︎

We first examine the simulation results without the NRH term and compare them with the analytical predictions of DSA. Fig. 3 shows the CR energy spectra for \(B_0=10^{-5}~\mathrm{G}\), with and without NRH instability. In both cases, the isotropic distribution function follows the analytical DSA prediction \(p^{-4}\) with good agreement.

The maximum CR energies without the NRH term are consistent with Eq. 27 . As shown in Tab. 2, the simulation results match the order-of-magnitude estimate of Eq. 27 within a factor of two. In addition, Tab. 3 shows that reducing the initial fluctuation amplitude \(\xi_{B,\mathrm{ini}}\) by an order of magnitude lowers the maximum energy by a similar factor, consistent with Eq. 27 .

3.1.3 Case of weak background magnetic field with NRH instability (\(B_0=10^{-5},\,10^{-4}~\mathrm{G}\))↩︎

The key result for \(B_0=10^{-5},\,10^{-4}~\mathrm{G}\) is that, once NRH instability is included, \(E_{\max}\) converges to nearly a same value for a given \(B_0\), even when the initial fluctuation amplitude \(\xi_{B,\mathrm{ini}}\) is varied (Tab. 3). In particular, when \(\eta=10^{-4}\), reducing \(\xi_{B,\mathrm{ini}}\) from \(0.1\) to \(0.01\) changes \(E_{\max}\) by only factor 2 or less. This behavior indicates that even when the initial field amplitude is small, the growth and saturation of NRH instability render the final value of \(E_{\max}\) insensitive to \(\xi_{B,\mathrm{ini}}\).

The same convergence trend is evident in the magnetic energy fraction \(\varepsilon_B\), as summarized in Tab. 4. When NRH instability is included, \(\varepsilon_B\) converges to similar values for \(\xi_{B,\mathrm{ini}}=0.1\) and \(\xi_{B,\mathrm{ini}}=0.01\). Moreover, in almost all cases, \(\varepsilon_B^{\mathrm{up}}\) is smaller than \(\varepsilon_B^{\mathrm{down}}\), implying that magnetic-field amplification upstream regulates the efficiency of particle acceleration.

Fig. 4 shows the spatial distribution of \(\delta B/B_0\) upstream of the shock at \(t_{\text{final}}=100~\mathrm{yr}\) for \((B_0,~\eta,~\xi_{B,\mathrm{ini}})=(10^{-5}~\mathrm{G},~10^{-4},~0.1)\). Without NRH instability (left), the fluctuations remain close to the initial value, while with NRH instability (right), the magnetic field is amplified by nearly an order of magnitude near the shock front. The amplification decreases with distance from the shock, as clarified later in the discussion of Fig. 5. A similar trend is obtained for \(B_0=10^{-4}~\mathrm{G}\). Meanwhile, in the downstream region, magnetic turbulence is enhanced by shock compression in both cases, with and without the NRH term.

Fig. 5 presents the upstream profiles of the CR current density \(j_{\text{CR}}\) (left), defined by Eq. 11 , and the e-folding number of NRH instability \(t_{\text{adv}}/t_{\text{NRH}}\) (right). The e-folding number peaks in a finite region ahead of the shock, and Fig. 4 confirms that magnetic-field amplification occurs predominantly inside this region. The inverse of the linear growth rate of the NRH instability (growth timescale) is represented by [33] \[\begin{align} t_{\text{NRH }}&=\sqrt{\frac{\rho}{\pi}} \frac{c}{j_{\text{CR}}}\\ &\sim 2.6 ~\mathrm{yr}~\left(\frac{\rho_{\text{wind }}}{8.78 \times 10^{-24} \mathrm{~g} / \mathrm{cm}^3}\right)^{\frac{1}{2}}\\ &\times\left(\frac{j_{\text{CR}}}{6.2\times10^{-10}~ \mathrm{esu~s}^{-1} \mathrm{~cm}^{-2}}\right)^{-1}, \end{align} \label{eq:NRH95linear95timescale}\tag{28}\] while the advection time before the shock overtakes that point is expressed as \[\begin{align} t_{\text{adv}}&\equiv \frac{x-x_{\text{sh}}}{v_{\text{sh}}}\\ &\sim 13~{\text{yr}}~\left(\frac{v_{\mathrm{sh}}}{5.0\times10^9~\text{cm/s}}\right)^{-1} \left(\frac{x-x_{\text{sh}}}{6.6\times10^{-1}~\text{pc}}\right). \end{align} \label{eq:definition95advection95time95Bell}\tag{29}\] Magnetic-field amplification proceeds efficiently in regions where \(t_{\text{NRH}} \ll t_{\text{adv}}\). In the case of \(B_0=10^{-5}~\mathrm{G}\), the e-folding number reaches a maximum of \(\sim5.0\) at \(x-x_{\text{sh}}\sim6.6\times10^{-1}~\text{pc}\), and significant amplification occurs inside this location⁹. A comparison of the maximum e-folding numbers for weak and strong background magnetic field cases is provided in Tab. 5. The case of strong background fields (\(B_0=10^{-2}~\mathrm{G}\)) is discussed in more depth in Sec. 3.1.4.

Fig. 6 presents the wavelength spectra of magnetic fluctuations extracted from the upstream region (\(0.1\)–\(1~\mathrm{pc}\) ahead of the shock) in the fiducial model with \((\eta,~\xi_{B,\mathrm{ini}})=(10^{-4},~0.1)\), including the NRH instability. The left panel corresponds to \(B_0=10^{-5}~\mathrm{G}\), and the right panel to \(B_0=10^{-4}~\mathrm{G}\). In this analysis, the Fourier components \(\tilde{B}_i(k)\), the power spectrum \(P_k(k)\), and the wavelength spectrum \(P_\lambda(\lambda)\) are defined as follows, \[\require{physics} \begin{align} \tilde{B}_i(k) & = \int B_i~\mathrm{e}^{-2 \pi \mathrm{i} k x} \dd x \quad (i=y, z), \\ P_k(k) & \equiv \left|\tilde{B}_y\right|^2 + \left|\tilde{B}_z\right|^2, \\ P_\lambda(\lambda) & \equiv P_k(k)\left|\frac{\dd k}{\dd \lambda}\right|. \end{align} \label{eq:magnetic95power95spectrum95UFO}\tag{30}\] For a Kolmogorov spectrum (\(\propto k^{-3/5}\)), the wavelength spectrum scales as \(P_\lambda(\lambda)\propto \lambda^{-1/3}\).

In the \(B_0=10^{-5}~\mathrm{G}\) case, the spectrum at the final time peaks at \(\lambda \sim 2.8\times10^{-4}~\mathrm{pc}\), consistent within a factor of two with the analytically expected NRH maximum growth scale \(\lambda_{\text{NRH}}\sim1.6\times10^{-4}~\mathrm{pc}\) (Eq. 22 ). The final spectrum shows rapid amplification from the shortest scales up to \(\lambda_{\text{NRH}}\), followed by a gradual decline toward longer wavelengths. This behavior is explained by the NRH linear growth rate, expressed as [33] \[t_{\text{NRH}}^{-1}(k) = \sqrt{k\left(\frac{B_0 j_{\text{CR}}}{c \rho} - k v_{\text{A}}^2\right)},\] where the Alfvén speed is defined by \[v_{\text{A}} \equiv \frac{B_0}{\sqrt{4\pi \rho}}. \label{eq:definition95Alfven95speed95UFO}\tag{31}\] Accordingly, the minimum unstable wavelength is represented by \[\lambda_{\text{NRH}}^{\min} = \frac{B_0 c}{2 j_{\text{CR}}} = \frac{1}{2}\lambda_{\text{NRH}}.\] Thus, instability grows rapidly from \(\lambda_{\text{NRH}}^{\min}\), reaches maximum growth at \(\lambda_{\text{NRH}}\), and decreases toward longer wavelengths gradually.

For \(B_0=10^{-4}~\mathrm{G}\) (right panel), the NRH instability remains effective, producing amplification near the maximum growth scale. However, the spectral peak shifts to longer wavelengths, reaching \(\sim2.9\times10^{-3}~\mathrm{pc}\). This shift is consistent with the scaling \(\lambda_{\text{NRH}}\propto B_0/j_{\text{CR}}\), derived from Eq. 22 and the growth rate above. Since the spatial structure of \(j_{\text{CR}}\) does not change significantly with \(B_0\), increasing the background field by a factor of ten results in a nearly proportional increase of the peak wavelength.

3.1.4 Cases with Strong Background Magnetic Fields (\(B_0=10^{-3},~10^{-2}~\mathrm{G}\))↩︎

As discussed in the previous section, when the background magnetic field is sufficiently weak, the NRH instability supplies magnetic fluctuations in a self-regulating manner. As a result, for each \(B_0\), \(E_{\max}\) converges to nearly the same value, almost independent of the initial fluctuation amplitude \(\xi_{B,\mathrm{ini}}\). In contrast, for \(B_0 \gtrsim 10^{-3}~\mathrm{G}\) this self-regulation breaks down, and \(E_{\max}\) transitions to a regime determined by the initial conditions (\(B_0,~\xi_{B,\mathrm{ini}}\)). For example, \(B_0=10^{-3}~\mathrm{G}\) corresponds to a boundary case, where the enhancement of \(E_{\max}\) due to NRH instability is limited to a factor of 1.7 even for \(\xi_{B,\mathrm{ini}}=0.01\) (Tab. 3). At \(B_0=10^{-2}~\mathrm{G}\) the NRH instability ceases to operate entirely, and the final value of \(E_{\max}\) is governed solely by both the background field strength \(B_0\) and initial turbulence amplitude \(\xi_{B,\mathrm{ini}}\).

The reduced efficiency of NRH instability arises from the significant decrease in the upstream CR current density \(j_{\text{CR}}\). A smaller \(j_{\text{CR}}\) increases the NRH growth timescale \(t_{\text{NRH}}\) in Eq. 28 , limiting the e-foldings achievable within the advection time \(t_{\text{adv}}\). Fig. 7 shows that, for \((B_0,~\eta,~\xi_{B,\mathrm{ini}})=(10^{-2}~\mathrm{G},~10^{-4},~0.1)\), \(j_{\text{CR}}\) at \(t=100~\mathrm{yr}\) is more than two orders of magnitude smaller than in the \(B_0=10^{-5}~\mathrm{G}\) case (Fig. 5), resulting in insufficient amplification as also summarized in Tab. 5. The decline in \(j_{\text{CR}}\) with increasing \(B_0\) is explained by particle transport: a stronger magnetic field reduces the gyroradius for a CR of given energy, making it harder to leak into the upstream region. In the strong background magnetic field regime, NRH instability is therefore suppressed, and \(E_{\max}\) is determined by the initial values of (\(B_0,~\xi_{B,\mathrm{ini}}\)).

Even for strong background magnetic fields (\(B_0=10^{-2}~\mathrm{G}\)), the acceleration efficiency agrees with the analytical DSA prediction (Eq. 27 ) within a factor of 2. On the fluid side, however, clear damping of short-wavelength magnetic fluctuations appears, depending on the initial amplitude. As shown in Fig. 8 (left), the case with \(\xi_{B,\mathrm{ini}}=0.1\) exhibits strong damping for \(\lambda \lesssim 10^{-2}~\mathrm{pc}\), whereas the case with \(\xi_{B,\mathrm{ini}}=0.01\) (right) does not. This behavior is consistent with parametric instabilities, particularly stimulated Brillouin scattering, where a parent Alfvén wave decays into a backward Alfvén wave and a sound wave. The growth timescale of stimulated Brillouin scattering is expressed as follows [82]–[85], \[\begin{align} t_\text{B} &= 2 \xi_{B}^{-\frac{1}{2}}(1+\theta) \sqrt{\frac{\theta}{1-\theta}} \frac{1}{k v_{\text{A}}} \\ &= 0.9~ \mathrm{yr}~ \left(\frac{\xi_{B}}{0.1}\right)^{-\frac{1}{2}} \left(\frac{P_{\text{wind }}}{4.74 \times 10^{-7} ~\mathrm{dyn/cm}^2}\right)^{\frac{1}{4}} \\ &\times\left(\frac{B_0}{10^{-2}~ \mathrm{G}}\right)^{-\frac{3}{2}} \left(\frac{\lambda}{10^{-2}~ \mathrm{pc}}\right) \left(\frac{\rho_{\text{wind }}}{8.78 \times 10^{-24} ~\mathrm{g/cm}^3}\right)^{\frac{1}{2}}, \end{align}\] where \(\theta\) is given by \[\begin{align} \theta &= \frac{c_\text{s}}{v_{\text{A}}}= 0.315 \left(\frac{P_{\text{wind }}}{4.74 \times 10^{-7}~ \mathrm{dyn/cm}^2}\right)^{\frac{1}{2}} \\ &\quad\times\left(\frac{B_0}{10^{-2}~ \mathrm{G}}\right)^{-1} \left(\frac{\gamma}{5/3}\right)^{\frac{1}{2}}. \end{align}\] The scaling \(t_\text{B}\propto \xi_{B}^{-1/2}(k v_{\mathrm{A}})^{-1}\) implies \(t_\text{B}\propto \delta B^{-1}\), so that stronger initial fluctuations (\(\xi_{B,\mathrm{ini}}=0.1\)) lead to faster damping at short wavelengths, while weaker fluctuations (\(\xi_{B,\mathrm{ini}}=0.01\)) suppress damping, consistent with the simulation results (Fig. 8).

In the right panel of Fig. 8, the short-wavelength magnetic energy spectra collapse to \(P_\lambda(\lambda)\simeq \mathrm{const}\), which from Eq. 30 corresponds to \(P_k(k)\propto k^{-2}\), consistent with both theoretical predictions and solar wind observations. Previous MHD simulations have shown that circularly polarized Alfvén waves with broadband spectra undergo time evolution such that, under magnetically dominated conditions (\(c_{\mathrm{s}}\ll v_{\mathrm{A}}\)), strong parametric decay into backward Alfvén waves and sound waves occurs, while the instability is suppressed under gas-pressure-dominated conditions (\(c_{\mathrm{s}}\gg v_{\mathrm{A}}\)) [86]–[88]. Furthermore, [89] derived, based on weak turbulence theory, that the scattered Alfvén-wave spectrum scales as \(\omega^{-2}\). Observations also report that regions dominated by backward-propagating (sunward) scattered waves exhibit magnetic energy spectra with slopes close to \(k^{-2}\) along the background field direction [90]–[92].

The decay of magnetic fluctuations is attributed to physical processes rather than numerical dissipation, supported by following three evidences. First, the shortest initial wavelength \(\lambda_{\text{min}}\sim1.4\times10^{-3}~\text{pc}\) is resolved by 256 grid cells, and Eq. 24 ensures that more than 95% of the Alfvén wave amplitude remains. Second, the measured damping time is consistent with the timescale of stimulated Brillouin scattering. Third, as shown in Fig. 9, density fluctuations appear in the density and pressure, and their amplitudes grow over time. This behavior can be interpreted as the growth of sound waves via stimulated Brillouin scattering.¹⁰

In the present simulations, the decay of magnetic fluctuations has only a limited impact on the efficiency of particle acceleration. This is because the Kolmogorov spectrum used as the initial condition places most of the energy at long wavelengths (\(\sim\)pc scale), which are largely unaffected by parametric instabilities. Even if the short-wavelength components decay, efficient particle acceleration can be maintained as long as the long-wavelength components persist.

However, this result depends on the idealized assumption that large-amplitude (\(\xi_{B,\mathrm{ini}}\sim0.1\)) long-wavelength fluctuations are always present near the UFO shock. In realistic environments, magnetic fluctuations generated near the black hole may undergo significant damping through parametric instabilities before propagating to pc scales, reducing the amplitude to \(\xi_{B}\ll0.1\). Therefore, achieving particle acceleration up to \(\sim10^{18}~\text{eV}\) in UFOs requires that sufficiently strong long-wavelength fluctuations survive to pc scales or that fresh turbulence is generated in situ.

3.2 Variation of ISM Density↩︎

This section examines the impact of ISM density on the efficiency of the NRH instability and the maximum CR acceleration energy (see the ISM low and ISM high models in Tab. 1). As shown in Tab. 6, compared to the case with \(n_{\text{ISM}}=10^2~\mathrm{cm}^{-3}\), the maximum acceleration energy \(E_{\max}\) decreases by approximately one order of magnitude at \(n_{\text{ISM}}=1~\mathrm{cm}^{-3}\), while it increases by a factor of a few at \(n_{\text{ISM}}=10^4~\mathrm{cm}^{-3}\). The NRH instability becomes more efficient at higher ISM densities, leading to larger \(E_{\max}\).

The maximum CR energy depends strongly on the reverse shock velocity (see Eq. 27 ). As the ISM density increases, the shock velocity in the wind rest frame also increases. In the extreme limit \(n_{\mathrm{ISM}}\rightarrow\infty\), corresponding to complete reflection of the reverse shock by the ISM, the upstream velocity in the shock frame becomes \(v_{\text{u}}=v_{\text{wind}}+v_{\text{sh}}\) and the downstream velocity is \(v_{\text{d}}=v_{\text{sh}}\). With \(v_{\text{wind}}=0.2c\) and a compression ratio of 4, one obtains \(v_{\text{u}}=4v_{\text{d}}\). Thus, the shock velocity in the lab frame is \(v^{\prime}_{\text{sh}}=1/3v_{\text{wind}}\), and in the wind rest frame it becomes \(v_{\text{sh}}=4/3v_{\text{wind}}=8\times10^9~\mathrm{cm~s}^{-1}\). Therefore, even as ISM density increases indefinitely, the shock velocity saturates at this value. Tab. 7 shows that as the ISM density increases, the rise in shock velocity becomes smaller, with little change between \(n_{\text{ISM}}=10^2\) and \(10^4~\mathrm{cm}^{-3}\). Consequently, \(E_{\max}\) also shows a saturation trend, consistent with the scaling \(E_{\max}\propto v_{\text{sh}}^2\) in Eq. 27 .

In our simulations, increasing \(n_{\mathrm{ISM}}\) leads to higher \(v_{\text{sh}}\), as shown in Tab. 7, which enhances the efficiency of the NRH instability. A larger \(v_{\text{sh}}\) reduces the CR acceleration time in Eq. 26 , since \(t_{\text{acc}}\propto v_{\text{sh}}^{-2}\) when \(B_0\), \(\xi_{B,\mathrm{ini}}\), and \(E_{\mathrm{CR}}\) are fixed. Consequently, \(t_{\text{acc}}\) becomes shorter, allowing high-energy CRs to be produced more rapidly. These CRs escape further upstream, increasing the CR current density \(j_{\text{CR}}\) ahead of the shock. Figure 10 illustrates the spatial distribution of the upstream CR current at \(t=100~\mathrm{yr}\) for different ISM densities. It clearly shows that a higher ISM density results in a larger CR current escaping into the upstream region. As \(j_{\text{CR}}\) increases, the NRH instability growth time \(t_{\text{NRH}}\) in Eq. 28 decreases, which results in more e-foldings \(t_{\text{adv}}/t_{\text{NRH}}\) and stronger magnetic amplification.

3.3 Effect of p\(\gamma\) Cooling↩︎

In this subsection, we examine the impact of \(\mathrm{p}\gamma\) cooling on the maximum CR acceleration energy with NRH instability. As shown in Tab. 8, \(\mathrm{p}\gamma\) cooling becomes effective when the background magnetic field increases to \(B_0 \sim 10^{-3}~\mathrm{G}\), corresponding to \(E_{\max}\sim10^{17}~\mathrm{eV}\). The suppression is most prominent for \(B_0=10^{-2}~\mathrm{G}\). Fig. 11 illustrates this behavior: without \(\mathrm{p}\gamma\) cooling, CRs reach \(\sim\mathrm{EeV}\) (\(10^{18}~\mathrm{eV}\)), whereas with cooling, the maximum energy is reduced to the sub-EeV regime. In this case, \(E_{\max}\) saturates at a nearly constant value after \(t \sim 60~\mathrm{yr}\).

This behavior can be interpreted from the comparison of timescales in Fig. 1. In standard DSA theory, the maximum energy is determined where the acceleration timescale \(t_{\mathrm{acc}}\) (Eq. 26 ) equals the simulation time \(t_{\mathrm{final}}\). In Fig. 1, this condition corresponds to the intersection of the red dot-dashed line (\(t_{\mathrm{acc}}\)) and the orange dashed line (\(t_{\mathrm{final}}\)). When the \(\mathrm{p}\gamma\) cooling timescale \(t_{\mathrm{p}\gamma}\) becomes shorter than \(t_{\mathrm{final}}\), however, \(E_{\max}\) is instead determined by the intersection of \(t_{\mathrm{acc}}\) with \(t_{\mathrm{p}\gamma}\) (blue solid line). For the \(B_0=10^{-2}~\mathrm{G}\) model, the maximum energy estimated from \(t_{\mathrm{acc}}=t_{\mathrm{final}}\) is \(\sim1.8\times10^{18}~\mathrm{eV}\), while that from \(t_{\mathrm{acc}}=t_{\mathrm{p}\gamma}\) is \(\sim7.1\times10^{17}~\mathrm{eV}\). Thus, \(\mathrm{p}\gamma\) cooling reduces \(E_{\max}\) by a factor of \(\sim0.4\), in good agreement with the numerical results in Fig. 11.

The condition for \(\mathrm{p}\gamma\) cooling to impose a limit on \(E_{\max}\) can be expressed in terms of a critical magnetic field strength \(B_{\mathrm{p}\gamma}\). This corresponds to the point where \(t_{\mathrm{acc}}\), \(t_{\mathrm{final}}\), and \(t_{\mathrm{p}\gamma}\) intersect. Using \(t_{\mathrm{final}}=100~\mathrm{yr}\) and the intersection energy \(E_{\mathrm{CR}}=1.3\times10^{16}~\mathrm{eV}\) between \(t_{\mathrm{final}}\) and \(t_{\mathrm{p}\gamma}\), we estimate \[\begin{align} B_{\mathrm{p}\gamma} &> 7.1\times10^{-5}~ \mathrm{G} \\ & \times \left(\frac{t_{\mathrm{final}}}{100~ \mathrm{yr}}\right)^{-1} \left(\frac{v_{\mathrm{sh}}}{5.0\times10^9~\text{cm/s}}\right)^{-2}\left(\frac{\xi_{B}}{0.1}\right)^{-1}. \end{align}\]

4 Conclusion↩︎

In this work, we evaluated the growth of NRH instability–driven magnetic turbulence and the maximum CR acceleration energy \(E_{\mathrm{max}}\) in reverse shocks of AGN UFOs. Using a numerical framework that self-consistently couples the CR diffusion–convection equation with nonlinear MHD evolution, we examined the dependence on \(B_0\), initial strength of magnetic turbulence \(\xi_{B,\mathrm{ini}}\) in Eq. 21 , and injection rate \(\eta\), while including the effects of NRH instability and \(\mathrm{p}\gamma\) cooling.

Previous PIC simulations at very low CR maximum energies demonstrated strong magnetic field amplification due to the NRH instability that saturates in nonlinear stage [38], [39]. By contrast, [44], [45] showed only moderate amplification in simulations of SNRs consistent with observations. In their simulations, the maximum CR energies were \(1~\mathrm{PeV}\) or less. Our study demonstrates that when the maximum energy is even higher (\(\lesssim1 ~\mathrm{EeV}\)), magnetic field amplification due to NRH instability becomes much weaker. Since the CR current originates from particles escaping upstream near \(E_{\max}\), it is natural that the current diminishes with increasing maximum energy, thereby suppressing NRH instability growth. Many earlier works did not explicitly examine the dependence of magnetic amplification on maximum energy, but our results suggest that realistic theoretical models must describe NRH amplification as a function of \(E_{\max}\).

To accelerate CRs up to \(\sim10^{18}~\mathrm{eV}\) in UFOs, the following conditions must be simultaneously satisfied, which is not as easy as previously thought:

1. Near the reverse shock (within \(\sim1~\mathrm{pc}\)), both the background magnetic field \(B_0\) and turbulent amplitude \(\xi_{B,\mathrm{ini}}\) must be sufficiently large, specifically \(B_0\geq10^{-2}~\mathrm{G}\) and \(\xi_{B,\mathrm{ini}}\geq0.1\). Under these conditions, the NRH instability is ineffective, and the local parameters of the acceleration region determine \(E_{\max}\). In contrast, when the background magnetic field is weak (\(B_0<10^{-4}~\mathrm{G}\)), the NRH instability operates efficiently, but the acceleration efficiency is insufficient to reach the EeV range.

2. The magnetic fluctuation spectrum must be dominated by long wavelengths (Kolmogorov-type). If short wavelengths dominate initially, parametric instabilities such as stimulated Brillouin scattering cause their rapid decay, reducing the acceleration efficiency (Fig. 8).

3. The \(\mathrm{p}\gamma\) cooling must remain inefficient (i.e., the AGN photon field must be weak). [20] showed that EeV-scale acceleration is possible even with \(\mathrm{p}\gamma\) cooling if the magnetic energy fraction is high (\(\epsilon_B \simeq 0.05\)). In our simulations, the downstream magnetic energy fraction eventually reaches \(\epsilon_B^{\mathrm{down}}\simeq0.04\) for the initial condition \((B_0,\,\xi_{B,\mathrm{ini}})=(10^{-2}~\mathrm{G},\,0.1)\) (see Tab. 4), similar to their assumed acceleration conditions. However, the maximum CR energy remains sub-EeV (Fig. 11). Thus, achieving EeV energies likely requires a photon field weaker than that in [20]. The difference is only a factor of a few and is therefore not so serious.

Whether these conditions are realized in UFO environments requires future observational confirmation.

In the weak magnetic field regime (\(B_0\leq10^{-4}~\mathrm{G}\)), the NRH instability self-consistently amplifies magnetic fluctuation \(\xi_{B}\) regardless of \(\xi_{B,\mathrm{ini}}\), and \(E_{\max}\) is automatically determined. For \(\eta=10^{-4}\), \(E_{\max}\) reaches \(\sim10^{16}~\mathrm{eV}\) at \(B_0=10^{-4}~\mathrm{G}\) and \(\sim10^{15}~\mathrm{eV}\) at \(B_0=10^{-5}~\mathrm{G}\). At a higher injection rate \(\eta=10^{-3}\), Tab. 2 shows that \(E_{\max}\) can increase further by a factor of a few.

A transition occurs at \(B_0\gtrsim10^{-3}~\mathrm{G}\), where the linear NRH instability growth time \(t_{\mathrm{NRH}}\) (Eq. 28 ) tends to exceed the upstream advection time \(t_{\mathrm{adv}}\) (Eq. 29 ), so that the NRH instability is ineffective. In this regime, \(E_{\max}\) is controlled by the initial conditions (\(B_0,\;\xi_{B,\mathrm{ini}}\)) and by \(\mathrm{p}\gamma\) cooling (Tabs. 2, 8). The key reason for \(t_{\mathrm{adv}}/t_{\mathrm{NRH}}<1\) is simply due to the suppression of the escaping CR current \(j_{\mathrm{CR}}\) at larger \(B_0\). A stronger magnetic field decreases the gyroradius, making it harder for CRs of a given energy to escape far upstream.

At \(B_0\sim10^{-2}~\mathrm{G}\), acceleration up to sub–\(\mathrm{EeV}\) is possible if the initial magnetic turbulence \(\xi_{B,\mathrm{ini}}\) is long-wavelength dominated. However, if, in realistic UFO environments, most of the magnetic fluctuation energy may resides at scales \(\lambda\lesssim10^{-2}~\mathrm{pc}\), the acceleration efficiency would be strongly reduced. Such short-wavelength fluctuations decay through parametric instabilities, including stimulated Brillouin scattering, and \(E_{\max}\) would then fall below \(10^{17}~\mathrm{eV}\).

Recent XRISM results indicate mass outflow rates of \(\sim100~M_{\odot}\,\mathrm{yr}^{-1}\) [50], and both theory and observations suggest that UFOs are clumpy [49], [50]. Future work must assess how the large mass outflow rates and clumpy structure of UFOs affect efficiency of NRH instability and CR acceleration. Higher mass outflow rates may increase the number of particles contributing to the CR current, thereby enhancing NRH instability. If UFOs are generally clumpy and inhomogeneous, shock geometry would become highly non-uniform, and the acceleration efficiency could differ substantially from the uniform-wind case (e.g., [93] for SNR case). Addressing these issues will require combined theoretical, numerical, and observational efforts.

References↩︎

---

1. Variability of UFOs has been observed on timescales from months down to days [46]–[48]. Moreover, both theoretical and observational studies suggest that UFOs are clumpy and inhomogeneous [49], [50]. These factors may affect particle acceleration and will be considered in future work.↩︎

2. For the fiducial model with \(B_0=10^{-5}~\mathrm{G}\), the box size is increased to \(L_{\text{box}} = 15~\text{pc}\) to prevent the anisotropic CR current (Eq. 11 ) from reaching the boundary.↩︎

3. AGN outflows including UFOs are suggested to form shocks in the narrow line region [75], [76].↩︎

4. The difference between imposing a zero boundary or a free boundary is negligible. Since the simulation domain is sufficiently large, the reverse shock remains far from the right boundary, and CRs never reach it within the simulation time.↩︎

5. At \(p=10^{12}~\text{eV}~c^{-1}\), the diffusion length cannot be resolved with more than ten grid cells given the spatial resolution of Eq. 25 . See Eq. 23 for details.↩︎

6. The shock velocity in the UFO frame is estimated from the time evolution of the position where the fluid velocity discontinuously changes from \(0~\mathrm{cm}/\mathrm{s}\) to beyond \(2\times10^{9}~\mathrm{cm}/\mathrm{s}\).↩︎

7. After \(v_{\text{A}}t_{\text{final}}/\lambda\) periods, the condition for the amplitude of an Alfvén wave with wavelength \(\lambda\) to remain above \(95\%\) can be expressed as \[\begin{align} \left\{1 - 0.20\left(\frac{8}{n_x}\right)^{2}\right\}^{\frac{v_{\text{A}}t_{\text{final}}}{5\lambda}} > 0.95, \end{align} \nonumber\] where the Alfvén speed \(v_{\text{A}}\propto B_0\) is defined by Eq. 31 . In this study, we choose \(n_x\) according to the strength of the background magnetic field. Specifically, we set \(n_x = 16\) for \(B_{0} = 10^{-5}~\mathrm{G}\); \(n_x = 64\) for \(B_{0} = 10^{-4}~\mathrm{G}\) and \(B_{0} = 10^{-3}~\mathrm{G}\); and \(n_x = 256\) for \(B_{0} = 10^{-2}~\mathrm{G}\).↩︎

8. This choice corresponds to the characteristic scale where downstream magnetic fluctuations remain nearly constant due to shock compression, and provides a consistent averaging length across different \(B_0\).↩︎

9. The e-folding number in Fig. 5 (right) represents a local measure of how many times the instability can grow before the shock arrival at each position. The actual number is larger inside the peak because growth initiated farther upstream continues to accumulate over time.↩︎

10. A more rigorous identification of the decay mechanism requires Fourier analysis in both space and time, in order to verify whether peaks associated with stimulated Brillouin scattering appear on the dispersion relation of the sound wave.↩︎