A Simulation Study of Ultra-relativistic Jets. III. Particle Acceleration in FR-II Jets

Since we are interested in FR-II radio galaxies ejected into the ICM of galaxy clusters, we employ the so-called King profile to emulate the density stratification in the typical cluster core region:

$$\begin{eqnarray}&&{\rho }_{\mathrm{ICM}}(r)={\rho }_{c}{\left[1+{\left(\displaystyle \frac{r}{{r}_{c}}\right)}^{2}\right]}^{-3{\beta }_{K}/2},\end{eqnarray} \tag{ 1 }$$

where r is the distance from the cluster center and r_c , ρ_c , and β_K are the core radius, the core density, and a model parameter, respectively. We adopt β_K = 0.5 and r_c = 50 kpc. The hydrogen number density at the cluster center is set to be n_H,ICM = 10⁻³ cm⁻³ as in Paper II, so ρ_c = 2.34 × 10⁻²⁷ g cm⁻³. Furthermore, the cluster core is assumed to be isothermal with T_ICM = 5 × 10⁷ K. In such a numerical setup, we concentrate on the relatively late evolution of the jet as it propagates far away from the much denser environment near the central AGN. To balance the pressure gradient due to the stratified background, an external gravity is imposed, and hydrostatic equilibrium is achieved without a jet.

The length and timescales of jet simulations are normalized with r₀ = r_j and t₀ = r_j /c, respectively. Below, the simulation results are expressed in units of the following normalization: ⁴ ρ₀ = ρ_c , u₀ = c, and P₀ = ρ₀ c² = 2.1 × 10⁻⁶ erg cm⁻³. Then, the pressure at the cluster core corresponds to P_c /P₀ = 7.64 × 10⁻⁶ in dimensionless units.

The simulation domain is represented by a rectangular box in the 3D Cartesian coordinate system. The cluster center is located at the middle of the bottom surface, defined as (x, y, z) = (0, 0, 0). A circular jet nozzle with a radius of r_j = 1 kpc, through which the jet material is injected along the z-direction, is positioned at the cluster center. The outflow boundary condition is imposed on the bottom surface except for the jet nozzle. The continuous boundary condition is applied to the other five sides of the simulation domain.

The jet model can be specified by the jet power, Q_j , the density ratio, η ≡ ρ_j /ρ_c , and the pressure ratio, P_j /P_c , in our simulations. We consider only models with η = 10⁻⁵ and P_j = P_c , and examine the acceleration of UHECRs in representative jets with different Q_j . The jet power is the rate of the energy inflow through the jet nozzle, excluding the mass energy:

$$\begin{eqnarray}&&{Q}_{j}=\pi {r}_{j}^{2}{u}_{j}\left({{\rm{\Gamma }}}_{j}^{2}{\rho }_{j}{h}_{j}-{{\rm{\Gamma }}}_{j}{\rho }_{j}{c}^{2}\right),\end{eqnarray} \tag{ 2 }$$

where h_j = (e_j + P_j )/ρ_j and e_j are the specific enthalpy and the sum of the internal and rest-mass energy densities of the jet inflow, respectively. With fixed r_j , ρ_j , and P_j , Q_j is determined by u_j or Γ_j .

Table 1 lists the models, showing the ranges of jet-flow quantities: Q_j ≈ 3.34 × (10⁴⁵ − 10⁴⁷) ergs⁻¹, u_j /c ≈ 0.9905–0.9999, and Γ_j ≈ 7.3–71. They intend to include the characteristic values inferred for observed FR-II radio jets (Ghisellini & Celotti 2001; Godfrey & Shabala 2013). The first column shows the model name, following the nomenclature of Paper II. The first two elements are the exponents of Q_j and η. For instance, Q45-η5-S is for the model with Q_j ≈ 3.34 × 10⁴⁵ erg s⁻¹ and η = 10⁻⁵. The character "S" appended in three models emphasizes the density "stratification" in the background ICM, while "H" and "L" denote "high" and "low" resolutions. The grid resolution is given by the number of grid zones in a jet radius of r_j = 1 kpc, N_j = r_j /Δx, in the tenth column. Here, Δx is the size of grid zones. The fourth column shows the momentum injection rate or the jet thrust in Equation (10) of Paper II.

The dynamics of relativistic jets are commonly described with the jet head speed, ${u}_{\mathrm{head}}^{* }={u}_{j}\cdot \sqrt{{\eta }_{r}}/(\sqrt{{\eta }_{r}}+1)$ (Martí et al. 1997), which represents the approximate advance speed of the jet head derived from the balance between the jet ram pressure and the background pressure, and the jet crossing time given as ${t}_{\mathrm{cross}}={r}_{j}/{u}_{\mathrm{head}}^{* }$. Here, ${\eta }_{r}=({\rho }_{j}{h}_{j}{{\rm{\Gamma }}}_{j}^{2})/({\rho }_{c}{h}_{c})$ is the relativistic density contrast. The seventh and eighth columns of Table 1 show ${u}_{\mathrm{head}}^{* }$ and t_cross, respectively. The end time of simulations, t_end/t_cross, is given in the last column. We note that simulations run up to t_end/t_cross ∼ 110–200, longer than those in Paper II, so the bow shock reaches up to a few hundred kiloparsecs (see Figure 1), intending to reproduce realistic jet flows for the study of UHECR acceleration.

As in Paper II, a slow, small-angle precession with a period of τ_prec = 10 t_cross and an angle of θ_prec = 05 is added to the jet inflow velocity in order to break the rotational symmetry of the system.

The typical morphology of relativistic jets that is realized by numerical simulations may include the following features: recollimation shocks formed in the jet-spine flow, a terminal shock at the head of the jet, a cocoon of the shocked jet material flowing backward, the shocked ICM, and a bow shock that encompasses the entire jet-induced flow, (e.g., English et al. 2016; Li et al. 2018; Perucho et al. 2019; Paper II). Figure 1 shows two-dimensional (2D) slice images of the density and the z-velocity along the jet direction for all the models considered, illustrating some of these features. Similar images for the models in the uniform ICM were presented in Figures 2 and 3 of Paper II. As can be seen in these figures, the jet morphology depends on the jet power Q_j , where higher Q_j s induce more elongated jets, while lower Q_j s result in more extended, turbulent cocoons.

Zoom In Zoom Out Reset image size

The dynamics of jets in the stratified background models are overall similar to those in the uniform background models. While the latter are described in detail in Paper II, we here briefly comment the effects of density stratification, by examining the evolution of the length ${ \mathcal L }$ and width ${ \mathcal W }$ of the cocoon as a function of t/t_cross, shown in panels (a) and (b) of Figure 2; ${ \mathcal L }/{ \mathcal W }$ and ${ \mathcal W }$ for the models in the uniform ICM are shown in Figure 5 of Paper II. We find that in the low-power model of Q45-η5-S, the jet advance is a little slower in the stratified background model, owing to enhanced lateral expansion near the jet head, than in the uniform background model. In contrast, for the high-power models of Q46-η5-S and Q47-η5-S, the jet head penetrates a little faster into the density-decreasing ICM in the stratified background models. On the other hand, as the cocoon expands, although its width fluctuates differently in the different models, ${ \mathcal W }$ is on average almost identical in the models with stratified and uniform backgrounds.

Zoom In Zoom Out Reset image size

We also find that, with the higher numerical resolution in Q46-η5-H, ${ \mathcal L }$ and ${ \mathcal W }$ are a bit larger than in Q46-η5-S. In addition, as we also showed in Paper II, nonlinear structures such as shocks, turbulence, and velocity shear are better captured in higher resolution simulations. As a consequence, the MC simulations of CR acceleration would be somewhat resolution-dependent; the acceleration would be more efficient in Q46-η5-H than in Q46-η5-S (see Section 5.2 for a discussion on the resolution dependence of particle acceleration).

While the magnetic field is one of the key physical ingredients that govern the particle APs, this paper is based on RHD simulations, partly because fully relativistic magnetohydrodynamic (RMHD) simulations are more challenging and require higher resolutions (see, e.g., Martí 2019, for a review). Thus, we here adopt a set of prescriptions for the strength of the magnetic field, B, which utilize the hydrodynamic properties of simulated jet flows, such as the internal energy, the turbulent energy, and the shock speed.

In the estimation of synchrotron emission in simulated radio jets, B was often modeled assuming that the magnetic energy density is a fixed fraction of the internal energy density (see, e.g., Wilson & Scheuer 1983; Gómez et al. 1995). Following this approach, we first parameterize B with the plasma beta, β_p = P/P_B as follows:

$$\begin{eqnarray}&&{P}_{B}=\displaystyle \frac{{B}_{p}^{2}}{8\pi }=\displaystyle \frac{P}{{\beta }_{p}}.\end{eqnarray} \tag{ 3 }$$

Here, B_p denotes the magnetic field strength derived from the β_p prescription. We adopt β_p = 100 as a fiducial value and also consider β_p = 10 as a comparison case. We point out that the ICM is observed to be weakly magnetized if it has a characteristic value of β_p ∼ 100 (see, e.g., Ryu et al. 2008; Porter et al. 2015).

In turbulent flows, a magnetic field is generated via the so-called small-scale turbulent dynamo, and the magnetic energy density approaches equipartition with the kinetic energy density (see, e.g., Cho et al. 2009; Zrake & MacFadyen 2012). Considering that turbulence is ubiquitous in jet flows as demonstrated in Paper II, we introduce the dynamo-amplified field strength, B_turb, defined as:

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\mathrm{turb}}^{2}}{8\pi }\approx {{ \mathcal E }}_{\mathrm{turb}},\end{eqnarray} \tag{ 4 }$$

where ${{ \mathcal E }}_{\mathrm{turb}}$ is the kinetic energy density of the turbulent flow. Here, ${{ \mathcal E }}_{\mathrm{turb}}={{\rm{\Gamma }}}_{\mathrm{turb}}({{\rm{\Gamma }}}_{\mathrm{turb}}-1)\rho {c}^{2}$ and ${{\rm{\Gamma }}}_{\mathrm{turb}}={(1-{({u}_{\mathrm{turb}}/c)}^{2})}^{-1/2};$ the turbulent flow velocity, u_turb, is estimated by filtering the large-scale flow motions as described in Paper II (see also Section 3.3).

Furthermore, numerous shocks arise in jet flows (see Figure 4(a) below), and the magnetic field is expected to be amplified via both Bell's resonant and nonresonant CR streaming instabilities near the shocks (see, e.g., Bell 2004; Caprioli & Spitkovsky 2014a, 2014b). So we calculate B_Bell at "shock zones" (grid zones with shocks) in the simulated jet-induced flows, defined as:

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\mathrm{Bell}}^{2}}{8\pi }\approx \displaystyle \frac{3}{2}\displaystyle \frac{{u}_{s}}{c}{P}_{\mathrm{CR}}.\end{eqnarray} \tag{ 5 }$$

In our model, the CR pressure is approximated as ${P}_{\mathrm{CR}}\approx 0.1{\rho }_{1}{u}_{s}^{2}$ with the pre-shock density ρ₁ and the shock speed u_s , which reflects the DSA simulation results for nonrelativistic shocks (Caprioli & Spitkovsky 2014a).

We then take a practical, yet physically motivated approach, in which the highest estimate is chosen among the three model values:

$$\begin{eqnarray}&&B=\max ({B}_{p},{B}_{\mathrm{turb}},{B}_{\mathrm{Bell}}).\end{eqnarray} \tag{ 6 }$$

This is done in the fluid frame, and this comoving magnetic field strength is used in the calculation of the particle MFP and the mean acceleration timescales described below.

Panels (a) and (b) of Figure 3 show 2D slice images of the comoving magnetic field in the fluid frame, B, and the observed magnetic field in the simulation frame (i.e., observer frame), B_obs, for the Q46-η5-S model with β_p = 100. To obtain B_obs, B has to be Lorentz transformed, since the magnetic field is a frame-dependent quantity. Without vector information of the magnetic field, however, we approximate the magnetic field strength in the observer frame as B_obs ≈ Γ_f B, where Γ_f is the fluid Lorentz factor calculated with the fluid speed, u, in the simulation frame.

Zoom In Zoom Out Reset image size

In the background ICM, the flow is static, and hence B = B_obs = B_p ; with β_p = 100, the ICM has B_obs ≈ 1 μG in our setup, which is the typical magnetic field strength observed in the ICM (e.g., Carilli & Taylor 2002; Govoni & Feretti 2004). Our adopted model results in B_obs ∼ 10–100 μG in the backflow and the jet-spine flow, which is in a good agreement with the magnetic field strength inferred from X-ray and radio observations of radio galaxies (e.g., Begelman et al. 1984; Kataoka & Stawarz 2005; Anderson et al. 2022).

In Figure 3(c), PDFs of the comoving magnetic field B with β_p = 100 in the jet-spine flow (red) and the backflow (blue) are plotted for the three "S" models (thick lines); the PDFs of B with β_p = 10 for the Q46-η5-S model are also presented (thin solid lines). The peaks of the PDFs represent relatively quiet zones where B ≈ B_p , while the broad distributions mainly include dynamically active zones where B ≈ B_turb. With β_p = 100, B_turb is likely to be larger than B_p in the turbulent jet-spine flow and backflow, while B_Bell is dominant at shock zones. If we adopt β_p = 10, however, B_p could be larger than B_turb even in some zones of turbulent flows. The PDFs for the three "S" models with β_p = 100 show that the magnetic field is stronger in more powerful jets, since flows with higher pressure P and higher Γ_turb are produced.

Figure 2(c) shows the time evolution of the volume-averaged, comoving magnetic field strength in the cocoon, 〈B(t)〉, estimated with β_p = 100 for some models considered. As expected, 〈B(t)〉 is larger in jets with higher Q_j . Also, 〈B(t)〉 decreases in time due to both the lateral and radial expansions of the cocoon. Although nonlinear structures are better resolved in high-resolution simulations, we find that 〈B(t)〉 is not very sensitive to the resolution and is almost identical in the Q46-η5-S and Q46-η5-H models.

We note that our modeling of the magnetic field, although it is physically motivated, is still somewhat arbitrary. The resulting B affects the MFP and hence the particle acceleration, as shown in the next section. In Section 5, we compare the energy spectra of UHECRs in the Q46-η5-H jet for a specific scattering model with β_p = 100 and 10. A stronger magnetic field results in a spectrum that shifts a little to higher energies. However, we find that once B is in the range observed for radio jets, the differences due to the different magnetic field modelings would not be substantial.

Scattering of particles off underlying MHD fluctuations is a key element that governs particle transport in both spatial and momentum spaces, acceleration, confinement, and escape from the system. The important measure is the gyroradius of particles with energy E, which is given as:

$$\begin{eqnarray}&&{r}_{g}\approx \displaystyle \frac{1.1\,\mathrm{kpc}}{{Z}_{i}}\left(\displaystyle \frac{E}{1\,\mathrm{EeV}}\right){\left(\displaystyle \frac{B}{1\,\mu {\rm{G}}}\right)}^{-1},\end{eqnarray} \tag{ 7 }$$

where Z_i is the charge of the CR nuclei. The maximum energy derived from the confinement condition that r_g is equal to the radius of the acceleration system, ${ \mathcal R }$, is given as:

$$\begin{eqnarray}&&{E}_{H,{ \mathcal R }}\approx 0.9\,\mathrm{EeV}\cdot {Z}_{i}\left(\displaystyle \frac{{\text{}}B}{1\,\mu {\rm{G}}}\right)\left(\displaystyle \frac{{ \mathcal R }}{1\,\mathrm{kpc}}\right).\end{eqnarray} \tag{ 8 }$$

This geometrical condition is known as the "Hillas criterion," and ${E}_{H,{ \mathcal R }}$ is referred as the "Hillas energy" (Hillas 1984), which provides a rough estimate of the energy with which particles are confined before escaping from the system.

In general, the MFP of CRs, λ_f (E), is thought to be momentum- or energy-dependent. In a magnetized, strongly turbulent, collisionless plasma, the diffusion of particles across the magnetic field is often conjectured to follow Bohm diffusion, leading to λ_f (E) ∼ r_g (Bohm 1949). It is known that at shocks, self-generated magnetic fluctuations via various microinstabilities, such as resonant and nonresonant CR streaming instabilities, can be described by the Bohm limit, and hence λ_f (E) ∝ E (e.g., Caprioli & Spitkovsky 2014a). On the other hand, it is argued that in Kolmogorov-type turbulence, resonant scattering results in λ_f (E) ∝ E^1/3 (e.g., Stawarz & Petrosian 2008), while on scales larger than the coherence length of turbulence, nonresonant scattering might result in λ_f (E) ∝ E² (e.g., Sironi et al. 2013). Hence, for instance, in Kimura et al. (2018), the MFP was assumed to be scaled with energy as λ_f (E) ∝ E^1/3 on small scales and λ_f (E) ∝ E² on large scales in the cocoon.

We here adopt a simple prescription for the MFP:

$$\begin{eqnarray}&&{\lambda }_{f}(E)={\left(\displaystyle \frac{E}{{E}_{H,{L}_{0}}}\right)}^{\delta }{L}_{0},\end{eqnarray} \tag{ 9 }$$

where L₀ ∼ r_j is the coherence length of turbulence in our jet simulations (see Paper II) and ${E}_{H,{L}_{0}}$ is the Hillas energy at the coherence length. Then, the mean scattering time is given as τ(p) ≈ λ_f /c ∝ p^δ , where p ≈ E/c is the momentum of CRs. Here, both λ_f (E) and τ(p) are defined in the local fluid frame (i.e., scattering frame).

In the fiducial model, we adopt δ = 1/3 for $E\lt {E}_{H,{L}_{0}}$ and δ = 1 for $E\gt {E}_{H,{L}_{0}}$ in both the jet-spine flow and the backflow. If the particle is located in a shock zone, however, δ = 1 is assigned, regardless of its energy. In order to explore the dependence of the acceleration of the highest energy CRs on particle scattering (see Section 5.3 for a discussion), we consider two additional models specified in Table 2. In Model A, resonant scattering in Kolmogorov-type turbulence (δ = 1/3) is assumed for high-energy particles; in Model B, nonresonant scattering (δ = 2) is assumed for high-energy particles. Note that Model B is closest to, but slightly different from, that of Kimura et al. (2018), in which Bohm-type diffusion is adopted in the jet-spine flow. In their simple geometrical setup, the jet-cocoon system consists of an upward-moving jet-spine flow and a laterally expanding cocoon, so they focused mainly on RSA and did not include DSA.

Table 2. MFP Models

MFP:	${\lambda }_{f}={(E/{E}_{H,{L}_{0}})}^{\delta }{L}_{0}$
Bohm scattering:	δ = 1
Kolmogorov scattering:	δ = 1/3
nonresonant scattering:	δ = 2
fiducial:	$E\lt {E}_{H,{L}_{0}}$: δ = 1/3, δ = 1 (at shocks)
	$E\gt {E}_{H,{L}_{0}}$: δ = 1
Model A:	$E\lt {E}_{H,{L}_{0}}$: δ = 1/3, δ = 1 (at shocks)
	$E\gt {E}_{H,{L}_{0}}$: δ = 1/3
Model B:	$E\lt {E}_{H,{L}_{0}}$: δ = 1/3, δ = 1 (at shocks)
	$E\gt {E}_{H,{L}_{0}}$: δ = 2

Download table as:  ASCIITypeset image

Matthews et al. (2019) suggested that nonrelativistic or mildly relativistic shocks in the lobe can effectively accelerate UHECRs via DSA, while relativistic shocks such as terminal shocks and recollimation shocks would be less efficient as particle accelerators (e.g., Bell et al. 2018). The maximum energy of particles achievable at astrophysical shocks can be estimated from the condition that the diffusion length of UHECRs in the Bohm limit, l_diff ∼ λ_f (c/u_s ), should be smaller than the shock size, r_s :

$$\begin{eqnarray}&&{E}_{\mathrm{shock},\max }\approx 0.9\,\mathrm{EeV}\cdot {Z}_{i}\left(\displaystyle \frac{{\text{}}B}{1\,\mu {\rm{G}}}\right)\left(\displaystyle \frac{{\text{}}{u}_{s}}{{\text{}}c}\right)\left(\displaystyle \frac{{\text{}}{r}_{s}}{1\,\mathrm{kpc}}\right),\end{eqnarray} \tag{ 10 }$$

where u_s is the shock velocity. If r_s ∼ 1–10 kpc and B ∼ 10–100 μG in the lobe, radio jets could be potential sources of UHECRs of up to E∼10²⁰ eV.

In the test particle regime of DSA, the energy spectrum of CRs accelerated by nonrelativistic shocks takes a simple power-law form, $d{ \mathcal N }/{dE}\propto {E}^{-\sigma }$, where the slope, σ = (χ + 2)/(χ − 1), is determined solely by the compression ratio across the shock jump, χ = ρ₂/ρ₁ (e.g., Bell 1978; Drury 1983). For instance, for strong nonrelativistic shocks with M_s ≫ 1, χ = 4 and σ = 2.

The acceleration physics at relativistic shocks is more complex, and they depend on the shock speed and the magnetic field obliquity, as well as the particle scattering laws, among other parameters, in addition to shock compression (see Sironi et al. 2015, and references therein). In the test particle regime for ultra-relativistic shocks, the power-law slope approaches σ ≈ 2.2, and hence the energy spectrum tends to be steeper than that of strong nonrelativistic shocks. (e.g., Kirk et al. 2000; Keshet & Waxman 2005; Ellison et al. 2013).

On the other hand, since numerous shocks form in turbulent, jet-induced flows, DSA by multiple shocks is highly pertinent in this situation. Reacceleration by multiple, nonrelativistic shocks is known to flatten the DSA power-law spectrum to ∝p⁻³, independent of the shock compression ratio (e.g., Melrose & Pope 1993; Casse & Marcowith 2003). If the effects of other APs such as TSA and GSA (see below) are included as well, the energy spectrum of such multi-shock accelerated particles could be even flatter than E⁻¹.

The mean DSA timescale for CR protons at nonrelativistic shocks is given as:

$$\begin{eqnarray}&&{t}_{\mathrm{DSA}}\approx 3.52\times {10}^{3}\,\mathrm{yr}\,\displaystyle \frac{\chi (\chi +1)}{\chi -1}{\left(\displaystyle \frac{{u}_{s}}{c}\right)}^{-2}\left(\displaystyle \frac{E}{1\,\mathrm{EeV}}\right){\left(\displaystyle \frac{B}{1\,\mu {\rm{G}}}\right)}^{-1},\end{eqnarray} \tag{ 11 }$$

where the Bohm diffusion coefficient, ${\kappa }_{B}\,\approx {(3.13\times {10}^{22}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1})(B/1\mu {\rm{G}})}^{-1}(p/{m}_{p}c)$ is adopted (Drury 1983). Here, m_p is the proton mass. As discussed in detail in Paper II, we identify "shock zones" in the simulated jet-induced flows, as shown in Figure 4(a). The shock properties, such as M_s and u_s , are calculated from the simulation data and used to estimate t_DSA.

The velocity shear appearing at the interfaces between the jet-spine flow and the backflow and between the backflow and the shocked ICM (Figure 4(c)) is subject to Kelvin–Helmholtz instabilities, resulting in turbulence all over the jet-induced structures. In Paper II, by filtering out the bulk jet motions on scales larger than the characteristic scale of the jet, L₀ ∼ r_j , we found that the jet-spine flow and the backflow exhibit Kolmogorov power-law turbulence of ∝k^−5/3, where the solenoidal mode dominates over the compressive mode. Hence, the so-called TSA caused by incompressible turbulence is expected to operate for particles whose MFP is smaller than ∼r_j (e.g., Ohira 2013). For particles with λ_f (p) ≳ r_j , on the other hand, the RSA would become important, which will be discussed in the next subsection.

As described in Paper II, the turbulent component of flow motions is extracted as (Vazza et al. 2017):

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{\mathrm{turb}}({\boldsymbol{r}})=\displaystyle \frac{{\boldsymbol{u}}({\boldsymbol{r}})-{\left\langle {\boldsymbol{u}}({\boldsymbol{r}})\right\rangle }_{{L}_{0}}}{1-{\boldsymbol{u}}({\boldsymbol{r}})\cdot {\left\langle {\boldsymbol{u}}({\boldsymbol{r}})\right\rangle }_{{L}_{0}}/{c}^{2}},\end{eqnarray} \tag{ 12 }$$

where ${\left\langle {\boldsymbol{u}}({\boldsymbol{r}})\right\rangle }_{{L}_{0}}={\sum }_{i}{w}_{i}{{\boldsymbol{u}}}_{i}/{\sum }_{i}{w}_{i}$ is the mean of the flow velocity averaged over the cubic box of size L₀. Here, we use a simple weight function, w_i = 1, and set L₀ = r_j . This method cannot perfectly separate the turbulent component from the strong laminar component in the direction of jet propagation. So assuming that the turbulent velocity is almost isotropic, i.e., u_turb,x ≈ u_turb,y ≈ u_turb,z , the turbulence speed is approximated as:

$$\begin{eqnarray}&&| {u}_{\mathrm{turb}}| \approx {\left[\displaystyle \frac{3}{2}({u}_{\mathrm{turb},{\text{}}{\rm{x}}}^{2}+{u}_{\mathrm{turb},{\text{}}{\rm{y}}}^{2})\right]}^{1/2}.\end{eqnarray} \tag{ 13 }$$

Figure 4(b) shows a 2D slice image of ∣u_turb∣.

Ohira (2013) considered TSA in nonrelativistic incompressible turbulence. He derived analytic solutions for the momentum diffusion coefficient, D_TSA, and the acceleration timescale, t_TSA = p²/D_TSA, when the turbulence is of Kolmogorov-type. We here adopt Equation (19) of Ohira (2013) as follows:

$$\begin{eqnarray}&&{t}_{\mathrm{TSA}}\approx 2.88\times {10}^{4}\,\mathrm{yr}{\left(\displaystyle \frac{{L}_{0}/{\rm{\Gamma }}}{1\,\mathrm{kpc}}\right)}^{2/3}{\left(\displaystyle \frac{| {u}_{\mathrm{turb}}| }{c}\right)}^{-2}{\left(\displaystyle \frac{{\lambda }_{f}(E)}{1\,\mathrm{kpc}}\right)}^{1/3},\end{eqnarray} \tag{ 14 }$$

where the relativistic length contraction effect is included in the L₀ term. The energy spectrum of CRs produced by TSA depends on λ_f and the characteristics of turbulence. In general, it does not take a simple power-law form.

Once shear, Ω_shear = ∣∂u_z /∂r∣, develops, particles can be energized by encountering a velocity difference due to the shear, Δu = Ω_shear λ_f , as they are elastically scattered off MHD fluctuations frozen into the flow (see, e.g., Rieger 2019, for a review). Particles with λ_f (E) ≲ Δr (the width of the shear layer) undergo a stochastic AP inside the shear layer, which is called GSA. On the other hand, particles with λ_f (E) > Δr experience the whole velocity discontinuity by crossing the entire shear layer on each scattering, and can undergo the so-called discrete or nGSA.

As shown in Figure 4(c), Ω_shear is the largest along the interface between the jet-spine flow and the backflow, and hence GSA and nGSA operate mostly across this interface. The shear along the interface, which is relativistic, Ω_shear r_j /c ∼ 1, has a width of the order of the jet radius, Δr ∼ r_j (see also Paper II). So particles with λ_f (E) ≲ r_j are expected to be accelerated by GSA, while a fraction of high-energy particles with λ_f (E) > r_j are accelerated via nGSA.

For relativistic GSA, the acceleration timescale was estimated, for instance, by Webb et al. (2018). We adopt the timescale in their Equation (21):

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{GSA}} & = & \displaystyle \frac{15}{(4+\delta ){{\rm{\Gamma }}}_{z}^{4}{{\rm{\Omega }}}_{\mathrm{shear}}^{2}\tau (p)}\\ & \approx & 4.90\times {10}^{4}\,\mathrm{yr}\displaystyle \frac{1}{(4+\delta ){{\rm{\Gamma }}}_{z}^{4}}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{shear}}}{c/{r}_{j}}\right)}^{-2}{\left(\displaystyle \frac{{\lambda }_{f}(E)}{1\,\mathrm{kpc}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 15 }$$

where ${{\rm{\Gamma }}}_{z}={(1-{({u}_{z}/c)}^{2})}^{-1/2}$ and λ_f (E) ∝ E^δ is used. Note that particles with longer λ_f (E) experience larger velocity differences in the shear flow, so t_GSA is inversely proportional to the particle MFP.

Particles injected to the nGSA process with λ_f (E) > Δr gain energy, on average 〈ΔE/E〉 ∼ (Γ_Δ − 1) per each crossing of the velocity discontinuity, Δu, across the shear layer, where ${{\rm{\Gamma }}}_{{\rm{\Delta }}}={(1-{({\rm{\Delta }}u/c)}^{2})}^{-1/2}$ (Rieger & Duffy 2004). The mean energy gain per each cycle of crossing and recrossing the shear layer is given as $\langle {\rm{\Delta }}E/E\rangle \sim ({{\rm{\Gamma }}}_{{\rm{\Delta }}}^{2}-1)$, if the particle momentum distribution is almost isotropic. However, this could be an overestimation since velocity anisotropy could be substantial in relativistic shear flows. We take the mean acceleration timescale per cycle given in Equation (1) of Kimura et al. (2018) only as an approximate measure:

$$\begin{eqnarray}&&{t}_{\mathrm{nGSA}}\sim \zeta \displaystyle \frac{{\lambda }_{f}(p)}{c{{\rm{\Gamma }}}_{z}^{2}{\beta }_{z}^{2}},\end{eqnarray} \tag{ 16 }$$

where β_z = u_z /c and ζ ∼ 1 is a numerical factor that mainly depends on the anisotropy of the particle distribution.

Ostrowski (1998) presented an analytic solution for the momentum distribution function, f(p), for nGSA at a tangential discontinuity in the case of continuous mono-energetic injection. If there is no intrinsic scale in the system (such as the jet radius or the escape boundary size) and if particles are scattered with a mean scattering time of τ(p) ∝ p^δ , the accelerated spectrum can be represented by f(p) ∝ p^−3+δ (see their Equation (B2)), resulting in an energy spectrum of $d{ \mathcal N }/{dE}=4\pi {p}^{2}f(p)\propto {E}^{-1+\delta }$. For Bohm diffusion (δ = 1), $d{ \mathcal N }/{dE}\propto {E}^{0}$, which is much flatter than the canonical DSA power law $d{ \mathcal N }/{dE}\propto {E}^{-2}$ for strong nonrelativistic shocks.

Kimura et al. (2018) performed MC simulations for a mildly relativistic jet of Γ_j ≈ 1.4, represented by a simplified jet-cocoon system in a cylindrical configuration. Seed galactic CRs are energized through large-angle scatterings in a manner of random walks. Adopting various scattering prescriptions for λ_f (E) and a uniform magnetic field in the cocoon and the jet-spine flow, they found that nGSA produces a power-law spectrum of $d{ \mathcal N }/{dE}\propto {E}^{-1}-{E}^{0}$ for escaping CRs. The high-energy end of the spectrum above the cutoff energy decreases more gradually than an exponential.

As mentioned in the introduction, Caprioli (2015) suggested the "espresso" scenario of particle acceleration, which is conceptually similar to nGSA being described here. If a particle experiences a "one-shoot boost" by crossing and recrossing the shear layer in an ultra-relativistic jet with Γ_j ≫ 1, the energy is enhanced by a factor of $\sim {{\rm{\Gamma }}}_{j}^{2}$ per cycle. In later studies, Mbarek & Caprioli (2019, 2021) performed MC simulations, in which seed CRs are injected into a self-consistent jet configuration from MHD simulations and their trajectories are followed. They found that particles could be accelerated to become UHECRs of E ≳ 10²⁰ eV via one or two espresso shots.

To assess the relative importance of different APs, we compare the acceleration timescales given in Equations (11), (14), (15), and (16) in Figure 5. The adopted characteristic parameters are specified in the figure caption. For particles with E ≲ 10¹⁸ eV, t_DSA is the shorter than t_TSA and t_GSA, and DSA would be the dominant AP. For higher energy particles with E ≳ 10¹⁸ eV, GSA would become important. As noted above, only a small fraction of particles, which are energized via other processes to have λ_f (E) > r_j , could cross the entire shear layer and are injected to the nGSA process. Hence, although t_nGSA is shortest even at low energies, nGSA operates only for the highest energy CRs. TSA, on the other hand, would be only marginally important around E ∼ 10¹⁸ eV.

Zoom In Zoom Out Reset image size

If the simple condition t_DSA ≈ t_GSA is adopted, the transition energy above which GSA becomes faster than DSA is roughly E_trans ≈ 4 EeV ${(\langle {u}_{s}\rangle /c)(\langle B\rangle /1\mu {\rm{G}})\langle {{\rm{\Gamma }}}_{j}{\rangle }^{-2}(\langle {{\rm{\Omega }}}_{\mathrm{shear}}\rangle {r}_{j}/c)}^{-1}$. For the jet models considered here, the volume-averaged quantities range approximately as follows: (〈u_s 〉/c) ∼ 0.3–0.5, 〈Γ_j 〉 ∼ 1.3–5, (〈Ω_shear〉r_j /c) ∼ 1–1.5, and 〈B〉 ∼ 3–15 μG (see also Paper II). All of them increase slowly with Q_j , resulting in E_trans ∼ 1 EeV. In Section 5.1, we will show that indeed, at low energies, DSA seems to be important, while at high energies, GSA and nGSA are dominant, in MC simulations.

In the early development stages of radio jets, the maximum energy, ${E}_{\max }$, that CRs can achieve is limited by age. In later stages, ${E}_{\max }$ is expected to be limited by the size of the cocoon where CRs are confined. Since the smallest acceleration timescales in Figure 5 are much shorter than the duration of our simulations, t_end ≈ 10⁶−10⁷ yr (see Table 1), the size-limited ${E}_{\max }$ would be relevant.

Most CRs are expected to escape diffusively from the cocoon. So the size-limited ${E}_{\max }$ can be estimated from the condition that λ_f (E) is equal to the radius of the cocoon. With the lateral width of the cocoon, ${ \mathcal W }$, shown in Figure 2(b), ${\lambda }_{f}(E)\sim { \mathcal W }/2$ gives:

$$\begin{eqnarray}\begin{array}{rcl}{E}_{\max } & \approx & {E}_{H,{L}_{0}}\left(\displaystyle \frac{{ \mathcal W }}{2{L}_{0}}\right)\\ & \approx & 0.9\,\mathrm{EeV}\cdot {Z}_{i}\left(\displaystyle \frac{B}{1\,\mu {\rm{G}}}\right)\left(\displaystyle \frac{{r}_{j}}{1\,\mathrm{kpc}}\right)\left(\displaystyle \frac{{ \mathcal W }}{2{r}_{j}}\right),\end{array}\end{eqnarray} \tag{ 17 }$$

where again L₀ = r_j is used. For $\langle B\rangle { \mathcal W }\sim 200$ (see Figure 2(d)), ${E}_{\max }$ is estimated to be ∼200 EeV for protons. Note that the above size-limited ${E}_{\max }$ does not directly depend on Γ_j ; hence, ${E}_{\max }$ could be lower even in more powerful jet models with higher Γ_j , if ${ \mathcal W }$ is smaller.

On the other hand, the length of the cocoon, ${ \mathcal L }$, is greater than ${ \mathcal W }$, as shown in Figure 2(a). So CRs even with $E\gtrsim {E}_{\max }$ could be confined in the cocoon, and they may escape through the longitudinal direction. In addition, a small fraction of CRs can be boosted to much higher energies via nGSA, and their escape may not be described as a diffusive process. In Section 5.2, we will discuss how such anisotropic escape affects the high-energy end of the spectrum of UHECRs.

To track the acceleration of CRs in simulated jet-induced flows, we perform MC simulations, where the trajectories of CRs are integrated, according to the following recipes: (1) in the jet simulations described in Section 2, the snapshot data of jet-flow quantities are stored with a time interval of Δt_s = (1/3)t_cross; (2) at every Δt_s , 200 seed CRs are injected from the jet nozzle with energy spectrum, $d{{ \mathcal N }}_{\mathrm{inj}}/{{dE}}_{\mathrm{inj}}\propto {E}_{\mathrm{inj}}^{-2.7}$, for E_inj = (0.01–1) PeV. Note that −2.7 is the slope of the Galactic CR energy spectrum (e.g., Thoudam et al. 2016); (3) CR particles are assumed to be scattered elastically off small-scale MHD fluctuations that are frozen in the background flow in the "restricted" random walk scheme (see Section 4.2); (4) large-angle scattering is imposed by hand, adopting the prescription for the energy-dependent MFP, λ_f (E), in Equation (9) with the magnetic field model described in Section 2.4; (5) the probability of particle displacement at each scattering is assumed to obey an exponential distribution with λ_f (E); (6) at every elastic scattering event, the Lorentz transformation from the local fluid frame to the simulation frame is performed, which results in the energy change, ΔE (either positive or negative); (7) the energy evolution of CR particles is calculated along the evolution of jet-induced flows using the snapshot data; and (8) the information on the particles escaping from the computational box is stored and used to analyze the properties of UHECRs, such as the energy spectrum.

Seed CRs with E_inj have gyroradii that are much smaller than the grid size, i.e., r_g ≪ Δx = 0.1–0.33 kpc (see the tenth column of Table 1), and hence go through scatterings within a grid zone. For the calculation of the energy change due to such subgrid scatterings, we employ variation of the 3D fluid velocity inside a grid cell, which is approximated using trilinear interpolation with the fluid velocities along the direction of particle trajectory. These subgrid scatterings would be applied mostly to CRs of E < 10¹⁸ eV (see Equation (8)).

In the calculation of conventional random walk transport, the components of the scattering angle, (δ μ, δ ϕ), are chosen randomly in the ranges −1 ≤ δ μ ≤ + 1 and 0 ≤ δ ϕ ≤ 2π, respectively, where $\mu =\cos \theta$. In real jet flows, however, magnetic field fluctuations may not be large enough to scatter high-energy particles fully isotropically with λ_f (E) ≳ L₀. Thus, to mimic roughly pitch-angle diffusion without knowing the magnetic field configuration, we adopt a random walk scheme, in which the scattering angle with respect to the incident direction is chosen from the "restricted" range of $\delta {\mu }_{\max }\leqslant \delta \mu \leqslant +1$. In addition, we employ an energy-dependence, with which the scattering angle with respect to the incident direction is forced to be smaller than ∼π[L₀/λ_f (E)]. Without a prior knowledge of the turbulent nature of magnetic field fluctuations, it may provide a crude model to account for pitch-angle scattering in an energy-dependent way.

In our random walk scheme, the maximum value of scattering angle is modeled specifically as:

$$\begin{eqnarray}&&\delta {\theta }_{\max }\approx \pi \cdot \min \left[\psi \left(\displaystyle \frac{{L}_{0}}{{\lambda }_{f}}\right),1\right].\end{eqnarray} \tag{ 18 }$$

Here, ψ ≲ 1 is a free parameter that is devised to reflect the strength of magnetic field fluctuations. Again, L₀ ∼ r_j is assumed for the coherence scale of turbulence. For high-energy particles with λ_f (E) ≫ L₀, this prescription leads to forward-beamed scattering with $\delta {\theta }_{\max }\ll 1$. For low-energy particles with λ_f (E) ≪ L₀, by contrast, it results in isotropic scattering with $\delta {\theta }_{\max }\approx \pi$. Adopting this scheme tends to reduce the energy gain of high-energy particles, especially for nGSA near the jet–backflow interface. We take ψ = 1 as a fiducial value, and also present the isotropic scattering case (ψ → ∞ ) for a demonstration of model dependence.

As illustrated in Figure 5, in jet-induced flows, DSA and TSA would be important for low-energy particles, while GSA and nGSA would become more important for higher energy particles near the jet–backflow interface. In MC simulations, however, in each scattering, the change in the particle energy often involves a combination of the different processes.

In an effort to evaluate the relative importance of the APs, we estimate the acceleration timescales, t_TSA and t_GSA, at each scattering. If the particle crosses a single or multiple shock zones in the displacement after scattering, t_DSA is also calculated. We then attempt to identify the primary acceleration process (PAP) by determining the shortest acceleration timescale among the two or three timescales. If DSA is chosen as the PAP, the scattering event is tagged as "shock;" if TSA is chosen, it is tagged as "turbulence;" if GSA is chosen, it is tagged as "shear." We do not include t_nGSA in the PAP selection, since only a very small fraction of high-energy particles undergo nGSA. Through this crude evaluation, we will see that DSA is indeed the PAP for E ≲ 1 EeV, while GSA becomes dominant above EeV, as presented in the next section.

We begin with a discussion on the relative contributions of different APs. For each scattering event, a fraction of ΔE is assigned to each AP, according to the weight function, ${\xi }_{\mathrm{AP}}={t}_{\mathrm{AP}}^{-1}/{\sum }_{\mathrm{AP}}^{}{t}_{\mathrm{AP}}^{-1}$, where the summation includes the three APs as described in Section 4.3. For all particles escaping from the system until t = t_end, whose final energies lay in the logarithmic bin of $[\mathrm{log}E,\mathrm{log}E+d\mathrm{log}E]$, the contributions of ξ_APΔE are summed to make the cumulative energy gains, ${{\mathscr{E}}}_{\mathrm{AP}}(E)$, for each AP.

Figure 6 shows the fraction, ${{\mathscr{E}}}_{\mathrm{AP}}(E)/{{\mathscr{E}}}_{\mathrm{tot}}(E)$, for the three APs as a function of the particle energy E. Here, AP stands for "shock" (red), "turbulence" (green), and "shear" (blue), and ${{\mathscr{E}}}_{\mathrm{tot}}(E)={{\mathscr{E}}}_{\mathrm{shock}}(E)+{{\mathscr{E}}}_{\mathrm{turb}}(E)+{{\mathscr{E}}}_{\mathrm{shear}}(E)$. While these fractions should be only rough estimates, the figure confirms once again that for E ≲ 1 EeV, particles are energized mainly by DSA, whereas GSA/nGSA become increasingly important above 1 EeV. TSA makes only a supplementary contribution.

Zoom In Zoom Out Reset image size

We note in Figure 6 that for E ≲ 1 EeV, the ratio of DSA to GSA/nGSA contributions, ${{\mathscr{E}}}_{\mathrm{shock}}/{{\mathscr{E}}}_{\mathrm{shear}}$, is higher for the jet models with higher Q_j . This seems contradictory to the simple expectation that GSA and nGSA would become more significant with higher Q_j (higher Γ_j ). In fact, the MFP is given as λ_f ∝ (E/B)^δ in Equation (9) and B scales approximately as $\propto {Q}_{j}^{0.3}$ (see Figure 2(c)) in our models, reducing the probability of GSA and nGSA across the shear layer, especially for low-energy particles. So the relative importance of the different APs illustrated in the figure would be regarded as being specific to the various models employed here.

We next calculate the total energy gains due to different APs, ${\tilde{{\mathscr{E}}}}_{\mathrm{AP}}=\sum {{\mathscr{E}}}_{\mathrm{AP}}(E)$, summed over all the particle energies. Table 3 lists the fractions, ${\tilde{{\mathscr{E}}}}_{\mathrm{AP}}/{\tilde{{\mathscr{E}}}}_{\mathrm{tot}}$, where ${\tilde{{\mathscr{E}}}}_{\mathrm{tot}}={\tilde{{\mathscr{E}}}}_{\mathrm{shock}}+{\tilde{{\mathscr{E}}}}_{\mathrm{turb}}+{\tilde{{\mathscr{E}}}}_{\mathrm{shear}}$. Overall, shear acceleration (including both GSA and nGSA) is the dominant process, contributing to ∼85% of the energization of UHECRs, while DSA and TSA together generate only ∼15% of the energy. A noted point is that although TSA is subdominant in the whole energy range, its total contribution is larger than DSA. In addition, the contributions of DSA and TSA tend to be a bit larger for less powerful jet models; this is a consequence of more extended cocoons filled with shocks and turbulence in less powerful jets.

Table 3. Energy Gains via the Different APs

Model Name	$\tfrac{{\tilde{{\mathscr{E}}}}_{\mathrm{shock}}}{{\tilde{{\mathscr{E}}}}_{\mathrm{tot}}}$(%)	$\tfrac{{\tilde{{\mathscr{E}}}}_{\mathrm{turb}}}{{\tilde{{\mathscr{E}}}}_{\mathrm{tot}}}$(%)	$\tfrac{{\tilde{{\mathscr{E}}}}_{\mathrm{shear}}}{{\tilde{{\mathscr{E}}}}_{\mathrm{tot}}}$(%)
Q45-η5-S	3.2	13.6	83.2
Q46-η5-S	2.2	12.5	85.3
Q47-η5-S	1.9	12.2	85.9

Download table as:  ASCIITypeset image

We next present the energy spectrum of escaping particles. The MC simulations for GSA/nGSA by Ostrowski (1998) and Kimura et al. (2018) showed that the energy spectrum behaves as ∝E⁻¹ − E⁰ below the "break energy," E_break, while it could be approximated by another steeper power law above E_break, instead of an exponential cutoff. ⁵ Here, E_break is introduced to designate the energy above which the spectrum rolls over from one power law to another power law, whereas Ostrowski (1998) and Kimura et al. (2018) used different terms.

As described in Section 4.1, the underlying flow profile is updated and seed particles of E_inj = 0.01−1 PeV are injected at every Δt_s = (1/3)t_cross. Hence, the energy spectrum should be time-dependent. We calculate the time-integrated, cumulative energy spectrum of all particles escaping from the system up to a given time t, $d{ \mathcal N }(t)/{dE}$. Note that the total number of escaped particles, ${{ \mathcal N }}_{\mathrm{tot}}=\int (d{ \mathcal N }/{dE}){dE}$, increases with time. As shown in Figures 5 and 6, the dominant AP switches from DSA to GSA/nGSA at E ∼ 1 EeV. Thus, as the seed population of $\propto {E}_{\mathrm{inj}}^{-2.7}$ is accelerated, the resulting spectrum is expected to flatten roughly from a DSA power law of E⁻² at early stages to a GSA/nGSA power law of E⁻¹ − E⁰.

Figure 7 shows the evolution of the time-integrated energy spectrum, $E[d{ \mathcal N }(t)/{dE}]$, for the three "S" models. Below E_break, the power-law portion continuously hardens over time, approaching $d{ \mathcal N }/{dE}\propto {E}^{-0.5}$. At early epochs it is somewhat flatter than the DSA spectrum, because TSA and GSA also operate even at early stages. The time-asymptotic spectrum, $d{ \mathcal N }/{dE}\propto {E}^{-0.5}$, is in agreement with the previous works of Ostrowski (1998) and Kimura et al. (2018), indicating that GSA/nGSA would be the dominant energization process for particles around E_break at late stages. The break shifts to higher energies with time during early stages, which is consistent with the age-limited maximum energy. It gradually approaches the size-limited maximum energy at late epochs (see below). For E > E_break, the spectrum becomes steeper as time goes on, but it is still not as steep as the exponential drop even at t_end. Overall, the time-integrated spectrum asymptotically saturates by the end of the simulations in all the models.

Zoom In Zoom Out Reset image size

Because acceleration is expected to be size-limited at late stages, ${E}_{\max }$ in Equation (17) is pertinent here, and E_break in the time-asymptotic spectrum would be similar to ${E}_{\max }$. The size-limited, maximum energy, ${E}_{\max }$, scales as $\propto B{ \mathcal W }$. Figure 2 shows that the cocoon width, ${ \mathcal W }$, increases in time, while the volume-averaged magnetic field strength, 〈B〉, decreases owing to the decreasing pressure in the laterally expanding cocoon. As a result, the value of $\langle B\rangle { \mathcal W }$ increases in time during the early stage and later approaches a time-asymptotic value for t/t_cross ≳ 50. The break of the energy spectrum in Figure 7 reflects these time-dependent behaviors of ${E}_{\max }(t)\propto \langle B\rangle { \mathcal W }$. In addition, the time-asymptotic value of $\langle B\rangle { \mathcal W }$ scales roughly as $\propto {Q}_{j}^{1/3}$ (Figure 2(d)), and so does ${E}_{\max }$, for the model parameters under consideration here. This is due to the fact that while ${ \mathcal W }$ does not differ much in the different jet models (see Figures 1 and 2(b)), 〈B〉 is larger in higher power jets (Figure 2(b)).

As described in Section 2.3, the jet power Q_j is the primary parameter that determines the properties of jet-induced flows, such as the Lorentz factor of the jet flow, Γ_j , the cocoon's shape, and the associated nonlinear structures, which in turn govern the ensuing particle acceleration via DSA, TSA, and RSA. Figure 8(a) shows the time-asymptotic energy spectrum, $E[d{ \mathcal N }(t)/{dE}]$ at t_end, for the models with different Q_j . The spectrum shifts to higher energies for higher Q_j , as expected. Otherwise, the overall shape of the spectrum is similar, except at the low-energy part of E ≲ 1 EeV where the spectra differ due to different energization histories via DSA and TSA.

Zoom In Zoom Out Reset image size

To quantify the jet-power dependence, we attempt to fit the time-asymptotic spectrum in Figure 8(a) to a functional form. As stated above, $d{ \mathcal N }/{dE}$ is close to ∝ E^−0.5 for E < E_break, while it drops roughly as another steeper power law for E > E_break. Hence, instead of an exponential function, we employ the following double power-law form:

$$\begin{eqnarray}&&\displaystyle \frac{{Ed}{ \mathcal N }}{{dE}}\propto {\left({\left(\displaystyle \frac{E}{{E}_{\mathrm{break}}}\right)}^{-a}+{\left(\displaystyle \frac{E}{{E}_{\mathrm{break}}}\right)}^{b}\right)}^{-1},\end{eqnarray} \tag{ 19 }$$

where a, b, and E_break are the fitting parameters. Table 4 lists these three fitting parameters and ${E}_{\max }$ in Equation (17) estimated with the time-asymptotic values of $\langle B\rangle { \mathcal W }$. Indeed, E_break is quite similar to ${E}_{\max }$. In the figure, the fitted double power-law spectra are overlaid with dotted lines. We also plot ${Ed}{ \mathcal N }/{dE}\propto {(E/{E}_{\mathrm{break}})}^{a}\exp (-E/{E}_{\mathrm{break}})$ for the Q45 model (dashed red) in the figure, to illustrate the difference between the power-low drop and the exponential cutoff beyond E_break.

Table 4. Fitting Parameters

Model Name	a	b	Ebreak(eV)	${E}_{\max }$(eV) a
Q45-η5-S	0.59	1.64	4.5E19	7.0E19
Q46-η5-S	0.51	1.58	1.3E20	1.6E20
Q47-η5-S	0.47	1.60	2.2E20	2.8E20

Note.

^a ${E}_{\max }$ is calculated with Equation (17), adopting the time-asymptotic values of $\langle B\rangle { \mathcal W }$.

Download table as:  ASCIITypeset image

We note that E_break shifts to higher values for higher Q_j . On the other hand, the power-law slopes, a and b, show only a weak dependence on Q_j . On average, a ∼ 0.5, so $d{ \mathcal N }/{dE}\propto {E}^{-0.5}$ below E_break, as expected.

For E > E_break, the power-law slope is, on average, close to b ∼ −1.6, meaning $d{ \mathcal N }/{dE}\propto {E}^{-2.6}$. To comprehend this slope, we examine the trajectories of particles and the ensuing energization in the MC simulations, focusing on RSA at the shear interface between the jet-spine flow and the backflow. As discussed in Section 3.5, most of particles are incrementally accelerated via GSA and other processes, whereas a small fraction of high-energy particles with λ_f (E) > r_j could be boosted in energy by a factor of $({{\rm{\Gamma }}}_{{\rm{\Delta }}}^{2}-1)$ via nGSA if they cross and recross the shear interface. Typically, CRs with E ≳ E_break cross the shear interface more than once, before they escape from the jet to the ICM. In particular, energization of the highest energy particles seems to be governed by the experience of nGSA episodes. Based on this picture, we categorize escaping particles roughly into three cases, as illustrated in Figure 9. In CASE 1, particles gain energy mainly via GSA and other processes and exit the cocoon to the ICM without experiencing a boost by nGSA. In CASE 2, after small incremental accelerations, particles cross into the jet-spine, and then cross out of the jet-spine into the cocoon, resulting in an $\sim {{\rm{\Gamma }}}_{{\rm{\Delta }}}^{2}$ boost via nGSA. They are confined in the cocoon, before they exit to the ICM. CASE 3 is the same as CASE 2, except that particles cross out of the jet-spine and exit directly to the ICM.

Zoom In Zoom Out Reset image size

Figure 9 shows the trajectories and energization of three sample particles, one from each case, as a function of time since injection. In the bottom panels, the red and blue lines follow the energy changes of the particles in the simulation and fluid frames, respectively. The two big jumps in the fluid frame energy (blue lines) for CASE 2 (at t − t_inj ≈ 0.19 and 0.27) and CASE 3 (at t − t_inj ≈ 0.21 and 0.32) are a consequence of nGSA. Given that the jet-spine flow at the big-jump scattering points has Γ_f ∼ several, the energy gains in the jumps match the expectation due to scatterings into or out of the jet-spine flow. As the result of the energy boosts via nGSA, the particles of CASE 2 and CASE 3 reach well above 10²⁰ eV. In contrast, the CASE 1 particle, which does not experience nGSA, fails to reach a very high energy.

In Figure 10, we plot the time-asymptotic energy spectra for the particles of CASE 1, CASE 2, and CASE 3, separately, in the case of the Q46-η5-H model. The total number of particles for each category is ${{ \mathcal N }}_{\mathrm{tot}}(\mathrm{CASE}1)\gg {{ \mathcal N }}_{\mathrm{tot}}(\mathrm{CASE}2)\gg {{ \mathcal N }}_{\mathrm{tot}}(\mathrm{CASE}3)$. In CASE 1, particles are confined by the Hillas condition before they escape from the cocoon, so the spectrum (red) peaks at ${E}_{\max }$, above which it drops almost exponentially.

Zoom In Zoom Out Reset image size

In CASE 2, particles that have already experienced nGSA scatter in the cocoon before they escape to the ICM. Note that the cocoon has an elongated shape (see Figure 1), so particles even with $E\gtrsim {E}_{\max }$ can be confined lengthwise along the vertical direction, as shown in the top-middle panel of Figure 9. The spectrum of the particles escaping from a cylinder with a finite radius and infinite length after isotropic scattering, is given as $d{ \mathcal N }/{dE}\propto {E}^{-2}$. The break in the spectrum of CASE 2 particles (green) shifts to a higher energy due to nGSA, and follows a power-law distribution of slope ∼ −1.3, or $d{ \mathcal N }/{dE}\propto {E}^{-2.3}$ above ${E}_{\max }$. This is somewhat steeper than expected for the idealized setup, since the cocoon is not an infinitely stretched cylinder. The CASE 2 spectrum is the most dominant component in the high-energy part of the total spectrum. So the CASE 2 scenario would explain the power-law spectrum above E_break.

CASE 3 particles are relatively rare, and hence their energy spectrum (blue) may not be accurately realized in our simulations. Nevertheless, the break in the spectrum shifts to a higher energy by about an order of magnitude, compared to that of CASE 1. Considering that one cycle of nGSA produces an $\sim {{\rm{\Gamma }}}_{f}^{2}$ boost in energy and the average Lorentz factor of the jet-spine flow is $\left\langle {{\rm{\Gamma }}}_{f}\right\rangle \approx 3$ (see the insert in Figure 10), the shift of the break is well explained as a consequence of nGSA.

The combined spectrum of CASE 1, CASE 2, and CASE 3 results in a slope of ∼−1.6 above the break, which is a bit steeper than that of the CASE 2 spectrum, as demonstrated in Figure 10.

In Figure 8(b), we compare the spectra of different resolution models for the Q46 jet. Compared to the fiducial case of Q46-η5-S (green), the spectrum of the higher resolution model, Q46-η5-H (cyan), shifts to higher energies, whereas the spectrum of the lower resolution model, Q46-η5-L (purple), shifts to lower energies. As mentioned in Section 2.3, with a higher grid resolution, more significant nonlinear structures are induced; then, both DSA and TSA are more efficient owing to more frequent shocks and better developed turbulence in the cocoon. On the other hand, although E_break is slightly higher in higher resolution models, the resolution dependence seems not to be large.

The MC simulation results presented above are based on a number of modelings, such as those for the magnetic field strength in the jet-induced flows, the MFP of CR particles, and the random walk scheme. In this subsection, we briefly examine the effects of those modelings on the energy spectrum of escaping CRs. For that purpose, we use the Q46-η5-H jet.

We first examine the dependence on the magnetic field model described in Section 2.4. In relatively quiet flows, the magnetic field is likely to be prescribed by Equation (3); therein, B_p is stronger for lower β_p , which in turn may lead to more efficient particle acceleration. In contrast, in regions with well-developed turbulence, B_turb would dominate and the magnetic field is unaffected by the adopted value of β_p . In Figure 11, the spectrum for the fiducial case (β_p = 100, black solid line) is compared to that for the β_p = 10 case (dotted–dashed). The figure demonstrates that adopting a stronger magnetic field with smaller β_p leads to slightly more efficient particle acceleration, but the difference is not large in the range of the magnetic field strengths we consider.

Zoom In Zoom Out Reset image size

As pointed out in Section 3.1, the nature of magnetic turbulence and also the diffusion/scattering of CR particles in jet-induced flows are not completely understood, especially on scales larger than the coherence length of turbulence. Hence, we consider the three MFP models listed in Table 2. In Figure 11, the energy spectra for the three models are displayed. Below ${E}_{H,{L}_{0}}\sim 1$ EeV, with the same λ_f , the three spectra are basically identical. On the other hand, for E > 1 EeV, the spectra differ for the different MFP models. In Model A with δ = 1/3 (smaller λ_f ), high-energy CRs go through more scatterings and tend to achieve higher energies, whereas in Model B with δ = 2 (larger λ_f ), it works in the opposite way. Hence, the Model A spectrum shifts to higher energies and is also harder with a flatter slope above the break. In Model B, the spectrum shifts slightly to lower energies, while the slope above the break is nearly the same as in the fiducial model.

In Section 4.2, we introduced a restricted random walk scenario through Equation (18). Here, we examine the dependence on this model by comparing the spectrum with isotropic scattering (ψ → ∞) to the spectrum with restricted scattering (ψ = 1, the fiducial case) in Figure 11. As noted, the restricted, forward-beamed scattering reduces the frequency of particles crossing across the shear interface, and hence suppresses the efficiency of nGSA. Hence, the spectrum of the fully isotropic scattering case shifts to higher energies, compared to that of the fiducial model, as expected.