Model Spectrum of Ultrahigh-energy Cosmic Rays Accelerated in FR-I Radio Galaxy Jets

Nearby radio galaxies (RGs) of Fanaroff–Riley Class I (FR-I) are considered possible sites for the production of observed ultrahigh-energy cosmic rays (UHECRs). Among those, some exhibit blazar-like inner jets, while others display plume-like structures. We reproduce the flow dynamics of FR-I jets using relativistic hydrodynamic simulations. Subsequently, we track the transport and energization of cosmic ray (CR) particles within the simulated jet flows using Monte Carlo simulations. The key determinant of flow dynamics is the mean Lorentz factor of the jet-spine flow, 〈Γ〉spine. When 〈Γ〉spine ≳ 6, the jet spine remains almost unimpeded, but for 〈Γ〉spine ≲ 3, substantial jet deceleration occurs. CRs gain energy mainly through diffusive shock acceleration for E ≲ 1 EeV and shear acceleration for E ≳ 1 EeV. The time-asymptotic energy spectrum of CRs escaping from the jet can be modeled by a double power law, transitioning from ∼E −0.6 to ∼E −2.6 around a break energy, E break, with an exponential cutoff at Ebreak〈Γ〉spine2 . E break is limited either by the Hillas confinement condition or by particle escape from the cocoon via fast spatial diffusion. The spectral slopes primarily arise from multiple episodes of shock and relativistic shear acceleration, and the confinement–escape processes within the cocoon. The exponential cutoff is determined by nongradual shear acceleration that boosts the energy of high-energy CRs by a factor of ∼〈Γ〉spine2 . We suggest that the model spectrum derived in this work could be employed to investigate the contribution of RGs to the observed population of UHECRs.


INTRODUCTION
Although the origin of ultra-high-energy cosmic rays (UHECRs) with E > 1 EeV remains somewhat elusive, relativistic kpc-scale jets of radio galaxies (RGs) are considered to be the most promising source candidates (see, e.g., Rieger 2019;Matthews et al. 2020, for reviews).RG jets are commonly classified into two Fanaroff-Riley (FR) types based on their morphological characteristics (Fanaroff & Riley 1974): center-brightened FR-I and edge-brightened FR-II.In general, FR-I jets are observed to have lower radio luminosities and are therefore known to have lower kinetic powers compared to FR-II jets (e.g., Godfrey & Shabala 2013).However, the morphological manifestation of RG jets looks more complicated (e.g., Ghisellini & Celotti 2001;Mingo et al. 2019); it seems to be related to the flow dynamics of jet, which depends not only on the jet power but also on the interaction with the ambient medium (see, e.g., Hardcastle & Croston 2020, for a review).
With higher power, on average, individual FR-II jets are expected to produce a larger amount of UHECRs than individual FR-I jets.However, lower-power FR-I jets are more common in the universe (e.g., Hardcastle et al. 2019;Mingo et al. 2019).Especially in the local universe within the so-called GZK horizon of ∼ 100 Mpc Seo et al. (Greisen 1966;Zatsepin & Kuz'min 1966), virtually all RGs are the FR-I type (e.g., van Velzen et al. 2012;Rachen & Eichmann 2019).Among those nearby FR-I jets, several, including Virgo A (M87) and Centaurus A, exhibit blazar-like features.The inner jet of Virgo A is one-sided, hinting at the existence of bulk relativistic motions; superluminal motions found in optical and X-ray observations indicate the mean Lorentz factor of the jet-spine flow, ⟨Γ⟩ spine ≳ 6 (e.g., Biretta et al. 1999;Snios et al. 2019a).The jet of Centaurus A is also one-sided, but the flow looks mildly relativistic with ⟨Γ⟩ spine ∼ 1.2 − 2 (e.g., Wykes et al. 2019;Snios et al. 2019b).Meanwhile, Fornax A exhibits two inner jets embedded in the core of the host galaxy, implying recent low-power jet activity, alongside large diffuse lobes in the outer region (e.g., Geldzahler & Fomalont 1984;Maccagni et al. 2020).Hence, comprehending the variety of FR-I jet structures and then understanding the processes of UHECR acceleration within jet-induced flows are crucial for modeling the observed UHECRs originating from nearby RGs.
In Seo et al. (2023) (hereafter Paper I), we studied these acceleration processes in FR-II jets by performing relativistic hydrodynamic (RHD) simulations for the flow dynamics of jets, and Monte Carlo (MC) simulations for the transport, scattering, and energy change of CR particles injected into simulated jet flows.The magnetic field is one of the key ingredients that govern physical processes in jet-induced flows; for instance, CR acceleration and synchrotron emission rely on its distribution, strength, and fluctuations.Since our numerical approach is only hydrodynamic with limited numerical resolutions, in paper I, we adopted several physically motivated prescriptions, including magnetohydrodynamic (MHD) turbulence, to model the distribution and fluctuations of magnetic fields.
Both theoretical and observational studies suggest that magnetic fields play a critical role in driving relativistic jet flow from the central engine; helical magnetic fields look to be dynamically important in inner parsec-scale jets (e.g., Pudritz et al. 2012).On the contrary, magnetic fields are not dominant in the dynamics of the jet-spine flow and lobes on larger scales ranging from a few to hundreds of kpc (e.g., Begelman et al. 1984;Laing & Bridle 2014).This provides certain justification for our hydrodynamic approach to reproducing the jet-induced flow.In addition, even in the most advanced relativistic magnetohydrodynamic (RMHD) simulations, it is extremely challenging to fully resolve MHD turbulence on all the important scales (see, e.g., Martí 2019).Thus, it would be practical, yet scientifically justifiable, to adopt models for the distribution and fluctuations of magnetic fields and apply particle scattering laws, assuming a random walk, in MC simulations that follow the acceleration of CRs in RGs (e.g., Ellison et al. 1990;Ostrowski 1998;Kimura et al. 2018).
With these caveats, in Paper I, we suggested that for E ≲ 1 EeV, CRs gain energy mostly through DSA.But for E ≳ 1 EeV, RSA becomes more important, which leads to CRs well above 10 20 eV.TSA makes a relatively minor contribution over the entire energy range (see Figure 5 below).The time-asymptotic energy spectrum of CRs escaping from FR-II jets was approximated by the double-power law: where s 1 ≈ −0.5 and s 2 ≈ −2.6.The break energy is given by the geometrical constraint known as the "Hillas criterion".
In this paper, we investigate the CR acceleration in FR-I jets.The flow structures in FR-I jets differ from those in FR-II jets (see, e.g., Perucho et al. 2014).For example, with lower powers, in FR-I jets, ⟨Γ⟩ spine would be lower, resulting in less relativistic velocity shear in jet-induced flows.Thus, GSA and nGSA may not be as efficient as in FR-II jets.Here, we employ the same numerical approaches as in Paper I to reproduce FR-I jets and derive the energy spectrum of CRs produced at FR-I jets.We then propose a model energy spectrum for UHECRs from FR-I RGs.
The paper is organized as follows.In Section 2 and 3, three-dimensional jet simulations and CR transport simulations are described, respectively.The results are presented in Section 4, followed by a brief summary in Section 5. We consider relativistic jets that propagate through either the uniform or stratified background medium.The jets can be specified with the parameters of jet inflow, such as the power, Q j , the density, ρ j , the pressure, P j , and the radius, r j .Then, the bulk Lorentz factor, Γ j = 1/ 1 − u 2 j /c 2 , is related to the jet power as Q j = πr 2 j u j Γ 2 j ρ j h j − Γ j ρ j c 2 , where u j , h j = (e j + P j )/ρ j , e j , and c are the jet velocity, the specific enthalpy, the sum of the internal and rest-mass energy densities, and the speed of light, respectively.The jet flow is injected through a circular nozzle located in the bottom (z = 0) center of the simulation box.Intending to reproduce the outer part of FR-I RGs where the jet is kinetically dominated rather than magnetically dominated, the bottom corresponds to a surface that is tens to hundreds of parsecs away from the central black hole.
Table 1 lists the model parameters of our simulations.The first column shows the model name; the numbers after Q and r are the exponent of Q j and r j in units of pc, respectively.The models with Q j = 3.5 × 10 42 − 3.5 × 10 45 erg s −1 and r j = 10 − 100 pc represent FR-I jets (see Godfrey & Shabala 2013) that propagate up to ∼ 10 kpc.In addition, Q46-r1000, which is the Q46-η5-H model in Paper I, is included as a FR-II jet to be compared with FR-I jets.
For FR-II jets in Paper I, the stratified, isothermal intracluster medium (ICM) was chosen as the background (see in the paper for the parameters).On the other hand, for FR-I jets in this paper, the background medium appropriate for galactic halos is employed; initially, it is assumed to be isothermal.In the models with r j = 10 pc, the jet propagates up to l jet ∼ 1 kpc, which is smaller than the size of typical galaxy halos.Here, l jet is the jet propagation length defined by the location of the bow shock in the z-axis (see Figure 1).So the background medium is set to be uniform with the parameters relevant to the galactic environment, ρ 0 = 3.0 × 10 −25 g cm −3 and P 0 = 2.4 × 10 −10 dyne cm −2 , corresponding to T = 6.0 × 10 6 K (e.g., Perucho et al. 2014).In the models with r j = 31.6− 100 pc, the socalled King profile is adopted for the background density, ρ b (r) = ρ 0 (1 + (r/r core ) 2 ) −3β K /2 , with r core = 12 r j and β K = 0.73 (King 1962).The core radius, r core , for each model is listed in the fifth column of Table 1.External gravity is imposed to balance the pressure gradient due to the stratification.The background setup may be somewhat arbitrary, but the results presented in the next section are not sensitive to it.The jet-tobackground density ratio is η = ρ j /ρ 0 = 10 −5 , and the pressure radio is P j /P 0 = 1; then, P j /ρ j c 2 = 9.3 × 10 −2 , so initially the injected jet flow is not thermally relativistic.On the other hand, with Γ j ∼ 4 − 35, the jet is kinetically driven.
As in Paper I, simulations are performed using the RHD code described in Seo et al. (2021b), which is based on the WENO and SSPRK schemes.To correctly produce the thermodynamics of relativistic fluids, the RC version of the equation of state is employed (Ryu et al. 2006).In all the simulations, the same grid  1 for the jet parameters.Note that the length scales of the images are different in the three panels.In Panel (d), uz along the z-axis is shown for the three jet models, as a function of z/ljet, where ljet is the jet propagation length.In Q46-r1000 (red) with ⟨Γ⟩spine ∼ 10.3 and Q44-r32 (green) with ⟨Γ⟩spine ∼ 5.8, the jet-spine flows remain relatively unimpeded.On the other hand, in Q44-r100 (blue) with ⟨Γ⟩spine ∼ 2.7, the jet-spine flow suffers more significant deceleration.The mean Lorentz factor, ⟨Γ⟩spine, is estimated by averaging Γ along the jet axis in the region of z ≤ 2/3 ljet.resolution of r j /∆x = 5 is used, where ∆x is the size of grid zones.The simulation results are presented in terms of the jet crossing time, t cross = r j /u head , where ) is the approximate advance speed of the jet head and η r = (ρ j h j Γ 2 j )/(ρ 0 h 0 ) is the relativistic density contrast (Martí et al. 1997).Here, h 0 is the enthalpy of the background medium.Simulations run up to t end = 50 − 75 t cross , and then the jet propagates typically up to l jet ∼ 100 r j .To break the rotational symmetry of the system, a small precession of period 10 t cross and angle 0.5 o is added.We note that this is not necessarily related to the real physical precession of the inner jet from the central engine.

Flow structures of FR-I jets
Figure 1 compares the overall jet-flow dynamics and the deceleration of the jet-spine in the Q46-r1000, Q44-r32, and Q44-r100 models.Among the three, the FR-II model, Q46-r1000 with Γ j = 22.6, suffers the least deceleration, and the jet-spine flow remains relativistic almost up to the jet head.The two Q44 models have the same jet power Q j but different injection Lorentz factor Γ j .In Q44-r32 with Γ j = 11.2, the jet-spine remains relativistic for most of its path, resulting in ⟨Γ⟩ spine ∼ 5.8.The mean Lorentz factor, ⟨Γ⟩ spine , is estimated by averaging Γ along the jet axis in the region of z ≤ 2/3 l jet where the jet-spine is well preserved.The values of ⟨Γ⟩ spine at t end in our jet models are listed in the last column of Table 1.In the Q44-r100 model with Γ j = 3.9, the jet flow suffers more significant de- celeration, leading to ⟨Γ⟩ spine ∼ 2.7.As a result, this jet has a broader cocoon with a larger ratio of width to length.This comparison demonstrates that, in addition to the jet power Q j , the Lorentz factor Γ j plays a crucial role in determining the dynamics of jet flow, in particular the deceleration of the jet spine.
Figure 2 shows the two-dimensional (2D) distributions of nonlinear structures in Q44-r100, which are expected to play key rolls in CR acceleration.These include the shock speed (u s ), the turbulent velocity (u turb ), and the velocity shear (Ω shear ).This FR-I jet initially exhibits relativistic speeds, producing strong shear flow only near the jet nozzle (z ≲ 3 kpc).But it soon transitions to more or less turbulent flow for z ≳ 3 kpc, which subsequently generates a diffusive jet head and a broad cocoon.Shocks induced in the jet-spine and the backflow (cocoon) are subrelativistic, with u s ≲ 0.2 − 0.3c.
The Q46-r1000 jet, on the other hand, induces relativistic shocks with u s ∼ 0.4-0.5calong the jet-spine flow all the way to the termination region, as shown in Figure 4 of Paper I. Compared to the Q44-r100 model, there is also stronger shear along the interface between the jet-spine and the backflow, which causes more turbulence through the Kelvin-Helmholtz instability.Hence, CRs are expected to be energized to higher energies in higher-power FR-II jets like the Q46-r1000 jet, since DSA and RSA are the dominant processes for CR acceleration in jet-induced flows (see Section 4.1 below).

MONTE CARLO SIMULATIONS FOR CR TRANSPORT
The acceleration of UHECRs in FR-I jets is evaluated using Monte Carlo (MC) simulations that follow the transport and energization of particles in simulated jet flows.The interactions of particles with underlying MHD fluctuations are fundamental in determining the acceleration of CRs as well as their confinement and escape from the system.Thus, the distribution and strength of magnetic fields and the particle scattering laws are the key elements of our MC simulations.Here, we brief those, leaving the details to Paper I. A crucial measure of particle scattering and transport is the gyroradius, which depends on the local magnetic field strength B: where Z i and E are the charge and the energy of CR nuclei.

Magnetic Field Model
Since our RHD simulations are hydrodynamic, we adopt physically motivated model for the distribution of magnetic field.Our recipes incorporate magnetic field amplification via small-scale turbulent dynamo, E B,turb ≡ B 2 turb /8π ≈ E turb (e.g., Cho et al. 2009), and CR streaming instabilities at shocks, P B,Bell ≡ B 2 Bell /8π ≈ (3/2)(u s /c)P CR (e.g., Bell 2004).Here, E turb ≈ Γ turb (Γ turb − 1) ρc 2 with Γ turb = 1/ 1 − (u turb /c) 2 is the kinetic energy density of turbulent flow; u s is the shock speed, and P CR is the CR pressure at shocks approximated as ∼ 0.1ρ 1 u 2 s with the pre-shock density ρ 1 .In the regions of shock-free and weak turbulence, P B,P ≡ B 2 P /8π ≈ P/β p is adopted with the plasma beta β p ≈ 100.We get u turb , u s , ρ 1 , and P from the simulated jet flows (see Seo et al. Figure 3 (a) and (b) display the ratio of the magnetic energy density to the kinetic energy density, E B /E kin , and the ratio of the magnetic pressure to the gas pressure, P B /P , respectively.Here, E kin = Γ f (Γ f − 1)ρc 2 , and Γ f = 1/ 1 − (u/c) 2 is the fluid Lorentz factor calculated with the fluid speed u.The figures show that the magnetic field is subdominant, that is, either E B < E kin or P B < P in most regions.In our magnetic field model, B is set mostly by B turb in the jet spine as well as in the dynamically active parts of the cocoon, such as the interfaces between the jet-spine and the cocoon and between the cocoon and the shocked background medium.On the other hand, in relatively quiet regions of the cocoon, B is fixed mostly by B P .Figure 3 (c) plots the probability distribution functions (PDFs) of the magnetic field strength in the jet-spine flow (red lines) and the cocoon (blue lines) for three models, Q46-r1000 (FR-II), Q44-r100, and Q44-r32.The peaks represent quiet regions, while the broad distributions mainly include dynamically active parts.The magnetic field is stronger if E turb and P are larger.
The magnetic field strength in the observer frame can be approximated as

Particle Scattering Laws
CR particles are assumed to scatter elastically off underlying fluctuating magnetic fields on all relevant scales that are frozen in the local flow.The scattering is modeled with Bohm diffusion (e.g., Bohm 1949), Kolmogorov-type resonant scattering (e.g., Stawarz & Petrosian 2008), or nonresonant scattering (e.g., Sironi et al. 2013).In general, the scattering mean free path (MFP) can be formulated as where L 0 is the coherence length scale of fluctuating magnetic fields.The local Hillas energy is derived from the confinement condition, r g ≈ L 0 , where B(x, y, z) is the magnetic field strength in the local fluid frame (i.e., scattering frame).We set the charge Z i = 1 since only proton is considered in this work.We approximate L 0 ≈ r j , as the coherence length scale of induced turbulence would be comparable to the jet radius (Seo et al. 2021a).
As the fiducial model, we adopt δ = 1/3 (Kolmogorov) for E ≤ E H,L0 and δ = 1 (Bohm) for E > E H,L0 , but δ = 1 for all energy around shocks.We also consider δ = 1/3 and 2 (nonresonant) for E > E H,L0 as comparison cases.In addition, as another comparison case, we consider the diffusive transport model in turbulent magnetic fields presented by Harari et al. (2014): where the coefficients are set as a H = 4, a I = 0.9, and a L = 0.23 for the Kolmogorov spectrum.Here, the first term ∝ E 2 represents the nonresonant scattering at high energies in the quasirectilinear propagation regime, while the third term ∝ E 1/3 corresponds to the resonant scattering at low energies in the diffusive scattering regime.The middle term ∝ E is introduced to describe a smooth transition from low to high energies.This MFP model is similar to our model with δ = 2 for E > E H,L0 , except the smooth transition around E H,L0 .
After scattering, "restricted" random walks are applied; high-energy particles with λ f (E) ≳ L 0 move in mostly forward directions, while particles with λ f (E) ≲ L 0 scatter almost isotropically (see Paper I for details).We also consider fully isotropic scattering as a comparison case.See Section 4.5 for the differences due to different particle scattering models.

Injection of Galactic CRs
Seed CRs are injected into the jet inflow with an initial spectrum of dN inj /dE inj ∝ E −2.7 inj for E inj = 10 13 − 10 15 eV.They represent a population of galactic CRs accelerated within the host galaxy.While this injection model including the slope and range is somewhat arbitrary, we found that the spectrum of escaping UHECRs doesn't depend on it at all.This is because, in our MC simulations, the particle energy increases by a factor of ∼ 10 5 ; once we inject a sufficient number of galactic CRs, we obtain a smooth spectrum of escaping UHECRs up to ∼ 10 21 eV.Hence, in effect, the injected spectrum behaves almost like a delta function in energy.Since CRs are treated as test particles, the amplitude of the resulting energy spectrum scales linearly with the number of injected CRs.The trajectories of CRs are followed according to the prescribed models for scattering MFP and random walks.Figure 4 illustrates the trajectory and energization history of a sample particle.Upon each scattering, a net energy change, ∆E, arises as a result of the Lorentz transformation between the moving fluid frame and the simulation (laboratory) frame.In this example, the particle gains energy mainly through multiple episodes of DSA and GSA up to ∼ 1 EeV, advecting along the jet-spine flow.Then, it undergoes diffusive scatterings in the cocoon (between the triangle and star symbols in Panel (a)) and eventually escapes from the turbulent head region.As shown in the shaded box in Panel (b), the acceleration is rather inefficient for 0.01 ≲ t − t inj ≲ 0.06, during which the particle ex-

Relative Importance of Acceleration Processes
Similar to FR-II jets, various CR acceleration processes, including DSA, TSA, GSA, and nGSA, are at play within FR-I jets.We first utilize the acceleration timescales, t DSA , t TSA , t GSA , and t nGSA (see Eqs. [11], [14], [15], and [16] in Paper I) to compare the relative efficiencies of these acceleration processes (APs).For illustrative purposes, these timescales are displayed in Figure 5(a) for the Q44-r100 model with the following representative parameters: the shock compression ratio, χ ∼ 4, u s /c ∼ 0.1, B ∼ 30 µG, the turbulence speed, |u turb |/c ∼ 0.1, L 0 ∼ 0.1 kpc, the velocity shear, Ω shear = ∂u z /∂r ∼ 0.05 c/r j with the z-velocity along the jet direction, u z , Γ z = 1/ 1 − u 2 z /c 2 ∼ 1, and δ = 1 for MFP.The figure indicates that seed CRs would gain energies at first via DSA, and then RSA (GSA and nGSA) becomes increasingly significant above E ∼ 1 EeV.The particles with λ f (E) ≲ r j stay within the shear layer spread over the jet-spine and the backflow, and gain energy mainly via GSA.In contrast, those with λ f ≳ r j can be scattered across the jet-backflow interface and boosted to higher energies via nGSA.
The relative contributions of APs in our jet models are estimated as follows.For each scattering event, a fraction of ∆E is distributed among different APs with weight ξ AP = t −1 AP / AP t −1 AP , where the summation includes DSA (shock), TSA (turbulence), and GSA (shear); nGSA is not considered separately since only a small fraction of high-energy particles undergo it.For all particles escaping from the system up to t end , whose final energy lies in the logarithmic bin of [log E, log E + d log E], ξ AP ∆E is summed to make the cumulative energy gain, E AP (E), and E tot = E shock + E turb + E shear .Figures 5(b)-(c) show E AP (E)/E tot for the two Q44 models presented in Figure 1.The figures confirm that whereas DSA is dominant for E ≲ 1 EeV, shear acceleration becomes important beyond E ≳ 1 EeV.TSA plays a secondary role.

Energy spectrum of CRs from FR-I Jets
The energy spectrum of the CRs that are accelerated and escape from the jet system is typically timedependent, transitioning from an "age-limited" regime (bluish lines) to a "size-limited" regime (reddish lines), as shown in Figure 6(a).During the early stage, all particles are accelerated mainly by encountering multiple shocks and gradual velocity shear in the jet-spine flow (see Figures 4), and the spectrum extends to higher energies over time as the maximum energy of CRs increases with the jet's age.In the late stage, particles are confined in the jet system and continuously energized until they escape.The shape of the time-integrated, cumulative energy spectrum of all particles escaping from the jet converges to its time-asymptotic form.
Figures 1 and 2 show that the jet spine typically disperses over a width of ∼ several r j , and the corrugated interface between the upward-moving jet flow and the downward-moving backflow spreads over of order ∼ r j .Only UHECRs energized to E ≳ 1 EeV have large enough MFPs to cross the shear layer between the jetspine and the backflow, as mentioned above.Hence, as CRs are initially energized mainly via DSA, GSA first becomes important, and then UHECRs with E ≳ 1 EeV cross and recross the shear interface and are energized via nGSA; each cycle of cross-recross increases the energy by a factor of (Γ 2 ∆ − 1) (e.g., Rieger 2019; Caprioli 2015; Mbarek & Caprioli 2019).Here, Γ ∆ is the Lorentz factor of the velocity jump across the shear layer.
In Paper I, we categorized UHECRs produced at RGs into three types, depending on whether they experience an additional boost via nGSA before escape and how they escape from the jet structure (see Figure 9 in Paper I).In CASE 1, particles are energized mainly via DSA and RSA, and then escape from the cocoon without undergoing an additional boost via nGSA.In CASE 2, after gaining energies via DSA and RSA, particles undergo an additional boost and then escape from the cocoon.CASE 3 is the same as CASE 2, except that UHECRs escape directly from the jet-spine into the background medium.In Figure 6(b), the time-asymptotic, cumulative energy spectra, N (E), of the particles categorized into different CASEs are separately plotted for the Q44-r32 model at t end .Most particles belong to the CASE 1 category, while only a small fraction becomes CASE 2 and 3 by undergoing an additional boost to higher energies via nGSA before escape.
Figure 7 shows N (E) of all escaping particles for different FR-I jet models at t end with color-coded solid lines.The spectra shift to higher energies for models with higher Q j , while the models with the same Q j but larger Γ j have spectra that extend to higher energies.In contrast to the double-power law in Equation ( 1) for FR-II jets, the fitting of the CR energy spectra for FR-I jets requires an exponential cutoff: Table 1 lists the four fitting parameters, s 1 , s 2 , E break , and Γ fit , with which Equation ( 6) could approximate the energy spectra of escaping CRs obtained from our MC simulations; the fittings are shown with dot-dashed lines in Figure 7.

Spectral Slopes and Exponential Cutoff
The slopes, s 1 and s 2 , in Table 1 exhibit relatively weak dependence on jet models; they are, on average, s 1 ≈ −0.6 below the break energy, E break , and s 2 ≈ −2.6 above the break energy.These values for FR-I jets are almost identical to those for FR-II jets presented in Paper I, since the two types of RG jets share similar flow dynamics and hence the same acceleration processes operate.
The power law for E < E break , ∼ E −0.6 , is fixed mostly by CASE 1 particles, as shown in Figure 6(b); we interpret it as a consequence of multiple episodes of DSA and RSA in the jet-spine flow.Multiple shock passages are known to lead to a E −1 energy spectrum for a large number of shocks, independent of the shock compression ratio (e.g., Melrose & Pope 1993;Kang 2021).Moreover, MC simulation studies by Ostrowski (1998) and Kimura et al. (2018) showed that the energy spectrum of particles accelerated via nGSA scales as E −1 − E 0 .They used a simpler, cylinder-shaped jet-cocoon system with discrete velocity shear as the background flow for their MC simulations, instead of time-evolving, complex jetinduced flows with lots of shocks and turbulence, as seen in Figure 2. The slope −0.6 is close to that of multiple  6) are overlaid with dot-dashed lines.See Table 1 for the jet and fitting parameters.Note that the abscissa in Panel (c) stretches to a higher energy than in other panels.
DSA and also in the range of the values of Ostrowski (1998) and Kimura et al. (2018).Our energy spectra are harder than the canonical single DSA spectrum of The behavior of the spectra for E > E break , on the other hand, is the combined result of particles of all CASEs (Figure 6(b)).In the conventional Hillas picture in which CRs are confined by an isotropic system with mean magnetic field strength ⟨B⟩ and size L, the spectrum of escaping particles would decrease exponentially above the Hillas energy (Equation (4) with ⟨B⟩ and L).In the case of jet-induced flows under consideration, the escape scenario is modified by the following two factors: particles are confined in the elongated cocoon, and some UHECRs receive the ⟨Γ⟩ 2 spine boost via nGSA and reach above the Hillas energy before escaping from the jet.This picture is consistent with the so-called espresso acceleration model proposed by Caprioli (2015) and Mbarek & Caprioli (2019).As we mentioned in Paper I, the particle escape from an idealized, infinite cylindrical volume leads to an energy spectrum that scales as E −2 .Somewhat steeper spectra with ∼ E −2.6 are obtained in our simulations since our jets have finite-sized, barrel-shaped cocoons.
At the highest energy end, an exponential cutoff would appear at E break ⟨Γ⟩ 2 spine , as a consequence of the nGSA boost.Indeed, the values of Γ fit in the exponential cutoff of Equation ( 6) are close to those of ⟨Γ⟩ spine , as can be seen in Table 1.Hence, it is ⟨Γ⟩ spine that primarily governs the highest energy end of the spectrum.When ⟨Γ⟩ spine ≳ 10, as in the case of FR-II jets, the cutoff energy is sufficiently high, so the highest energy part of the energy spectrum would resemble a power law.In contrast, with ⟨Γ⟩ spine ≲ 6 in our simulations, the energy spectrum for FR-I jets needs an exponential cutoff.

Break Energy of Power-Law
The break energy, E break , in the model spectrum is mainly governed by the CASE 1 particles that are confined and escape from the cocoon (Figure 6(b)).If the acceleration is sufficiently efficient, it would appear at the Hillas energy, estimated with the mean magnetic field strength, ⟨B⟩, and the transverse size of the cocoon, W. Note that E H is different from the local Hillas Energy, E H,L0 , in Equation ( 4), which is defined by the local magnetic field strength, B(x, y, z), and the coherent scale of turbulence, L 0 .On the contrary, if the acceleration is slow, CRs may escape diffusively from the turbulent cocoon before they reach the Hillas energy, E H ; then the maximum energy is limited by escape via fast spatial diffusion.The sample trajectory shown in Figure 4 is for a CASE 1 particle that undergoes inefficient acceleration in the diffusive head region.The timescale for particles to diffuse across the cocoon is t diff ≈ (W/2) 2 /(cλ f ).The particles around E break have λ f ≳ r j and hence could cross the jet-backflow interface.Then, nGSA would be operative, and its acceleration timescale can be approximated as t nGSA ≈ λ f /c⟨Γβ⟩ 2 acc (e.g.Kimura et al. 2018), where β ≡ u z /c and ⟨Γβ⟩ acc is the volume-averaged value over the region of nGSA. Figure 5(a) shows that t D is smaller than t nGSA for E ≳ E break in Q44-r100, in which the jet flow is significantly decelerated.Thus, in the case of slow acceleration, the maximum energy may be estimated from the condition t D ≈ t nGSA .For the fiducial model with λ f (p) ∝ E, this condition results in The twelfth column of Table 1 lists the values of E H , which are calculated using the volume-averaged magnetic field strength in the cocoon and the width of the cocoon at z = l jet /2 in simulated jets at t end .The thirteenth column lists the values of E D at t end using ⟨Γβ⟩ acc , the volume-averaged value at the jet-spine in the upper half of the jet (z = (1/2 − 1)l jet ) where the majority of CASE 1 particles are confined before they escape to the background (see Figure 4).Both E H and E D approach to time-asymptotic values in our jet models, although their time evolution is not explicitly presented.In the FR-I models with large injection Lorentz factors, Γ j = 11.1 and 34.5 (⟨Γ⟩ spine ≳ 6), E break ≈ E H , indicating that RSA is efficient and the majority of particles reach E H before they escape.In contrast, in the models with Γ j = 3.9 (⟨Γ⟩ spine ≲ 3), ⟨Γβ⟩ acc ∼ 0.4 − 0.5, and hence E break ≈ E D ≲ E H , implying that RSA may not be sufficiently efficient because the jet flow is only mildly relativistic.Therefore, E break ≈ E H may be adopted in our model spectrum for jets with ⟨Γ⟩ spine ≳ 6, whereas for jets with ⟨Γ⟩ spine ≲ 3, E break ≈ E D needs to be employed.
The break energy, E break , is expected to depend on the properties of the jet and hence the jet model parameters such as Q j , r j , and the background density ρ b , in our simulations.In Paper I, we found that the time-asymptotic value of E break ≈ E H scales as Q 1/3 j for FR-II jets.
To elicit the approximate dependence of E H on Q j as well as r j and ρ b , we here consider a relativistic jet injected into a uniform background medium, and assume an idealized cylindrical volume with radius R c and vertical length l c for the jet-induced cocoon.If a significant fraction of the kinetic energy of the jet dissipates into the internal energy of the cocoon, the gas pressure of the cocoon can be approximated as P c ≈ (Q j t)/(πR 2 c l c ).We take R c ≈ u r t and l c ≈ u head t; the radial expansion speed is given as u r ≈ (P c /ρ b ) 1/2 from the ram pressure balance, and the head advance speed is given as u head ≈ √ η r c with u j ≈ c and η r ≈ Q j /(πr 2 j c 3 ρ b ) ≪ 1 (see Section 2.1).Using these relations, we obtain the scaling for the cocoon pressure: Then, the width of the cocoon can be estimated as Moreover, in our prescription for B, the magnetic field is stronger if the turbulence energy is larger and the gas pressure is larger (see Section 3.1); hence, we may approximate the mean magnetic field in the cocoon as Then, E H scales as We note that while the width increases in time as W ∝ t 1/2 , the magnetic field strength decreases as ⟨B⟩ ∝ t −1/2 .As a consequence, E H is independent of time.As a matter of fact, this is one of the reasons why the shape of the energy spectra in our simulations approaches a time-asymptotic form, as mentioned above.
In Table 1, the Hillas energy, E H , scales approximately as Q 1/4 j , rather than as Q 1/3 j .This could be because the FR-I jets considered here are more compact with smaller spatial dimensions than the FR-II jets of Paper I. On the other hand, the dependence on r j is not obvious, possibly because the jets propagate through the background of different stratifications (that is, different r core ) in models with different r j .
Observed RG jets are commonly quantified with radio or X-ray luminosities that may be related to Q j , but other properties such as r j and ρ b are not easily estimated from observations.Hence, we propose an expression for E break for the model spectrum of UHECRs from FR-I jets as where Q 0 = 3.5 × 10 44 ergs −1 and α = 1/4.The first fuzzy factor, φ ∼ min[1, ⟨Γβ⟩ acc ], becomes φ ≈ 1 for high-power jets with ⟨Γ⟩ spine ≳ 6, while φ ∼ ⟨Γβ⟩ acc ∼ 0.4 − 1 for low-power jets with ⟨Γ⟩ spine ≲ 6, based on the jet simulations considered here.This factor could be roughly estimated if the jet-spine structure is resolved in observations.The second fuzzy factor, ξ ≈ (r j /r 0 ) 1/2 (ρ b /ρ 0 ) 1/4 , where r 0 = 100 pc, ρ 0 = 3 × 10 −25 gcm −3 .Although it would be difficult to constrain ξ from observations of RG jets, ξ ∼ 0.4 − 1 for the jet models considered here.The above expression may be applied to FR-II jets as well, with α = 1/3 and φ ∼ 1.

Model Dependence of Energy Spectrum
The spectrum in Equation ( 6) is the result of MC simulations based on a number of modelings described in Section 3. Here, we examine the dependence on particle scattering and magnetic field models, as those could affect the resulting spectrum.
In Figure 8(a), the energy spectra with δ = 1/3 (yellow dotted) and 2 (light green dotted) for E > E H,L0 are compared to the fiducial spectrum with δ = 1 (black solid) for the Q44-r100 model.Below E H,L0 ∼ 1 EeV, the three spectra are basically identical.Above E H,L0 , the spectrum with δ = 1/3 extends to higher energies, while the spectrum with δ = 2 falls off at lower energies.This is because with δ = 1/3, high-energy particles suffer more scatterings compared to the Bohm model with δ = 1, and hence reach higher energies.In contrast, with δ = 2, it works in the opposite way.The cyan dot-dashed line shows the spectrum for the diffusion model in Equation (5) proposed by Harari et al. (2014), which almost overlaps the spectrum with δ = 2, as expected.In the same panel, the energy spectrum of the fully isotropic scattering model is also shown with a purple dashed line.With isotropic scatterings rather than restricted random walks, the chance for high-energy particles to cross the shear interface increases, and nGSA is enhanced.Hence, the spectrum extends to higher energies.
Figure 8(b) shows how the magnetic field model affects the resulting spectrum.The black solid line draws the spectrum of the fiducial model with E B,turb = E turb and P B,P = 0.01P ; the cyan dotted line plots the spectrum with E B,turb = E turb and P B,P = 0.1P , and the purple dashed line plots the spectrum with E B,turb = (1/2)E turb and P B,P = 0.01P .This comparison demonstrates that the spectrum of escaping UHECRs would shift to higher energies if the magnetic field is stronger (cyan dotted), while it works in the opposite way if the magnetic field is weaker (purple dashed).
Figure 8 implies that the predicted energy spectrum generated in RG jets is, to a certain extent, contingent upon the modeling details of physical processes, including the characteristics of magnetic turbulence and particle scattering.

SUMMARY
In paper I, we performed RHD simulations to reproduce the dynamics of relativistic jets in FR-II RGs and MC simulations to track the transport and energization of CR particles in the simulated jet-induced flows.Following the same numerical approach of RHD and MC simulations, we here investigate the jet flows and the acceleration of UHECRs in FR-I RGs.
In FR-I jets with smaller power Q j and smaller Lorentz factor Γ j , the jet-spine gets decelerated more significantly.Compared to high-power FR-II jets, these low-power jets develop broader cocoons and more diffusive head regions, containing more numerous shocks and turbulence but weaker velocity shear.
The main results related to particle acceleration are summarized as follows: 1. DSA at shocks, TSA by turbulence, and RSA (both GSA and nGSA) at shear operate for the acceleration of UHECRs.While DSA is dominant in energizing CRs up to E ∼ 1 EeV, RSA becomes increasingly important beyond E ∼ 1 EeV.TSA plays a secondary role.
2. The time-integrated, cumulative energy spectrum of all particles escaping from the jet system approaches a time-asymptotic form.It can be modeled by a double power law with an exponential cutoff as given in Equation (6).
3. The break energy of power laws, E break , is either fixed by the Hillas confinement condition or limited by particle escape from the cocoon via fast spatial diffusion.It depends on jet parameters and may be modeled as Equation ( 13).
4. The power-law slopes are on average s 1 ≈ −0.6 below E break and s 2 ≈ −2.6 above E break , and show only a weak dependence on jet models.The hard spectrum, E −0.6 , for E < E break is thought to result from multiple episodes of DSA and RSA in the jet system.In contrast, the steep spectrum, E −2.6 , for E > E break is caused by the particle confinement and escape from the elongated cocoon.
5. The exponential cutoff energy is given as ∼ E break ⟨Γ⟩ 2 spine , where ⟨Γ⟩ spine is the mean Lorentz factor of the jet-spine.It is governed by nGSA at the shear interface, which increases the CR energy by a factor of ∼ ⟨Γ⟩ 2 spine .The model spectrum derived for FR-I jets in this work, along with the one for FR-II jets in Paper I, may be utilized to investigate the contributions of nearby RGs as well as cosmological populations of RGs to UHECRs observed at Earth.This would help elucidate the origins and mechanisms behind the production of UHECRs.We leave such a study as a future work.

Figure 1 .
Figure1.2D slice images of the density, log ρ, and the z-velocity along the jet direction, uz, for three models, (a) Q46-r1000 (FR-II), (b) Q44-r32, and (c) Q44-r100, at t = t end .See Table1for the jet parameters.Note that the length scales of the images are different in the three panels.In Panel (d), uz along the z-axis is shown for the three jet models, as a function of z/ljet, where ljet is the jet propagation length.In Q46-r1000 (red) with ⟨Γ⟩spine ∼ 10.3 and Q44-r32 (green) with ⟨Γ⟩spine ∼ 5.8, the jet-spine flows remain relatively unimpeded.On the other hand, in Q44-r100 (blue) with ⟨Γ⟩spine ∼ 2.7, the jet-spine flow suffers more significant deceleration.The mean Lorentz factor, ⟨Γ⟩spine, is estimated by averaging Γ along the jet axis in the region of z ≤ 2/3 ljet.

Figure 2 .
Figure 2. 2D slice images of the quantities that exhibit nonlinear flow dynamics in the Q44-r100 model at t = t end : (a) shock speed, us, (b) turbulent flow velocity, u turb , and (c) velocity shear, Ω shear .

Figure 3 .
Figure 3. Panel (a): Ratio of the magnetic energy density to the kinetic energy density, EB/E kin , for the Q44-r100 model, where E kin = Γ f (Γ f − 1)ρc 2 with the fluid Lorentz factor Γ f .Panel (b): Ratio of the magnetic pressure to the gas pressure, PB/P , for the same model.Here, the "comoving" magnetic field strength B = max(Bp, B turb , B Bell ) in the fluid frame is used.The two ratios are shown for the region inside the bow shock, whereas the background medium (BM) is colored white.Panel (c): Probability distribution functions of B in the jet-spine flow (red) and the backflow (blue) for the FR-II (dotted lines), Q44-r100 (solid lines), and Q44-r32 (dashed lines) models.All are shown at t end .2021a, for the evaluation of u turb ).The highest estimate is chosen among the three model values, B(x, y, z) = max(B turb , B Bell , B P ), as the local comoving magnetic field strength.Figure3(a) and (b) display the ratio of the magnetic energy density to the kinetic energy density, E B /E kin , and the ratio of the magnetic pressure to the gas pressure, P B /P , respectively.Here, E kin = Γ f (Γ f − 1)ρc 2 , and Γ f = 1/ 1 − (u/c) 2 is the fluid Lorentz factor calculated with the fluid speed u.The figures show that the magnetic field is subdominant, that is, either E B < E kin or P B < P in most regions.In our magnetic field model, B is set mostly by B turb in the jet spine as well as in the dynamically active parts of the cocoon, such as the interfaces between the jet-spine and the cocoon and between the cocoon and the shocked background medium.On the other hand, in relatively quiet regions of the cocoon, B is fixed mostly by B P .Figure3 (c) plots the probability distribution functions (PDFs) of the magnetic field strength in the jet-spine flow (red lines) and the cocoon (blue lines) for three models, Q46-r1000 (FR-II), Q44-r100, and Q44-r32.The peaks represent quiet regions, while the broad distributions mainly include dynamically active parts.The magnetic field is stronger if E turb and P are larger.The magnetic field strength in the observer frame can be approximated as B obs ≈ Γ f • B(x, y, z); then typical values range B obs ∼ 10 − 30µG in the backflow and B obs ∼ 10 − 100µG in the jet-spine flow in our model.These estimates lie within the range of values from Xray and radio observations (e.g.,Begelman et al. 1984;Kataoka & Stawarz 2005;Anderson et al. 2022).The y, z); then typical values range B obs ∼ 10 − 30µG in the backflow and B obs ∼ 10 − 100µG in the jet-spine flow in our model.These estimates lie within the range of values from Xray and radio observations (e.g., Begelman et al. 1984; Kataoka & Stawarz 2005; Anderson et al. 2022).The general features of the spatial distribution of B obs for a FR-II jet are shown in Figure 3(b) of Paper I. These features also apply to lower-power FR-I jets.Our magnetic field model employs E B,turb /E turb ≈ 1 and β p ≈ 100.While they result in B obs consistent with observations, these values are somewhat arbitrary.Hence, as comparison models, we consider the cases of E B,turb /E turb ≈ 1/2 and β p ≈ 10; the former gives rise to a weaker B than the fiducial model, while the latter gives a stronger B. With these, we examine the effects of magnetic field modeling on the acceleration of UHECRs in Section 4.5.

Figure 4 .
Figure 4. Panel (a): Trajectory of a sample particle since the injection into the jet flow, superimposed on the distribution of log ρ in the Q44-r100 model.The trajectory is color-coded by the energy of the particle.Panel (b): Energization history of the same particle along the trajectory.The shaded region indicates the period during which the particle scatters within the diffusive head region.The starting and ending locations of that period are marked by the white triangle and black star symbols, respectively, in panel (a).

Figure 5 .
Figure 5. Panel (a): Timescales of the acceleration processes (APs) that energize CRs in jet-induced flows as a function of the energy of particles for the parameters relevant to the Q44-r100 model.nGSA is operative only when λ f is large enough to cross the jet-backflow interface, and hence tnGSA (blue dashed line) is drawn for E > 0.1 EeV.The diffusion timescale, tD, (black dotted line) is also shown for comparison.Panels (b-c): Fractions of the cumulative energy gains due to different APs as a function of the energy of particles escaping from the jet system for two models, (b) Q44-r32 and (c) Q44-r100.See the text for details.periences diffusive scatterings within the upper cocoon (z ≳ 3 kpc).

Figure 6 .
Figure 6.Panel (a): Time-integrated energy spectrum, E[dN /dE/Ntot], of all particles escaping from the system up to a given time, t, in the Q44-r32 model.Here, Ntot = (dN /dE)dE.The lines are color-coded from blue to brown, based on the time, t/tcross.Panel (b): Time-asymptotic, cumulative energy spectra, drawn separately for particles classified as CASE 1, CASE 2, and CASE 3, in Q44-r32.See the main text for the description of CASEs.The back solid line shows the total spectrum, and the black dotted line is the fitting to it.

Figure 7 .
Figure7.Time-asymptotic, cumulative energy spectra of all particles escaping from the jet system up to t end , E[dN /dE/Ntot], (a) for the Q42 and Q43 models, (b) for the Q44 models, and (c) for the Q45 models.Each spectrum is normalized with Ntot = (dN /dE)dE.The results of MC simulations are shown with color-coded solid lines, while the fittings to Equation (6) are overlaid with dot-dashed lines.See Table1for the jet and fitting parameters.Note that the abscissa in Panel (c) stretches to a higher energy than in other panels.

Figure 8 .
Figure 8. Panel (a): Energy spectra with different modelings for the scattering of CR particles in the Q44-r100 model.The diffusive scattering model by Harari et al. (2014) is given in Equation (5).The purple dashed line shows the case where fully isotropic scattering is adopted.Panel (b): Energy spectra with different modelings for the magnetic field distribution in the Q44-r100 model.Here, the spectrum of the fiducial case (black solid line) is the same as that in Figure 7(b).

Table 1 .
Model Parameters of Simulated Jets a and Fitting Parameters for UHECR Energy Spectrum b Q46-r1000 is a FR-II jet model, presented as Q46-η5-H in Paper I.
b Columns 8 − 11: fitting parameters for the energy spectrum of escaping CRs given in Equation (6).cIn the model name, the numbers after Q and r are the exponent of Qj and rj in units of pc, respectively.dThe Hillas energy, EH , the diffusion-limited maximum energy, ED, and the mean Lorentz factor in the jet spine, ⟨Γ⟩spine, are the values at t end .e