Acceleration and Spectral Redistribution of Cosmic Rays in Radio-jet Shear Flows

A steady-state, semi-analytical model of energetic particle acceleration in radio-jet shear flows due to cosmic-ray viscosity obtained by Webb et al. is generalized to take into account more general cosmic-ray boundary spectra. This involves solving a mixed Dirichlet–Von Neumann boundary value problem at the edge of the jet. The energetic particle distribution function f 0(r, p) at cylindrical radius r from the jet axis (assumed to lie along the z-axis) is given by convolving the particle momentum spectrum f0(∞,p′) with the Green’s function G(r,p;p′) , which describes the monoenergetic spectrum solution in which f0→δ(p−p′) as r → ∞ . Previous work by Webb et al. studied only the Green’s function solution for G(r,p;p′) . In this paper, we explore for the first time, solutions for more general and realistic forms for f0(∞,p′) . The flow velocity u = u(r) e z is along the axis of the jet (the z-axis). u is independent of z, and u(r) is a monotonic decreasing function of r. The scattering time τ(r,p)=τ0(p/p0)α in the shear flow region 0 < r < r 2, and τ(r,p)=τ0(p/p0)α(r/r2)s , where s > 0 in the region r > r 2 is outside the jet. Other original aspects of the analysis are (i) the use of cosmic ray flow lines in (r, p) space to clarify the particle spatial transport and momentum changes and (ii) the determination of the probability distribution ψp(r,p;p′) that particles observed at (r, p) originated from r → ∞ with momentum p′ . The acceleration of ultrahigh-energy cosmic rays in active galactic nuclei jet sources is discussed. Leaky box models for electron acceleration are described.

. In this paper, we explore for the first time, solutions for more general and realistic forms for ¥ ¢ ( ) f p , 0 .The flow velocity u = u(r)e z is along the axis of the jet (the z-axis).u is independent of z, and u(r) is a monotonic decreasing function of r. , where s > 0 in the region r > r 2 is outside the jet.Other original aspects of the analysis are (i) the use of cosmic ray flow lines in (r, p) space to clarify the particle spatial transport and momentum changes and (ii) the determination of the probability distribution y ¢ ( ) r p p , ;
Active galactic nuclei (AGNs) jets are observed at a large range of scales, extending from the black hole event horizon to megaparsec distances (e.g., Blandford et al. 2019).
First-order Fermi acceleration at relativistic shocks using the pitch angle focusing transport equation was investigated by Kirk & Schneider (1987a, 1987b), Kirk & Duffy (1999), Achterberg et al. (2001), and Pelletier et al. (2017).Selfconsistent studies of particle acceleration at relativistic shocks using particle-in-cell (PIC) codes were studied by Spitkovsky (2008aSpitkovsky ( , 2008b)), Martins et al. (2009), Sironi & Spitkovsky (2009, 2011), and Sironi et al. (2013).A PIC code was used by Sironi et al. (2021) to investigate particle acceleration in relativistic shear flows.The advantage of PIC code simulations is that a consistent solution for both the particles and electromagnetic fields can be obtained by this method.However, because of computational constraints, the ratio of the ion to electron mass m i /m e imposes constraints on the computation time.This is not a constraint for electron-positron plasmas for which m p /m e = 1, where m p is the positron mass.However, it is an important constraint for an electron-proton plasma m p /m e = 1836.Sironi & Spitkovsky (2009) used m p /m e = 25 for their electron-proton plasma shock simulations.In relativistic shock simulations, the Weibel instability (WI ;Fried 1959;Weibel 1959;Medvedev & Loeb 1999; Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.Wiersma & Achterberg 2004;Achterberg & Wiersma 2007;Achterberg et al. 2007) is similar to the two-stream instability in which the two fluids consist of a cold ion beam upstream of the shock that interacts with the hot shocked downstream plasma, which can also be thought of as a beam in the shock frame.The original formulation of Weibel (1959) involved perturbation analysis of the collisionless Liouville equation, but Fried (1959) applied the analysis to beam plasmas.Perucho (2019) gives an overview of the extensive literature on linear and nonlinear instabilities of hydrodynamical and magnetized radio jets.
PIC simulations of the Kelvin-Helmholtz instability and magnetic reconnection for radio jets have been studied by Sironi et al. (2021).Mizuno et al. (2014) investigated the current-driven kink instability for relativistic jets.Nishikawa et al. (2016) investigated the evolution of global relativistic jet collimations and expansion effects on the kinetic Kelvin-Helmholtz instability (kKHI) and WI.Nishikawa et al. (2020) studied rapid particle acceleration due to recollimation shocks and turbulent magnetic fields in injected jets with helical magnetic fields.
Particle acceleration in Faranoff-Riley II (FR II) radio jets has been studied by Matthews et al. (2019), Bell et al. (2019), and Araudo et al. (2016Araudo et al. ( , 2018)).Lemoine (2019) developed a theory of generalized Fermi acceleration in which the affine connection coefficients describing non-inertial and gravitational forces are included in the Lorentz force equation.In this development, non-inertial forces also arise from the nonuniform background flow because the particle momentum is measured in the fluid frame (see also Webb 1985Webb , 1989;;Achterberg & Norman 2018a, 2018b).
Scenarios for the acceleration of UHECRs in AGN jets have included particle acceleration at sub-relativistic shocks in the jet-induced backflow in FR II jets (e.g., Matthews et al. 2019), second-order Fermi acceleration in turbulent flows in the jet cocoon (e.g., Bicknell & Melrose 1982;Hardcastle 2010), particle acceleration due to cosmic-ray viscosity in relativistic shear flows (e.g., Rieger & Duffy 2004, 2005b;Webb et al. 2018Webb et al. , 2019Webb et al. , 2020;;Rieger 2019;Wang et al. 2021Wang et al. , 2023)), discrete shear acceleration at interfaces between the jet spine and the cocoon backflow (Ostrowski 1998;Kimura et al. 2018), turbulent shear acceleration (TSA; Ohira 2013), and the socalled expresso mechanism (Caprioli et al. 2015;Mbarek & Caprioli 2019).He et al. (2023) applied the shear acceleration model of electrons in radio-jet sources of Rieger & Duffy (2019) and Wang et al. (2021), taking into account synchrotron and inverse Compton energy losses.They used the shear acceleration, leaky box model to explain the observations of X-ray spectra in largescale FR II radio sources (e.g.,3C 273,C 403,3C 17,Pictor A,3C 111,and other sources).Wang et al. (2023) carry out 3D MHD simulations of particle acceleration in relativistic jet shear flows, with sufficient resolution to investigate the formation of the spine-sheath structure and the development of turbulence for a relativistic jet propagating into a static cocoon.Different spine velocities to match both FR type I and type II jets were investigated.The simulations illustrate the growth and time development of the Kelvin-Helmholtz instability (KHI; e.g., Ferrari et al. 1978;Hardee 1979;Mizuno et al. 2007;Sironi et al. 2021) and current-driven instabilities (CDI or kink instabilities; e.g., Lyubarskii 1999;Hardee 2007;Mizuno et al. 2014) in the jet.In general, CDI dominates in high-magnetization jets with toroidal magnetic fields, and KHI dominates in lowmagnetization jets.Wang et al. (2023) find that a sheath is generated on the interface of the spine and the cocoon as a result of the KHI.The large-scale velocity in the sheath is close to linear.Seo et al. (2023) carry out relativistic hydrodynamic simulations of FR II radio jets and assess the relative importance of the different mechanisms using Monte Carlo simulations for the particle acceleration.
In this paper, we develop the model of Webb et al. (2020) for energetic particle acceleration by cosmic-ray viscosity in relativistic jets in which the jet velocity has the form of u = u(r)e z along the jet axis (the z-axis), where u(r) is a monotonically decreasing function of cylindrical radius r about the jet axis.The particles are energized by scattering back and forth across the shear flow, in the region 0 < r < r 2 , where r 2 is the outer radius of the shear flow.The particles also scatter in r > r 2 , but undergo no changes in energy in the r > r 2 region.This leads to a Dirichlet-Von Neumann boundary value problem for the energetic particle distribution function f 0 (r, p) at the edge of the jet at r = r 2 (Webb et al. 2020).The solutions are obtained for a scattering mean free time of the form t t = a ( ) p p 0 0 in the region 0 < r < r 2 , and t = t a ( ) ( ) p p r r s 0 0 2 in the outer region r > r 2 (here, s is a positive constant).
In the region r > r 2 we set b = const.The simplest choice is to set β = β 2 .ξ is known as the relativistic rapidity of the flow.More generally, it is possible to obtain other solution profiles for u(r) or β(r), but there are constraints imposed on τ(r, p) depending on u(r) in that case (see Webb et al. 2020).Rieger & Duffy (2022) have investigated the spectral slope of the accelerated particles in relativistic shear flows for a variety of fluid velocity profiles both by semi-analytical and numerical solutions of the cosmic-ray transport equation (which is an elliptic partial differential equation in r and p).In the semianalytical approach, they use a leaky box model to describe the spatial diffusive transport.They investigate convex velocity profiles for k > 1, Gaussian velocity profiles (exponentially decaying), and linear decreasing velocity profiles (for k = 1), and the effect of the velocity profile u(r) on the spectral index of shear accelerated particles.
Section 2 introduces the model of Webb et al. (2020) for the steady-state transport equation of cosmic rays in a relativistic shear flow.The boundary conditions at r = r 2 are specified in terms of a Dirichlet-Von Neumann boundary condition connecting the particle transport in the shear-free flow in r > r 2 to the inner shear flow in 0 < r < r 2 , where the particles are energized by cosmic-ray viscosity.The cosmic-ray continuity equation is written in conservation form, which leads to the concept of cosmic-ray flow lines in (r, p) space.This idea is further developed in Section 4.3.
Section 3 provides the shear acceleration solution of Webb et al. (2020) as r → ∞ .The solution for f 0 (r, p) has the form is the effective Green's function, which describes particle acceleration in 0 < r < r 2 and diffusive particle transport with no changes in energy in r > r 2 .
In Section 4, ¢ ( ) G r p p , ; is used to investigate the origin of particles in the asymptotic spectrum as r → ∞ , ¥ ¢ ( ) f p , 0 that are observed at position r with momentum p.The probability distribution y ¢ ( ) r p p , , p that particles observed at (r, p) originated at r → ∞ with momentum ¢ p is determined (see Gleeson & Webb 1975, 1980, who studied the solar modulation of galactic cosmic rays).Solutions for f 0 (r, p) for boundary spectra ¥ ¢ ( ) f p , 0 are obtained using the convolution integral (Equation ( 2)) for distributions ¥ ¢ ( ) f p , 0 consisting of a negative power law in momentum at momenta ¢ > p p b .¥ ¢ ( ) f p , 0 is continuous as p = p b and constant for < ¢ < p p 0 b .The solution for f 0 (r, p) splits into accelerated particles with > ¢ p p and decelerated particles with < ¢ p p .Section 4.3 determines cosmic-ray flow lines in (r, p) space.
In Section 5, observations and constraints on particle acceleration in radio jets are discussed.The Hillas (1984) constraints on particle acceleration, relativistic beaming in jets (Forman 1970;Lind & Blandford 1985), numerical simulations of particle acceleration in radio jets (e.g., Seo et al. 2023), and the role of affine connection coefficients to describe particle acceleration (e.g., Lemoine 2019) are discussed.Radio galaxies, microquasars and nonthermal emissions and jets, stellar mass, and supermassive black holes from mergers of stellar-mass black holes are mentioned.Particle acceleration by first-order Fermi mechanism at shocks, for example, at the head shock where the jet impinges on the interstellar medium, at helical shocks, at re-confinement shocks, the effect of stellar ablation, stellar-mass black hole jets, and shear acceleration are described.
The work of Webb et al. (2018Webb et al. ( , 2019Webb et al. ( , 2020) ) is extended in this paper to include (i) the solutions of Webb et al. (2020) for more general boundary spectra as r → ∞.The convolution Equation (2) gives the general solution, where The Hillas (1984) constraints on the width and magnetic field strength of the jet required to contain a UHECR with a given energy in the acceleration region are described.The flow lines show how monoenergetic particles with initial momentum ¢ p as r → ∞ split into accelerated and decelerated particles.Flow lines for the case of a cutoff power-law boundary spectrum (Equation ( 81)) are obtained.In the latter case, the flow lines show particle energy gains.
The main aim of the present paper is to explore the largescale implications of the acceleration of energetic charged particles in radio-jet shear flows.In particular, the probability y ¢ ( ) r p p , ; p that particles observed at (r, p) originated from r → ∞ with momentum ¢ p determined for both monoenergetic momentum spectra as r → ∞ and for more general spectra ¥ ¢ ( ) f p , 0 (e.g., Figure 8).This allows one to determine the mean value of the in moving from r → ∞ to the observation point r.A second main contribution is the determination of the cosmic-ray flow lines (Equation ( 3)) for both monoenergetic spectra and more general spectra as r → ∞ (e.g., Figures 13-15).
It is important to note that our analysis does not deal with all the complications that play a role in numerical fluid, PIC, or kinetic simulations.For example, we do not delve into the role of the magnetic field, KHIs, Rayleigh-Taylor instabilities, WI (which may be thought of as two fluid beam instabilities), and self-consistent wave or turbulence and particle interactions, which can be addressed by PIC codes.The virtue of our relatively simple mathematical approach is to reveal the gross overall implications of particle transport and acceleration (deceleration) in relativistic shear flows.The flow velocity profile u = u(r)e z is assumed to be such that u(r) is a monotonic decreasing function of r.In this case, there is a 1-1 map between u and r in the shear acceleration region 0 < r < r 2 .The diffusion coefficient κ(r, p) depends on the fluid velocity profile u(r) for 0 < r < r 2 .
Appendix A derives the solution (Equation ( 2)) for f 0 (r, p).The solution is derived using the Fourier transform of f 0 (r, p) with respect to = ( ) t p ln , and uses the adjoint Green's function G ¢ ¢ ( ) r p r p , ; , and Green's formula.In Appendix A, we discuss the asymptotic forms for ¢ ( ) G r p p , ; in the limit as ¢  p p. The solution (Equation ( 2)) for f 0 (r, p) was derived in Webb et al. (2020) using a real Fourier-Bessel eigenfunction expansion.Appendix B concerns the convergence of the eigenfunction expansion and integral representation of the solution as ¢  p p, by using convergence criteria for large N, where N is the number of terms in the series (e.g., Hardy 1946;Bromwich 1947) Rieger & Duffy (2019), Webb et al. (2020), and Wang et al. (2021).In Appendix D, g ˜denotes the relativistic γ of the particle, and γ denotes the relativistic γ of the fluid velocity of the jet.The main result is an estimate of g ãnd total particle energy g = Ẽ m c 0 2 at which the systematic particle energy gains due to shear acceleration balance the particle energy losses.At this point g g = ˜max .The leaky box model includes a source term of −f 0 /τ esc representing the rate of particle escape from the leaky box.For radio jets with a characteristic width of ∼100 pc, and with a magnetic field of 10-30 μG, we obtain Ẽ 10 eV max 15 for the electrons.The existence of an energy or Lorentz gamma g max occurs in models with the turbulence power spectrum of the turbulence µ -( ˜) P k k xx q in the inertial range, where 1 < q < 2 and the scattering mean free time t g µ a ˜, where α = 2 − q. g max does not exist for the leaky box model for the case of Bohm diffusion for the case of q = 1.Nontrivial solutions of the leaky box model in terms of confluent hypergeometric functions were obtained by Wang et al. (2021).
Section 6 concludes with a summary and discussion.

Model and Equations
The relativistic diffusive transport equation for cosmic rays was derived by Webb (1989) and alternative derivations were obtained by Achterberg & Norman (2018a, 2018b).A simplified version of the transport equation was used by Webb (1990) and Webb et al. (2018Webb et al. ( , 2019Webb et al. ( , 2020) ) to describe cosmicray acceleration in relativistic radio-jet shear flows from Webb et al. (2020), the steady-state diffusive cosmic-ray transport equation for relativistic shear flows with u(r) = u(r)e z directed along the axis of the jet reduced to the form is the viscous shear acceleration coefficient due to cosmic-ray viscosity.Here, are the relativistic gamma (γ) and Lorentz β = u(r)/c of the flow, where c is the speed of light.For the case of monoenergetic injection of particles with momentum p 0 at radius r 1 , The transport equation (Equation ( 5)) applies for the shear flow region 0 < r < r 2 , where u = u(r)e z .
Outside the jet shear flow region in r > r 2 , there is no particle acceleration due to cosmic-ray viscosity, and the transport Equation (5) reduces to the diffusion equation where we assume there are no particle sources in r > r 2 at finite r.
In Equations ( 5) and (9) we assume, for the sake of analytical simplicity, that the energetic particle diffusion coefficient κ has the form where τ(r, p) is the particle scattering time.This approximation is reasonable for relativistic particles with speeds of v ≈ c, but is not very good for nonrelativistic particles.For nonrelativistic particles, a more acceptable form for the diffusion coefficient is The scattering time τ(r, p) in the solutions for f 0 (r, p) obtained by Webb et al. (2020) is assumed to have the form where H(x) is the Heaviside step function.Solutions are also possible for different forms for τ(r, p) than Equation 12, but the model Equation (12) for τ(r, p) will be used in the present paper.

Mixed Boundary Conditions
In this section, we briefly discuss the mixed Dirichlet-Von Neumann boundary condition applied at the edge of the jet at radius r = r 2 used by Webb et al. (2020) in their model of particle acceleration by cosmic-ray viscosity in relativistic jet shear flows.
The general solution of the diffusive transport equation (Equation ( 9)) in the region r > r 2 has the form where C(p) and D(p) are integration constants, and we assume that the integral in Equation (13) is finite.By appropriate choice of the integration constants (i.e., C(p) and D(p)), in Equation (13) we obtain the solution for f 0 (r, p) in r r 2 as Setting r = r 2 in Equation ( 14) gives the boundary condition at r = r 2 as The boundary condition (Equation ( 15)) is a mixed Dirichlet-Von Neumann boundary condition on the solution for f 0 (r, p) describing particle acceleration by shear in the jet.In Webb et al. (2018Webb et al. ( , 2019) ) it was assumed that κ → ∞ in the region r > r 2 .In that case, the boundary condition becomes a Dirichlet boundary condition in which f 0 is specified on the boundary at r = r 2 .Webb et al. (2020) showed that the modified boundary condition Equation (15) gives rise to harder accelerated particle spectra at momenta p ? p 0 because particles that have exited the system at r = r 2 could reenter the inner region and become further accelerated by the shear flow in r < r 2 .Taking κ(r, p) in r > r 2 of the form the solution (Equation ( 14)) takes the form The boundary condition (Equation( 15)) then becomes This is the mixed Dirichlet-Von Neumann boundary condition that was applied by Webb et al. (2020) at r = r 2 in their study of particle acceleration by cosmic-ray viscosity in relativistic radio-jet shear flows.This boundary condition is applied to the solution for f 0 (r, p) in 0 < r < r 2 .

Boundary Condition as r → 0
The boundary condition on the solution for f 0 (r, p) as r → 0 (see Webb et al. 2020) is This corresponds to no particle sources as r → 0.

Cosmic-Ray Continuity Equation Perspective
Note that the viscous energization term in the transport Equation (5) can be written in the Fokker-Planck form (Earl et al. 1988;Webb 1989;Achterberg & Norman 2018b), where is the second-order Fermi acceleration momentum dispersion coefficient due to cosmic-ray viscosity, and for the number density of particles in the momentum interval (p, p + dp) about p, the transport equation (Equation ( 5)) can be written in the form Equation (24) can be written more formally in the conservation form where is the radial diffusive particle flux in cylindrical coordinates, and is the mean spatial flow velocity in (r, p)-space.Similarly, defines the mean rate of change of momentum á ñ  p c for the cosmic-ray continuity equation (Equation ( 25)), where the subscript c denotes the continuity equation.Equation (25) describes the balance between particle spatial transport due to changes in radial diffusion and particle momentum due to shear acceleration and cosmic-ray viscosity.
The mean rate of change of momentum á ñ  p c from Equation (28) may be written in the form is the local momentum space spectral index for the cosmic-ray distribution function f 0 .The rate á ñ  p c clearly depends both on the particle momentum space diffusion coefficient D pp as well as on the shape of the particle momentum spectrum via the spectral exponent γ f .
In previous papers (e.g., Rieger & Duffy 2004, 2005a, 2005b;Webb et al. 2018Webb et al. , 2019Webb et al. , 2020)), the rate of the systematic change of momentum 〈Δp〉/Δt has been taken as the representative rate of change of particle momentum, which clearly gives a bad approximation in cases where |γ f | ? 1 (e.g., near a cutoff in the spectrum) and also in cases where |γ f | = 1 where the spectrum is nearly flat.
Note that á ñ  p c is negative for f 0 an increasing function of p and is positive for f 0 a decreasing function of p (i.e., for γ f > 0).Thus, if f 0 is a Gaussian in p about p = p 0 , then the particles with p < p 0 are decelerated and the particles with p > p 0 are accelerated.In the limit that the width of the Gaussian decreases to zero, but the area under the Gaussian is preserved, i.e., for is the derivative of a delta function then the profile of á ñ  p c resembles an N-wave in momentum space, with á ñ <  p 0 c for p < p 0 and á ñ >  p 0 c for p > p 0 .This means that the use of 〈Δp〉/Δt to estimate the particle acceleration time is in general difficult to use, unless there is a good estimate for γ f .In particular the use of 〈Δp〉/Δt to predict the time to accelerate particles from momentum p = p 0 up to momentum p gives a formula in which the particles attain an infinite energy in a finite time (e.g., Webb et al. 2019Webb et al. , 2020) ) which clearly violates causality.Note that the γ f factor is important at late times after a diffusive equilibrium has been established.Corrections to the acceleration time will also occur at early times, when telegrapher equation inertial times should be taken into account (Webb et al. 2018(Webb et al. , 2019(Webb et al. , 2020)).
In more general models including energy losses of particles due to synchrotron losses and second-order Fermi acceleration á ñ  p c in Equations ( 25), (28), and (29) is replaced by the formula represents particle momentum losses due to synchrotron losses (e.g., Blandford 1979;Webb et al. 1984).A similar formula to Equation (33) also applies to particle energy losses due to the inverse Compton effect, which has the same momentum dependence as synchrotron losses.The coefficient D s for synchrotron losses has the form Here, m 0 is the rest mass of the particle species of interest (e.g., m 0 = m e for electrons and m 0 = m p for protons), τ r is the timescale of synchrotron loss, B is the magnetic field strength, and σ T is the Thomson cross section.The basic model in the present paper is mainly concentrated on the steady-state solutions of Equation (5) describing changes in particle momentum due to cosmic-ray viscosity in a cylindrically symmetric model of radio-jet shear flows.A convenient way to visualize the particle transport and acceleration in (r, p) space is by means of integration of the flow line differential equations where á ñ  r c and á ñ  p c are given by Equations ( 27) and (29).
3. Solutions for f 0 (r, p) in Radio-jet Shear Flows From Webb et al. (2020), the distribution function solution of the transport equation for f 0 (r, p), satisfying Equation (5) for 0 < r < r 2 and satisfying the diffusion equation (Equation ( 9)) in the shear-free region r > r 2 is given by an application of Green's formula (90) of Webb et al. (2020).The solutions for shear accelerated particles depend on carrying out transformations of the independent variables (r, p) of the system, to more convenient variables, that lead to analytical solutions of the viscous shear and diffusive transport equation (Equation ( 5)) for f 0 (r, p).We use the variables Here, it is assumed that the fluid velocity parameter β = u(r)/c is a monotonic decreasing function of cylindrical radius r about the axis of the jet.In Equation (36), β 02 is the relativistic relative velocity of the two velocities β 0 and β 2 , where β 0 is the value of β(r) at r = 0 on the axis of the jet, and β 2 is the value of β(r) at the edge of the shear flow region at r = r 2 .The variable ξ(r) in Equation (36) is the relativistic rapidity of the flow β = β(r).The scattering time τ(r, p) used in the model has the functional form (Equation ( 12)) in which there is a change in the radial dependence of τ(r, p) across the boundary at r = r 2 .
Using the variables η = η(r) and ¢ ( ) T p p , defined in Equation (36) the solution for f 0 (r, p) of the transport equation (Equation ( 5)) in the shear flow region 0 < r < r 2 , and in the shear-free flow region r > r 2 , where f 0 satisfies the diffusion equation (Equation ( 9)) and satisfying the boundary conditions (Equation ( 18)) at r = r 2 and the zero source boundary condition (Equation ( 19)) at r = 0, Webb et al. (2020) obtained the solution for f 0 (r, p) in the region 0 < r r 2 of the form is the appropriate Green's function of the system discussed in Equation (38) et seq., below and also in Appendix A.
(i) Case ò ≠ 0 Consider the case where the particles undergo diffusion and scattering in the region r > r 2 outside the shear flow energization region.In this case, the parameter ò = k/s ≠ 0. The solution for f 0 (r, p) for r r 2 is given by Equation (17).Using the boundary condition (Equation ( 18)), one can rewrite Equation (17) in the form  , which represents particles that have scattered both in 0 < r < r 2 and in r > r 2 .This latter component of f 0 (r, p) in r > r 2 decays as ( ) r r s 2 and goes to zero as r → ∞ .The non-singular term ¢ ( )( ) G r p p r r , ; in Equation (39) consists of both accelerated particles with > ¢ p p and decelerated particles with < ¢ p p .(ii) Case ò = 0 and/or ò → 0 In the special case that ò = k/s → 0 or for ò = 0, κ 2 (r, p) → ∞ in r > r 2 .In this case, there is no particle scattering in r > r 2, and hence, The boundary condition at r = r 2 implies that outgoing particles from within 0 < r < r 2 free escape the inner shear flow region and never return after they cross the boundary at r = r 2 because there is no scattering in r > r 2 (i.e., κ → ∞ in r > r 2 ).In this case, for .41 2 These features of the ò = 0 solution can be verified by letting ò → 0 in the general solution (Equations ( 37)-( 39)).
3.1.Solutions for 0 < r r 2 The Green's function ¢ ( ) G r p p , ; in Equation (37) in the shear flow region 0 < r r 2 has the form where the {λ n } and the {j n } satisfy the eigenvalue equation The parameters {χ n } in Equation (42) are given by In the limit as p → ∞ , ~m -¥ G p , where which is the spectral index of viscous shear accelerated particles obtained by Webb et al. (2020).In the limit as ò → 0 (weak scattering in r > r 2 ), j 1 = 2.4048, which corresponds to the first zero of J 0 (x).This limit corresponds to no scattering in r > r 2 , which leads to relatively soft accelerated spectra and weak particle acceleration.In the opposite limit ò → ∞ , j 1 → 0 (see, e.g., Figure 2 of Webb et al. 2020), and μ ∞ → 2a = (3 + α), which is the spectral index of viscous shear accelerated particles obtained by Berezhko (1982) and Rieger & Duffy (2006).
Alternative forms for ¢ ( ) G r p p , ; equivalent to Equation (42) can be obtained by using the eigenvalue equation which is useful in evaluating the solution for G in the limit as ò → 0. This limit corresponds to the case where particles free escape the system in r > r 2 (i.e., κ → ∞ in r > r 2 ).One can also write the solution (Equation ( 42)) for G in the form which is useful for ò ? 1.This corresponds to the case of strong scattering in r > r 2 .In this case, it is difficult for particles to penetrate to r = r 2 from r → ∞ .A hard momentum spectrum of particles is accelerated by the shear flow with μ ∞ → 2a = (3 + α) (because j 1 → 0 in this limit).However, the particle number density is very small at r ≈ r 2 .

Characteristics of G(r, p; p′)
Below we present sample plots of ¢ ( ) pG r p p , ; versus ¢ p p, which are useful in describing the acceleration, deceleration, and redistribution of energetic particles in the shear flow, in which particles injected with momentum ¢ p at r → ∞ end up at position r with momentum p.These results are analogous to plots of ¢ ( ) pG r p p , ; versus ¢ p p studied by Gleeson & Webb (1975) in describing the solar modulation of galactic cosmic rays in the solar cavity, for models of cosmic-ray transport in the solar wind, in which the effect of particle energization at the solar wind termination shock was neglected.
There are some delicate numerical issues in computing ¢ ( ) pG r p p , ; versus ¢ p p for momenta ~¢ p p , where the convergence of the eigenfunction expansions (Equations ( 42)-( 47)) is problematical.A discussion of these issues is provided in Appendix A. A description of the convergence properties of the Fourier-Bessel series solutions (Equations ( 42)-( 47)) and techniques to obtain accurate numerical solutions is provided in Appendix B. Our main purpose in this section is to explore the properties of ¢ ( ) G r p p , , on the physical parameters.Numerical details such as the convergence characteristics of the Fourier-Bessel expansions for ¢ ( ) G r p p , ; are discussed in Appendices A and B. The slow speed of convergence of the Fourier-Bessel series for ¢ p p is implemented using the Van Wijngaarden transformation (van Wijngaarden 1953).
The figure shows that the particles at (r, p) split into two groups, namely, (i) accelerated particles with < ¢ < p p 0 1 , i.e., > ¢ p p and (ii) decelerated particles with ¢ > p p 1 (i.e., < ¢ p p ). Detailed inspection of the solution for ¢ ( ) pG r p p , ; versus ¢ p p in the expanded inset for < ¢ < p p 0.95 1.08 shows that G has a single maximum at roughly ¢ = p p 1.02.From Equation (42), the dominant term in the series ¢ ( ) pG r p p , ; as ¢  ¥ p p is the n = 1 term, which implies 0 as , 49 where Similarly, where a = (3 + α)/2, and μ ∞ is the power-law spectral index of the accelerated particles as p → ∞ .The results (Equations ( 49)-( 52)) are evident in Figure 1.
Figure 2 show plots of ¢ ( ) pG r p p , , versus ¢ p p for different radial positions r/r 2 = 0, 0.6, 1.0, and 2.0.The inset shows plot of ¢ ( ) pG r p p , ; near ¢ p p 1 ( < ¢ < p p 0.9 1.1).The main points to note are (i) ¢ ( ) pG r p p , ; has a single maximum as a function of ¢ p p; (ii) for the case of r = r 2 , pG → ∞ as ¢  p p 1; and (iii) pG is a positive power law in ¢ p p for ¢ < p p 1 and a negative power law for ¢ > p p 1.If one interprets the figure with ¢ p fixed, then the accelerated particles correspond to the region ¢ < p p 1 and the decelerated particles correspond to the region ¢ > p p 1.In Figure 2, we have included a curve for r/r 2 = 2, which corresponds to the region outside the shear flow region 0 < r < r 2 .The figure shows the second term in Equation (39 ) and does not include the d -¢ ( ) p p term in Equation (39).The second term in Equation (43) vanishes as r → ∞ and the first term has the same spectral shape as that at r = r 2 .
Panel (a) in Figure 3 shows how the solution for ¢ ( ) pG r p p , ; varies with ò, for ò = 0.5, 1, 2, 5, 10, and 100 for the case of r/r 2 = 0.6, s = 1, α = 1/3, β 0 = 0.5, and β 2 = 0.The main point to note is the spectral hardening of the accelerated particles as ò increases (see Webb et al. 2020).Again recall that the accelerated particles in Figure 3 correspond to the region ¢ < p p 1. Note that the parameter k in the solution is given by the relation k = òs.
Panel (b) in Figure 3 illustrates essentially the same features of the solution for ¢ ( ) G r p p , ; as panel (a), namely, the hardening of the asymptotic momentum spectrum of the particles at large p for a fixed ¢ p (i.e., the particles with ¢ < p p 1 in the plot).However, in panel (b), k is fixed at k = 10, and the different values of s = k/ò correspond to varying ò.Recall from Webb et al. (2020) that increasing the value of ò corresponds to increasing the role of particle scattering outside the shear layer in the region r > r 2 .
Figure 4 shows plots of ¢ ( ) pG r p p , ; versus ¢ p p for a range of 1, 0.2, 0.3, 0.5, 0.7 and 0.9 48) and p is chosen to be 1 as an example.The black dotted and dashed lines represent the asymptotic power-law spectra given by Equations ( 49) and (51), respectively.And the vertical black solid line at ¢ = p p separates particles into the accelerated (left) and decelerated (right) groups.We see good agreement between the exact and asymptotic solutions.
The other parameters are k = s = 1 (so that ò = 1), α = 1/3 and β 2 = 0.The main point to note is the hardening of the accelerated particle spectrum as β 0 increases (i.e., for particles with < ¢ < p p 0 1 ).Similarly, the spectrum flattens with increasing β 0 for the decelerated particles with ¢ > p p 1.

The Origin of Particles at (r, p)
The Green's function ¢ ( ) G r p p , ; can be used to determine where particles observed at (r, p) originate in the asymptotic momentum spectrum ¥ ¢ ( ) f p , 0 (see Gleeson & Webb 1975 for a similar investigation of cosmic-ray modulation in the solar wind, where particles observed at Earth say, at (r, p), originate with momentum ¢ p as r → ∞, where r corresponds to heliocentric radius).
Following the approach of Gleeson & Webb (1975, 1980), we specify the asymptotic spectrum as r → ∞ , which we refer to as the extragalactic spectrum in the present exposition.We describe the particle spectrum by the differential number density U T as the number of particles with kinetic energy T in the range of (T, T + dT) per unit volume, which is related to the differential number density of particles with momentum p (denoted by U p ) and the differential intensity j T by the formulas gives U T and j T per unit of kinetic energy T in terms of the distribution function f 0 (r, p) at cylindrical radius r and momentum p. From relativistic kinematics and mechanics, the total particle energy E is given by where E 0 = m 0 c 2 is the particle rest mass energy.
Figure 5 shows the cosmic-ray spectrum observed by airshower array experiments.We take ¥ ¢ ( ) f p , 0 , which in part resembles this spectrum.
In the solution (Equation ( 2)) for f 0 (r, p), we need to specify the extragalactic spectrum Before studying realistic and semi-realistic cases, it is instructive to study idealized schematic examples, which illustrate the principles involved.First note that the solution for f 0 (r, p) in Equation (2) can be split up into accelerated . 5 9 We use the notation

G r p p G r p p H p p G r p p G r p p H p p
, ; , ; , , ; , ; , 60 where H(x) is the Heaviside step function, which splits ¢ ( ) G r p p , ; up into accelerated (G A ) and decelerated (G D ) particles.We use the notation Green's functions in Equation (46) observed at different locations: r/r 2 = 0, 0.6, 1.0, and 2.0 (i.e., out of the jet), which are represented by the blue, red, green, and violet solid curves, respectively.Again, the asymptotic spectra for accelerated and decelerated particles are shown in the black dashed and dotted lines, respectively.See the text for more details.
Using Equation (46) for Using Equations ( 59)-( 64) in Equation (58), we obtain the formulae , , e x p , , , e x p .If the asymptotic particle momentum spectrum as particles per unit volume at position r originating in the momentum range of

If one uses kinetic energy ¢
T instead of momentum ¢ p in the above analysis, then the distribution y ¢ T is given by , ; 1 , ; .70 The mean values of ¢ T and ¢ p for fixed (r, p), are given by , depend on both r and p.However, in applications of Equation (71) it is simplest to take r as a constant, in which case á ¢ñ p and á ¢ñ T are functions of p or T and r is a constant parameter.This was the approach taken by Gleeson & Webb (1975) in their study of cosmic-ray modulation in the heliosphere.
The mean value á ¢ ñ ( ) p r p , for particles observed with momentum p at radius r using Equation (71) can be written in the form , 72 . 7 4 0 0 0 In the shear flow region 0 < r < r 2, use of Equations ( 58)-(66) and similar results for f 1 we obtain 75 for the zeroth and first moments in Equations ( 73) and (74).
Using Equations ( 73)-( 75) gives . 76 Similarly, one can determine higher-order moments of ψ p with respect to ¢ p .For example, á ¢ ñ p 2 may be written in the form 77 One can also use Equations ( 76) and (77) to compute the standard deviation of y In the region r > r 2 we find Equations (78)-( 80) can be used to determine á ¢ñ = p f f 1 0 in the region r > r 2 .Similarly, it is possible to calculate á ¢ ñ p 2 and s = á ¢ ñ -á ¢ñ p p 2 2 2 by similar methods.

Examples
As a simple example consider the case where For the spectrum (Equation ( 81)  , in which case we obtain .A similar resonant phenomenon occurs in diffusive shock acceleration (DSA) of energetic particles at astrophysical shocks, when the shock runs over a preexisting power-law spectrum of seed particles with µ m f p 0 c , which have the characteristic spectral index of DSA, namely, μ c = 3r c /(r c − 1), where r c is the shock compression ratio (e.g., Axford 1981Axford , 1994;;Drury 1983;Forman & Webb 1985).In the latter case, the accelerated spectrum of diffusive shock accelerated particles has a spectrum with Axford 1981;Forman & Webb 1985).
It is instructive to note in the solution for f 0 (r, p) in Equation (75) that the nth term in the series (Equation ( 75)) for f 0 may be written in the form

A r p a p a T r p a H p p T r p a H p p
, , , , , , .85 Asymptotics as  ¥ p Our main interest here is the solution for f 0 (r, p) for p > p b and in particular the form of the solution for p → ∞ , which is described by ( ) , namely, where One can also write down expressions for ( ) In the limit as  ¥ p , the asymptotic spectral index μ ∞ for f 0 (r, p) is given by the dominant n = 1 term in Equation (85).Thus, In the examples below, we assume μ ∞ = a + χ 1 < γ s + 2. All particle spectrum as a function of E (energy per nucleus) from air-shower measurements: the high-energy cosmic-ray spectrum scaled by E 2.6 , i.e., E 2.6 × j T vs. E, where E is the total particle energy for the energy range of 10 13 eV < E < 5 × 10 20 eV from a variety of experiments (from Patrignani et al. 2016).The figure shows the spectrum for j T consists of three basic power laws, (a) j T ∼ E −2.6 for 10 13 eV < E < 3 × 10 15 eV, (b) j T ∼ E −3.0 for 3 × 10 15 eV < E < 10 17 eV and (c) j T ∼ E −3.4 for energies 10 17 eV < E < 5 × 10 18 eV followed by (d) the ankle in the spectrum for E > 5 × 10 18 eV.
Analysis of the expression (Equation ( 76)) for á ¢ñ p p gives 89 For p ? p b , the denominator D in Equation ( 89) is dominated by the n = 1 term µ and the numerator N in Equation ( 89) is approximated by Using the approximations (Equations ( 90) and ( 91) , .9 3 For r/r 2 = 0.5 and using the same parameters as in Figure 6, we obtain g c g m The results are consistent with the direct numerical evaluation of á ¢ñ p p from Equation (76).
The solutions for the accelerated particles f A (r, p) and the decelerated particles f D (r, p) are given by Equation (65).The solution for f 0 (r, p) is given by f 0 (r, p) = f A (r, p) + f D (r, p).In principle, there could also be a population of particles f N with no change in momentum, in which case f 0 = f A + f D + f N .We can show that f N = 0 for the particles in 0 < r < r 2 .However, f N ≠ 0 in the outer region r 2 < r < ∞ .In the outer region Some basic characteristics of the solution (Equations ( 58), ( 59 The boundary spectrum as r → ∞ is the cutoff power-law spectrum (Equation ( 81)).For galactic cosmic rays, we take γ s = 2.65 implying f 0 ∝ p −4.65 and j T ∝ p −2.65 .We have in mind the cosmic-ray spectrum in the momentum range of 1 GeV < pc < 10 6 GeV, with p b c = 1 GeV.Thus, the spectrum roughly represents galactic cosmic rays above pc = 10 9 eV and below the cosmic-ray knee at pc = 10 15 eV.Thus, Equation (81) roughly describes the galactic cosmic-ray spectrum below the knee.
The figure shows the accelerated particles ( f A ), the decelerated particles ( f D ), and the total particle distribution function f 0 = f A + f D , as well as the boundary spectrum (Equation ( 81)) as r → ∞ .The inset shows the spectra near the bend at p = p b .At small momenta (p < p b ), f D (the red curve) dominates, but the accelerated particles f A (the blue curve) dominate at large p.The spectra are shown at the fixed radius r = 0.5r 2 where r 2 is the cylindrical radius of the jet.There is significant acceleration of the particles at the higher Figure 6.The spectra of accelerated particles f A (r, p), decelerated particles f D (r. p), and total distribution f 0 (r, p) = f A (r, p) + f D (r, p) of accelerated plus decelerated particles for the case of an asymptotic power-law spectrum f 0 ( ∞ , p) described by Equation (81).The particles f A and f D are defined by Equation (65).See the text for more details.momentum, and the accelerated particle spectrum at large p is flatter than the input spectrum f 0 ( ∞ , p) as r → ∞ .
From Equation (82), the asymptotic momentum spectral index at large momenta (p ?p b ), with µ m For the parameters used in Figure 6, where j 1 is the first root for n = 1 of the eigenvalue equation (Equation ( 45)).For general n, the {j n } depend on ò = k/s, which describes the coupling of the particles in the region 0 < r < r 2 to those in the region r > r 2 .The parameter k describes the shear flow in 0 < r < r 2 and the parameter s describes the radial dependence of the diffusion coefficient in the region r > r 2 .Note that μ D = γ s + 2 is the spectral index of the decelerated particles f D (r, p) for p ? p b , whereas μ A = a + χ 1 is the spectral index for f A for p ? p b .
Figure 7 shows f 0 (r, p) versus p/p b for a range of radii r/r 2 = 0, 0.5, 1.0, and 2.0, for the jet shear flow (Equation ( 95)).The solution for f 0 (r, p) in r > r 2 is given in terms of f 0 (r 2 , p) and f 0 ( ∞ , p) by Equation (38).Note that in r > r 2 , there is a component of the solution , which represents particles that have scattered only in r > r 2 but have experienced no changes of energy.Scattering of particles in r > r 2 and in 0 < r < r 2 leads to a more effective acceleration process compared to solutions with a free escape boundary at r = r 2 (see, e.g., Webb et al. 2020).This effect is incorporated into the model, by means of the mixed Dirichlet-Von Neumann boundary condition at the edge of the jet (see Webb et al. 2020).
Figure 8 shows sample plots of y ¢ ( ) r p p , ; In case (a), there is a particle energy loss (i.e., á ¢ñ > p p). Figure 9 shows a plot of á ¢ñ p versus p for the shear flow model in Figures 7 and 8 for the case of r/r 2 = 0.5.In Figure 9, á ¢ñ ~p p 0.78 0.55 for large enough p.This corresponds to á ¢ñ ~n -( ¯) p p p 0.78 , where ν = (γ s + 2) − μ ∞ and μ ∞ = a + χ 1 is the spectral index of shear accelerated particles as described by the asymptotic formula (Equation ( 92)).In the present case, ν = 0.45, and K = 0.78 in Equation (92).Note that Figure 9 depends on the values of ò and η 2 in the model.
Figure 10 81).Inside the shear flow region, 0 < r < r 2 the solution consists of accelerated and decelerated particles given by Equation (58) et seq.Outside the shear flow region in r > r 2 , f 0 (r, p) is given by Equations (37)-(39).Note that at low momenta, ∂f 0 /∂r is positive, but after the crossover point in the spectra in p > p b the gradient ∂f 0 /∂r is negative.
p > p b (red curve).However, for p < p b , the particles lose energy and -á ¢ñ < p p 0 (blue curve) in this momentum regime.There are particle energy losses to the flow below the breakpoint p b = 1 GeV c −1 for the asymptotic spectrum f 0 ( ∞ , p) as r → ∞ in Figure 7.
Figure 11 shows the zero contour level of ¶ ¶ ( ) .The zero contour is at about p ∼ 3p b at small r at r/r 2 = 10 −2 , but is at p ∼ 2p b at r/r 2 = 10 3 at large r.
Figure 12 shows a contour plot of á ¢ñ ( ) ) p p log 10 in the (r, p) plane.At large p á ¢ ñ ( ) p r p , has contours with á ¢ñ < p p but á ¢ñ > p p at small p at p < p b , where p b = 1 GeV c −1 .The contours are shown using the color code on the right-hand side of the plot.

Flow Lines
There is a temptation to regard the difference in momentum -á ¢ñ p p as the change of the momentum of particles in moving from r → ∞ with momentum ¢ p to the point (r, p), where p is the momentum at position r.However, this interpretation does not account for the local cosmic-ray continuity equation perspective in Section 2.3.From Equations ( 25) and (28) we can write the steady-state continuity equation in the form in regions where there are no sources of particles (i.e., Q = 0).Here, is the momentum spectrum spectral index for f 0 .From Equations (99) and (100), one can visualize p vs. ¢ p for the same model described in Figure 7, for a radio jet with β 0 = 0.9, β 2 = 0.The vertical dashed lines show the value of the mean momenta á ¢ ñ ( ) p r p , and the actual values of á ¢ñ p are given in each panel, after the horizontal dashed key.
the local transport in (r, p) space in terms of the flow line The continuity equation (Equation ( 99)) may be written in the form or alternatively in the form The conservation law (Equation ( 103)) will be automatically satisfied if there exists a potential Λ(r, p) such that In that case, the continuity equation (Equation ( 103)) reduces to the integrability condition which is satisfied if Λ has continuous second-order partial derivatives.From Equations (101) and (104), we obtain . 106 From Equation (106), it follows that the flow line equations possess integrals of the form where const. is an integration constant that has a different value for each flow line.To obtain Λ(r, p) requires further work.Webb & Gleeson (1980) and Gleeson & Webb (1980) investigated cosmic-ray transport in the solar wind (without a termination shock) using flow line differential equations.However, in their analysis, they used particle momentum in an inertial frame fixed in the solar system, rather than the comoving frame particle momentum used in the present analysis.Webb et al. (1983Webb et al. ( , 1985) ) indicated how the flow lines of Webb & Gleeson (1980) are modified in the presence of the solar wind termination shock, where the particles are accelerated by the first-order Fermi mechanism, as well as being decelerated by adiabatic energy losses inside the termination shock.
Flow Lines for r > r 2 .In the model discussed in Section 2, there is no particle acceleration in the region r > r 2 .Thus, in this region The solution for f 0 in r > r 2 is given by Equation (38).Using Equation (38), we find is the diffusion coefficient in r > r 2 and κ 0 = c 2 τ 0 /3 is the value of κ at r = r 2 and p = p 0 .Here, f 0 (r 2 , p) is the distribution function solution at r = r 2 .For f 0 ( ∞ , p) given by Equation (81) f 0 (r, p) is given by Equation (38).From Equations (108) and (109), it follows that = ( ) p const.110 on the flow lines.The direction of the flow lines is parallel to the r-axis and its direction is determined by the sign of f 0 ( ∞ , p) − f 0 (r 2 , p) in Equation (109).Flow lines for 0 < r < r 2 .The flow lines in Equation (101) can be written in the form 112 Using η and ζ as independent variables, the flow line differential equations reduce to the form Flow lines for the Green's function solution ¢ ( ) G r p p , ; in Equation (42).where is the momentum spectral index based on the momentum spectrum for f 0 (r, p).
For the galactic cosmic-ray spectrum for kinetic energies 10 9 eV < T < 10 15 eV, we take f 0 ∝ p −4.65 and taking α = 1/3 corresponding to a Kolmogorov power-law spectrum for the turbulence, gives γ f = 4. 65  in the above kinetic energy range to within 10%.This estimate will in general depend on the cylindrical radius r from the jet axis and also on the value of γ f at different observation points (r, p).
The above estimate of á ñ  p c does not apply to the monoenergetic spectrum solution of Equation (5) for the case in which In this case, one finds that γ f < 0 for p < p 0 (decelerated particles) and γ f > 0 for p > p 0 (accelerated particles).
The result (Equation ( 118)) would of course be different if one included other particle acceleration mechanisms (e.g., second-order Fermi acceleration of particles due to Alfvénic turbulence as in Equation (31)), and particle momentum losses due to synchrotron and inverse Compton losses.

Discussion
An important orientation to the production of UHECRs may be traced back to Hillas (1984), based on the idea that to accelerate particles to a prespecified energy, e.g., E = 10 18 eV for protons say, one needs an acceleration region large enough in which to contain the particle based on the particle gyroradius or mean free path.Also one needs to accelerate the particles up to the required energy in a timescale consistent with the acceleration mechanism invoked, taking into account the timescales of the escape time and energy loss (see, e.g., Rieger 2022 for a recent review of UHECRs in AGN).Here, we reiterate in part the arguments presented by Liu et al. (2017), Webb et al. (2018Webb et al. ( , 2019Webb et al. ( , 2020)), and Rieger & Duffy (2019).
The constraints on the particle acceleration in radio-jet shear flows are (i) the width of the jet ΔL must be greater than the particle mean free path λ and gyroradius r g ; (ii) the particle acceleration time t acc must be less than the timescale of synchrotron loss t synch (i.e., t acc < t synch and other loss timescales); and (iii) the timescale of acceleration t acc < t dyn, where t dyn is the dynamical advection time along the jet, which limits the allowed values of the magnetic field B and the width of the jet ΔL.
The constraints (i)-(iii) on particle acceleration in radio-jet shear flows were investigated in detail by Liu et , where β j is the relativistic β on the axis of the jet, i.e., β j ≡ β 0 = u 0 /c.The bottom panel corresponds to a jet with γ j = 2.0 on the axis of the jet.The top and bottom panels show the lines ΔL = r g by the blue straight line.The condition (i) requires ΔL > r g , which corresponds to the area above the blue line.The condition (ii) that t acc < t synch corresponds to the area below the long dashed red line t acc = t synch .The condition (iii) t acc < t dyn corresponds to the area below the short dashed straight line t acc = t dyn .
One of the main results of this paper is the á ¢ ñ ( ) p r p , relation, which describes how particles observed at (r, p) originated with momentum ¢ p in the asymptotic boundary spectrum as r → ∞ .This is different than the analysis of Webb et al. (2018Webb et al. ( , 2019Webb et al. ( , 2020)), in that the averaging is now over the asymptotic spectrum ¥ ¢ ( ) f p , 0 as well as taking into account that there are many paths in phase space that start at ¥ ¢ ( ) p , and end up at (r, p).There is a probability folded into the trajectory associated with the Green ' based on detailed observations could be used.However, this exercise is beyond the scope of the present paper.
In the present paper, we explicitly studied cases where the particles sources are specified as r → ∞ (i.e., we specify ¥ ¢ ( ) f p , 0 ).It is also possible to develop solutions for f 0 (r, p) in which particle sources are specified at finite values of r at r = r s say, which could be thought of as the primaries.In that case, there is the possibility of investigating secondaries produced by nuclear collisions and spallation in the medium.This investigation is beyond the scope of the present paper.Webb et al. (2019Webb et al. ( , 2020) ) mention that a telegrapher-type transport equation is useful in overcoming the unrealistic noncausal nature of time-dependent momentum space diffusion models because, for example, in the work by Berezhko (1983) on the time-dependent Green's function for particle acceleration by cosmic-ray viscosity, the solution predicts the production of infinite energy particles at a finite time.The same type of behavior is obtained if one naively integrates the equation dp/dt = 〈Δp/Δt〉, as a differential equation in p and t (i.e., there is a blowup in p at a finite time).

Instabilities
In principle, the timescales for various instabilities should be taken into account to see if they are consistent with the particle acceleration process.For example, one would like to know if the KHI growth rate γ KHI is sufficiently large to provide the wave scattering field which is parameterized by the particle scattering time τ(r, p).Here, τ is thought of as a characteristic scattering time for the particles in plasma kinetic theory, and λ = vτ is taken as the mean free path of the particles (e.g., Jokipii 1966).However, in general τ and λ would be expected to be much less than characteristic length and timescales for the KHI or other large-scale fluid instabilities caused by stellar winds and supernova remnants in the interstellar mediums.Ratkiewicz et al. (1994) and Axford (1981Axford ( , 1994) ) argued that the cosmic-ray diffusion coefficient to be used in supernova remnant shock models of cosmic-ray acceleration by the first-order Fermi mechanism should be based on the longer timescales and length scales characteristic of the KHI and Rayleigh-Taylor instability, which are much greater than plasma kinetic scales (the cosmic-ray diffusion coefficient used in simple hydrodynamical cosmic-ray models is usually taken to be κ ∼ 10 28 cm 2 s −1 , in order to take into account the average lifetime of cosmic rays in the galaxy as deduced from 10 Be abundance).The alternative possibility of using the cosmic-ray diffusion tensor from kinetic theory to describe particle acceleration at the supernova forward shock was used by Zakharian et al. (2001).In the latter case, the particle acceleration was much faster at the quasi-perpendicular shock at the equator than at the quasi-parallel shock at the poles.It is perhaps problematical to adequately treat these large-scale phenomena using PIC codes, which concentrate on plasma kinetic instabilities.
Plasma instabilities that could contribute to the particle scattering are the WI (e.g., Wiersma &Achterberg 2004 andAchterberg &Wiersma 2007;Achterberg et al. 2007) and the kink instability (Mizuno et al. 2014).Nishikawa et al. (2020) investigated the role of the WI and the kKHi in relativistic jets in which there is a helical background magnetic field.The second-order Fermi-type dispersion term D pp due to cosmic-ray viscosity in the relativistic diffusive cosmic-ray transport equation of Webb (1989) depends on whether the particle scattering is weak Ωτ ? 1 or strong Ωτ = 1, where Ω is the particle gyrofrequency.The cosmic-ray streaming instability associated with wave-particle interactions (e.g., Lee 1971;Skilling 1971Skilling , 1975;;Bell 1978aBell , 1978b) is invoked in cosmicray acceleration by the first-order Fermi mechanism at nonrelativistic shocks.Achterberg (1979) in a study of Fermi acceleration uses the concept of a Fermi reservoir for both the particles and the turbulence (Burn 1975).This implies that the energy of waves plus particles must be conserved.This approach gives rise to two evolution equations for the energy densities of the turbulence U t and the energy density of the particles U p .The turbulence energy can contribute to accelerating the particles, but it can also go into heating the background thermal plasma.Achterberg invoked the weak turbulence theory of three wave resonant interactions to describe the turbulent cascade.Both Alfvén waves and magneto-acoustic waves are involved in the turbulent cascade.
Bell & Lucek (2001) considered a nonlinear generalization of the cosmic-ray streaming instability in which there are both backward and forward Alfvén wave growth due to work done on the waves due to the cosmic-ray pressure gradient.It is an approximate scheme to describe cosmic-ray modified fluid dynamics and particle acceleration in supernova remnant shocks.This approach neglects wave mixing in which the backward and forward Alfvén disturbances interact due to gradients in the background flow.The Bell & Lucek (2001) model has some connection to the work of McKenzie & Völk (1982) on three fluid models for cosmic-ray modified shocks (see also Ko 1992).Bell (2004) made a more careful analysis of wave instabilities for cosmic rays in supernova remnant shocks.He found that there are both Alfvénic and non-Alfvénic nonresonant modes generated in between the forward and backward shocks in supernova remnant shocks.This instability analysis is sometimes referred to as the Bell instability.These complicated processes are probably best dealt with using PIC  42)) for a monoenergetic spectrum of particles with momentum ¢ p as r → ∞.The arrows indicate the tangent vectors to the flow lines.The dark gray line at r = r 2 is the boundary of the shear flow.The black dashed curve represents á ñ =  r 0 inside the jet.The black solid line starting at r → ∞ with = ¢ p p represents the incoming monoenergetic particles with = ¢ p p .This curve extends into the region 0 < r < r 2 .On this curve á ñ =  p 0. Note that á ñ =  p 0 in the region r > r 2 as there are no changes of energy in this region.
code simulations or fluid simulations (although PIC simulations do not always incorporate all the physics at all scales).Many treatments of instability growth rates are concerned with the wave growth at early times and do not attempt to study wave growth saturation.The instabilities depend on the state of the background medium (e.g., magnetic field strength and orientation; Mach number of the flow; density, etc.).
Rieger & Duffy (2021) explore the constraints imposed by shear-driven instabilities on particle acceleration in relativistic shear flows.They argue that the shear layers in large-scale AGN jets can encompass a sizeable fraction of the jet radius (;0.1), which requires a sizeable injection of electrons for effective acceleration of electrons.They invoke the analysis of Urpin (2002) to characterize the generation of Kelvin-Helmholtz-like modes (shear-driven modes) that can contribute to the turbulence wave power scattering the cosmic rays.Without going into the details, the Rieger & Duffy (2021) analysis was based on Equations ( 26) and (41) in Urpin (2002).They obtained a growth rate of unstable modes ofmath-top="10pt" mathbottom="10pt" where c sj is the thermal gas sound of speed, M sj = V j /c sj is the Mach number of the jet flow, and .Urpin (2002) studied the dispersion equations for the shear flow unstable waves both at short wavelengths using the WKBJ expansion, and at long wavelengths of the order of the jet radius.Urpin (2002) studied long wavelength waves for a jet with a Lorentz Gamma profile of Γ j (r) = Γ 0 [b/((b + r)] 1/3 .The most unstable mode growth rates at small, but nonzero and large enough k have a form similar to that of Equation (119) (see Urpin 2002, Equation (41)).It turns out there are other modes for smaller k in the limit as k → 0, described by Urpin (2002), Equation (38).The unstable waves as k → 0 would scatter UHECRs (see Urpin 2002 and Rieger & Duffy 2021 for further details).

Other Models
Our model for particle acceleration in relativistic radio jets is clearly not the whole story.A more complete model of the acceleration of UHECRs in FR II radio jets has been developed by Seo et al. (2023).They carry out numerical simulations of relativistic hydrodynamical models of FR II radio jets.In their model, they inject cosmic rays continuously into the timedependent radio-jet flow and simulate the transport and acceleration of the cosmic-ray particles using Monte Carlo methods.They identify the principal cosmic-ray acceleration mechanisms, including DSA, relativistic shear acceleration (RSA), which can be split up into acceleration at shear  13.Note that the green curves cut across the á ñ =  p 0 curve horizontally and join continuously with the blue curves.These flow lines correspond to particles that first lose energy, but then gain energy as they cross the á ñ =  p 0 curve to become the red curves.There are also curves on which the particles continuously lose energy (the blue curves).discontinuities in the flow and (RGSA) acceleration in gradual shear flows, and TSA.They find that particles with energy E < 1 EeV gain energy mainly through DSA in the jet spine flow and in the jet backflow (cocoon) containing shocks and turbulence.They hypothesize that particles with energies E greater than a few exaelectronvolts are energized mainly by RSA at the discontinuous backflow interface, reaching energies of the order of 10 20 eV (see also Ostrowski 1998).The socalled expresso mechanism (Caprioli et al. 2015;Mbarek & Caprioli 2019), in which large discontinuous jumps in the particle energy are modeled by Monte Carlo methods is most likely operative.In the present paper, we only study particle acceleration in smooth relativistic shear flows (RGSA).Spitkovsky (2008aSpitkovsky ( , 2008b) ) and Sironi et al. (2013;and references therein) used PIC simulations to study particle acceleration in relativistic shocks.They found that particle acceleration in relativistic shocks is not very efficient in environments with magnetization s , where is the nonrelativistic Alfvén speed.Sironi et al. (2021) used PIC simulations to study particle acceleration and the KHI in relativistic shear flows.Drury & Völk (1981) and Axford et al. (1982) developed two fluid hydrodynamical models of cosmic-ray modified shocks which took into account the reaction of the accelerated cosmic rays on the background flow.These models satisfied the constraints of total momentum and energy conservation, but did not take into account the detailed momentum spectrum of the accelerated particles.Later work by McKenzie & Völk (1982) developed a three fluid model of cosmic-ray modified shocks, consisting of Alfvén waves propagating upstream of the shock that are excited by the cosmic-ray streaming instability, the thermal plasma, and the cosmic-ray hydrodynamical fluid (see also Ko 1992 for a model including both backward and forward Alfvén wave fields and second-order Fermi acceleration of the particles).A recent overview of cosmic-ray hydrodynamical models is given by Zweibel (2017).Zweibel studies both models with (a) self-excited waves and models (b) that have both an external source of waves and self-excited waves.Lemoine (2019) uses the affine connection coefficients of the relativistic transport equation or continuity equation to describe non-inertial forces where the particle momentum is measured in the local fluid frame.This approach can incorporate general relativistic effects due to non-Euclidean metrics in the vicinity of Schwarzschild and Kerr black holes (e.g., Webb 1991Webb , 1992;;Webb et al. 1994).Liodakis et al. (2022) suggest that polarized blazar X-rays are due to particle acceleration in relativistic shocks in radio jets.

Observational Constraints
In principle, one can determine the particle distribution function in the fixed inertial reference frame Σ of the jet by using the Lorentz invariance of the distribution function between the local fluid reference frame S¢ and the fixed reference frame Σ (e.g., Forman 1970).This clearly results in a large anisotropy in the fixed inertial frame, due to the Lorentz boosting between frames (e.g., Lind & Blandford 1985).Lind & Blandford (1985) studied the Lorentz boosting between the local fluid frame and the fixed reference frame, in the interpretation of particle acceleration in the flow.If the flow speed u is much less than the particle speed, the Lorentz boosting effect is referred to as the Compton-Getting effect (Gleeson & Axford 1968;Forman 1970).Lind & Blandford (1985) assumed that to first-order the particle distribution function and synchrotron radiation were isotropic in the fluid frame.In principle, a similar study could be carried out for the present model of particle acceleration by cosmic-ray viscosity in radio-jet shear flows, in which the zeroth, first, and second moments of the distribution function with respect to the solid angle in momentum space in the comoving frame could be determined.This calculation is clearly of interest to observations of radio jets.However, this calculation is beyond the scope of the present paper.Lind et al. (1989) and many other authors have carried out MHD numerical simulations of radio jets.The simulations show marked differences between strongly magnetized and weakly magnetized jets that in general must be considered in more general models of particle acceleration in radio jets (see also Seo et al. 2023).The anisotropy of the distribution function of accelerated particles measured in the fixed inertial reference frame due to relativistic beaming investigated by Lind & Blandford (1985) is clearly of observational interest.
Radio galaxies as well as microquasars emit nonthermal emission in their jets.This requires particle acceleration to high energies to account for the emission.Radio galaxies with jets that are pointed toward us often appear to us as blazars: these are the sources of high-energy neutrinos, implying that these radio galaxies accelerate protons (IceCube, Kun et al. 2014Kun et al. , 2021)).UHECRs also include nuclei such as oxygen and heavier ions (Mayotte et al. 2022;The Pierre Auger Collaboration et al. 2023).Radio galaxies with a selection of relativistic jets aimed at Earth (possibly up to four jets in the case of a supermassive black hole merger (e.g., Gergely & Biermann 2009), which appear as flat spectrum radio sources (Blandford 1979), are prime candidates as a source of UHECRs (Ginzburg & Syrovatskii 1963;Hillas 1984;Biermann & Strittmatter 1987;Jaroschewski et al. 2023).All their activity may be traced to mergers of two central supermassive black holes following a merger of two galaxies (Toomre & Toomre 1972;Press & Schechter 1974;Gergely & Biermann 2009).There is evidence that microquasars or stellar-mass black holes exhibit all these phenomena as well, including mergers and the acceleration of UHECR particles (Mirabel 2004(Mirabel , 2010;;Mirabel et al. 2011;Biermann et al. 2018).There is a general analogy between the activities of stellar mass black holes and supermassive black holes (e.g., Merloni et al. 2003Merloni et al. , 2006;;Falcke et al. 2004;Mirabel 2004;Markoff et al. 2015;Mościbrodzka et al. 2016).
But what is the precise mechanism and location of the accelerated particles?Six pathways of UHECR origin are well established:

The Head Shock of a Jet
Rotating black holes emit from their nearby accretion disk relativistic jets (Blandford & Znajek 1977;Falcke et al. 2004;Daly 2019).These jets are variable (e.g., Her A, Timmerman et al. 2022), and often exhibit a final small region of strong emission, considered to be the head shock (Cyg A, Carilli et al. 1991;Harris et al. 1994: PKS 215369, Gopal-Krishna et al. 2001;Young et al. 2005).The emission there exhibits evidence of weak relativistic motion (3C 179, Porcas 1981), and sometimes a flatter radio spectral index, suggesting fresh injection of newly accelerated particles (Kardashev 1962;McKean et al. 2016;Dabbech et al. 2018;Snios et al. 2018).The ubiquitous cutoff in nonthermal emission seen in such regions at about 3 × 10 14 Hz can be interpreted as being driven by ultrahigh-energy protons (Meisenheimer & Roeser 1986;Biermann & Strittmatter 1987;Roeser & Meisenheimer 1987;Meisenheimer et al. 1989).In this approach, ultrahigh-energy protons in their loss limit drive the wavelengths that in turn excites the irregular magnetic field cascade that governs the electron energy distribution and so defines their cutoff by losses (Kardashev 1962;Matthaeus & Goldstein 1982;Matthaeus et al. 1982;Drury 1983;Goldstein et al. 1995;Chhiber et al. 2021).This model gives a maximum emission cutoff frequency near 3 × 10 14 Hz independent of free parameters, such as the classical electron radius divided by the speed of light multiplied by the electron to proton mass squared to within a simple numerical factor.The corresponding proton energies are maximally around 10 21 eV.But this acceleration site may just reaccelerate energetic protons (or nuclei) injected elsewhere.

Helical Nonthermal Emission in Jets
Very commonly both on arc-second scales as well as on milliarcsecond scales relativistic jets show helical emission structures (like single or double-helix DNA, e.g., Kun et al. 2014;Britzen et al. 2017;Nishikawa et al. 2020;Meli et al. 2023).These structures can be interpreted as helical shock structures; if so, mass continuity requires that such shocks disappear at the central axis, so a spine-sheath structure is foreordained with a dark central spine or axis (M87 Owen et al. 1989Owen et al. , 2000)).Such shocks could explain acceleration and reacceleration along the entire length of a jet.

Re-confinement Shock Structures
Expanding supersonic jets can experience quasiperiodic reconfinement shock structures if the outside pressure exceeds the internal pressure (Sanders 1983;Norman & Winkler 1985).Such structures include conical oblique shocks and often a central Mach disk, plus slip discontinuities going down from the triple point (actually, in jets a triple point ring), where three shocks (two oblique and one perpendicular) meet with a slip discontinuity.As the outside pressure around jets is given by the hot interstellar material in a galaxy (usually an early Hubble-type galaxy, like an elliptical) the outside pressure drops off far faster than the internal pressure, and so jets are usually thought to be self-confined.We discuss the rare exceptions below, when outside ram pressure can kink a jet (Baan & McKee 1985).

Stellar Ablation
Occasionally stars can wander into a jet, and the ram pressure and radiation field can destroy them (e.g., Araudo et al. 2013).Clouds intersecting a jet can analogously provide lots of material (e.g., Gopal-Krishna et al. 2007;Gopal-Krishna & Wiita 2010;del Palacio et al. 2019).Destroying an entire star provides enough material for accelerating particles for quite a while.A single star of mass M å can provide all the material for energetic protons in a jet of power L jet for a time τ given by Obviously, all the associated electrons can also be provided, as well as the heavy elements from the destroyed star.If this were a strong contribution to UHECR composition, then heavy elements might abound, unless these stars were freshly formed.However, even in a starburst induced by the merger of two galaxies (see the radio galaxy Cen A) will not suddenly dominate all the old most common low-mass stars in a galaxy, like the Sun, or even much older stars.

Stellar-mass Black Hole Jets
The activity of stellar-mass black holes has long been studied (Mirabel 2003(Mirabel , 2004(Mirabel , 2010;;Mirabel et al. 2011), and these are usually referred to as microquasars because of their many similarities (Fender 2001;Fender et al. 2003;Falcke et al. 2004;Markoff et al. 2015).As has been long known, active stellar-mass black holes are quite common; since they derive from massive stars, with a new black hole born in our Galaxy about every 400 yr (Diehl et al. 2006(Diehl et al. , 2010)).Most of them are expected to be in binary systems (Chini et al. 2012(Chini et al. , 2013a(Chini et al. , 2013b)), allowing binary black hole mergers, now commonly detected in gravitational waves (LIGO/VIRGO/KAGRA lists GWTC-1, GWTC-2.0,GWTC-2.1, and GWTC-3).There is a proposal that we have seen one such merger in the radio structure of the starburst galaxy M82 (Kronberg et al. 1985;Bartel et al. 1987;Biermann et al. 2018) around the compact radio source 41.9+58 believed to have been a gamma-ray burst (Muxlow et al. 2005).A large fraction of massive stars are in triple or even quadruple systems, allowing second-generation mergers (i.e., three mergers in sequence: two mergers of four black holes to two black holes in a first generation, and the resulting merged black holes merge again in a second generation to form one final black hole), leading to quite high final masses of the black hole.

Shear Injection and Acceleration
As relativistic jets propagate through an environment, there is friction with particles both entering and escaping the jet.The outside medium is the halo gas of the galaxy, detectable in X-rays, and with a density and pressure rapidly decreasing toward the outside (Biermann et al. 1982;Biermann & Kronberg 1983;Roediger et al. 2015;Sheardown et al. 2018Sheardown et al. , 2019;;Ge et al. 2021;Geng et al. 2022); further outside there is the gas of the Local Group of galaxies or the cluster of galaxies with a much larger radial pressure decay length scale, but with a very low pressure.That implies that these jets are usually selfconfined, and do not interact much with the outside.Early-type galaxies are usually in groups and clusters of galaxies, so the jets easily propagate into the intergalactic medium, which is a hot tenuous gas.As these galaxies and their associated jets move around, there can be strong ram pressure for the jets.
However, there are several instances where such an interaction is directly visible, e.g., in the radio galaxies NGC 1265 and in 3C 449 (Baan & McKee 1985;Feretti et al. 1999;Croston et al. 2003).The jets show a kink, where they just bend into a new direction; the interpretation is that ram pressure in an intergalactic medium outside their host galaxy pushes them to the side, by inducing a highly oblique shock in the jet.Now inside the jet the medium is thought to be highly relativistic, and outside the medium is visibly hot, since we can detect it in X-rays.It follows that the density of the outside medium must be much higher than the inside medium, providing a very strong injection (see Equation (15)).However, there is no evidence of any enhanced emission downstream from that kink, so the injection of new particles must be weak.In this paper, we focus on very powerful radio galaxies, with a much higher jet power than those of either NGC 1265 or 3C 449; so there is even less motivation to expect a strong injection of particles.
However, the outside pressure could get quite large in the core of the host galaxy.Relativistic jets are often thought to start as pure Poynting flux jets, with an electron/positron plasma only (V404 Cyg, Siegert et al. 2016).Such jets need to add hadronic material early to explain the neutrino observations (Kun et al. 2021) that show that the emission is optically thick for energetic photons at the epoch of neutrino emission.Shear is a suitable candidate, and may be the only chance to inject hadronic material quickly.
Therefore, shear injection and acceleration may be strong near the foot of the jet.

Concluding Remarks
The acceleration of cosmic rays in astrophysical shear flows was initiated by Berezhko (1981Berezhko ( , 1982Berezhko ( , 1983) ) and by Berezhko & Krymskii (1981).Earl et al. (1988) re-derived the diffusive cosmic-ray transport of Parker (1965; see also Krymskiy 1964 for a constant solar wind speed V; Gleeson &Axford 1967 andDolginov &Toptygin 1967) by including the effects of cosmicray viscosity as an important mechanism for particle acceleration in nonuniform fluid flows.They also included the effects of advection, diffusion, adiabatic compression, and accelerating flows on the particle transport for the case of nonrelativistic fluid flows.The generalization of the diffusive transport equation for cosmic rays for relativistic flows was derived by Webb (1989; see also Achterberg & Norman 2018b).Pitch angle transport equations for cosmic rays for nonrelativistic flows were derived by Skilling (1971Skilling ( , 1975)).Webb (1985Webb ( , 1987) ) and Achterberg & Norman (2018a) derived the corresponding focused transport equations for relativistic flows.Particle transport and acceleration in relativistic shear flows were studied by Webb (1989Webb ( , 1990; see, e.g., Rieger 2019 for a recent review).
The basic shear acceleration model for particle acceleration in radio-jet shear flows was described in Section 2. In Section 2.3, we introduced a perspective on particle acceleration and transport based on writing the transport equation in conservation form.This allows one to consider the particle transport in (r, p) space in terms of flow lines, which are the solutions of the associated differential equation system = á ñ á ñ   dr dp r p in the (r, p) plane, where á ñ  r and á ñ  p are the mean particle advection velocity and rate of change of momentum.This approach was used by Gleeson & Webb (1974, 1980) and Webb & Gleeson (1980) to describe cosmic-ray transport in the heliosphere.However, in the latter papers, p corresponded to the particle momentum in an inertial reference frame fixed in the solar system.The changes of particle energy in these papers correspond to the work done on the cosmic rays by the solar wind due to the cosmic-ray pressure gradient (see also Jokipii & Parker 1967 for a hydrodynamical version of this result).In the present analysis, we use the particle momentum p in the comoving frame, moving with the background flow, and hence the changes of particle energy also contain non-inertial reference frame effects.
The application of flow line differential equations in (r, p) space to the steady-state particle transport equation in radio-jet shear flows was further explored in Section 4.3.The rate of change of particle momentum f derived from the cosmic-ray continuity equation involves the spectral slope parameter g = - ¶ ¶ f p ln ln f 0 of the particle distribution function f 0 (r, p).The rate of change of momentum á ñ  p was compared with the rate of change of the systematic momentum 〈Δp/Δt〉 in Section 4.4.If γ f > 0, the particles gain energy, but if γ f < 0 the particles lose energy.In other words, the changes in particle energy and momentum in this formulation depend on the sign of the spectral index γ f .Section 3 gives the detailed form of the Green's function ¢ ( ) G r p p , ; for the model.Plots of ¢ ( ) pG r p p , ; versus ¢ p p provide insight into the inverse particle transport problem, in which (r, p) are regarded as fixed, and ¢ p corresponds to the particle momentum as ¢  ¥ r .Section 4 discussed the origin of particles of momentum p at position r for representative asymptotic spectra r .Detailed analysis of the general solution (Equation ( 2)) and the Green's function ¢ ( ) G r p p , ; forms (Equations ( 42), ( 46), and (47)) are given in Appendix A. Appendix A develops a Fourier transform form of the solution, which can be investigated using Cauchy's theorem to obtain Equation ( 42).An alternative method to obtain ¢ ( ) G r p p , , was developed by Webb et al. (2020) in which a Fourier-Bessel expansion was developed, which is analogous to a Fourier series solution, except that the usual trigonometric functions are replaced by orthogonal Bessel eigenfunctions.The development in Webb et al. (2020) only uses real analysis.The work in Equation (A33) in Appendix A is useful in obtaining approximate asymptotic solutions in the limit that ¢  p p 1.In this limit, the usual eigenfunction expansion is difficult to evaluate.A detailed discussion of the convergence properties of the eigenfunction expansion (Equation ( 42)) and the equivalent expansions (Equations ( 46) and (47)) are described in Appendix B.
The main motivation for the present paper is the observation that the radio-jet shear acceleration model has steady-state solutions for f 0 (r, p) of the form is the effective Green's function and r is cylindrical radius about the jet axis (see Equation (2)).From Equation (122), one obtains 123 as the fraction of particles at (r, p) that originated in the momentum interval ¢ ¢ + ¢ ( ) p p dp , as ¢  ¥ r (see Figure 8).Gleeson & Webb (1975) and Li et al. (2009) studied a similar inverse problem for solar modulation of galactic cosmic rays in the solar cavity.Li et al. (2009) used a set of stochastic differential equations for representative particles starting at (r, p) and integrated the stochastic particle paths backward in time to the boundary at r = R.The spectrum of ¢ ( ) f R p , 0 versus ¢ p may be constructed from the ensemble of particle paths.Gleeson & Webb (1975) and Li et al. (2009) studied the energy losses of galactic cosmic rays with momentum ¢ p as ¢  ¥ r that are observed at radius r with momentum p in the heliosphere.
In this paper, we have studied the inverse transport problem of the shear acceleration model for particle acceleration in radio jets, in which one fixes the observation point (r, p) and determines y ¢ ¢ ( ) r p p dp , ; p defined by Equation (123), which in turn is used to determine the mean momentum á ¢ ñ ( ) p r p , .One of the main results of the paper is the dependence of á ¢ ñ ( ) p r p , on p for a fixed r illustrated in Figures 8 and 9, which apply to protons for which synchrotron and inverse Compton losses can be neglected.Figure 9 shows that á ¢ñ ~p p 0.78 0.55 for the semi-realistic, cutoff power-law spectrum f 0 ( ∞ , p) displayed in Figure 7, corresponding to a jet with β 0 = 0.9, , where p b is the break in f 0 ( ∞ , p) in Figure 7 (recall that ò = 0 corresponds to a free escape boundary at r = r 2 , but ò → ∞ corresponds to strong scattering in the shearfree region r > r 2 ).Representative values of á ¢ñ p are given in Figure 8.For example, p = 10 GeV c −1 gives á ¢ñ = p -3.07 GeV c 1 .Thus, p is substantially larger than á ¢ñ p in this example (i.e., substantial particle acceleration is obtained).For p = 10 6 GeV c −1 , p = 1553 GeV c −1 .Using the relation á ¢ñ ~p p 0.78 0.55 , we deduce that for E ∼ pc = 10 18 eV, protons correspond to á ¢ñ = ´p 6.95 10 eV c 13 1 as r → ∞ .These calculations give some idea of the original average momentum of the particles given by á ¢ñ p that gave rise to the particles observed at r = 0.5r 2 with p ∼ 10 18 eV c −1 .Thus, the acceleration process boosts the original average particle momentum á ¢ñ p by a factor of á ¢ñ ~ṕ p 1.44 10 4 in this case.The contours of á ¢ñ p in the (r, p) plane in Figure 12 give a more general view of á ¢ñ p p as a function of r and p for other values of the observation radius r and momentum p.These results should be consistent with the constraints of Figure 16 for an acceptable source.However, Figure 16 implicitly assumes that the particles are constrained to lie in the shear acceleration region 0 < r < r 2 .We identify ΔL = r 2 as the width of the jet shear acceleration region.This constraint should be relaxed somewhat in the case where there is significant particle scattering in the region r > r 2 , because these particles can be scattered back into the shear acceleration region 0 < r < r 2 .Figure 16, shows that to constrain and accelerate an E = 10 18 eV particle in the radio-jet shear flow, only the gray-shaded region of (ΔL, B) space is allowed.For example, to obtain E = 10 18 eV protons requires that 1 kpc < ΔL < 10 5 kpc for B = 1 μG and for a radio jet with γ j = 1.1 (i.e., β j = 0.4166).For B = 100 μG, ΔL lies in the range of 10 pc < ΔL < 10 3 kpc.The length of the jet d = 100ΔL was assumed.AGN jets with these characteristics for ΔL and B (Liu et al. 2017;Webb et al. 2019) include the possible sources MKN 501 and MKN 421 (Dermer 2007;Sahayanathan 2009;Abbasi et al. 2014;Caprioli et al. 2015).
Rieger & Duffy (2019), Webb et al. (2020), and Wang et al. (2021; see also Liu 2015 andLiu et al. 2017) developed a leaky box model for particle acceleration in radio-jet shear flows, that takes into account (i) particle acceleration by the shear flow, (ii) particle energy losses due to synchrotron and inverse Compton losses, and (iii) particle escape from the acceleration region due to spatial transport.Rieger & Duffy (2019) showed that for a spectrum of turbulence of the form P xx (k) ∝ k − q in the inertial range, for which the particle scattering time is t g µ a ˜, where α = 2 − q, there exists a particle energy g g = ˜max at which the particle energy gains balance the particle energy losses due to synchrotron and inverse Compton processes.At energies g g < ˜max the particle energy gains dominate whereas at energies g g > ˜max the particle energy losses dominate.For electrons, simple estimates (Rieger & Duffy 2019;Wang et al. 2021) imply that g = Ẽ m c e max 2 10 eV 15 for a range of jet models, whereas for protons efficient acceleration to several exaelectronvolts appears feasible.An overview of these leaky box models for particle acceleration in radio-jet shear flows is given in Appendix D.
The effect of synchrotron losses on the spectral index μ ∞ of particle acceleration in a leaky box model of radio-jet shear flows for Bohm diffusion was studied by Webb et al. (2020).For Bohm diffusion, there is no maximum energy g max for which synchrotron losses balance particle energy gains due to shear acceleration.They obtained the formula where the χ e and χ p are the values of χ for electrons (e) and protons (p).Thus, the timescale of proton synchrotron loss is much larger than the timescale of electron loss by a factor of ( ) m m p e 4 .Synchrotron losses can be neglected for protons but must be taken into account for electrons.
Explicit expressions for t acc , t synch , and t esc are given by Webb et al. (2020).Webb et al. (2020) Equations ( 163) and (167) can be written as The escape time is given by Webb et al. (2020), Equation Equations ( 127) and (128) apply to both electrons and protons.The timescale of the synchrotron loss for electrons in the model, given by Webb et al. (2020), (Equation 162) is equivalent to ´- In Equation (127), κ(p) = c 2 τ(p)/3 is the particle diffusion coefficient in the shear acceleration region 0 < r < r 2 .
The accelerated particles in the radio-jet shear flow should in general have a gyroradius of r g < ΔL, in order to be contained in the shear flow (however, this constraint can be relaxed somewhat if the particles in r > r 2 = ΔL can scatter back into the shear flow region (i.e., for the case of ò ≠ 0)).For the above parameters used by Webb et al. (2020), an E = 10 16 eV electron has r g = 3 × 10 19 cm, for B = 1 μG but ΔL = 3 × 10 20 cm, which is about 100 pc.In this example, r g /ΔL ≡ r g /r 2 = 0.1 and the electron is well confined by the jet.However, for E e = 10 18 eV, r g = 3 × 10 21 cm, then r g /r 2 = 10, and in that case, the electron is not well confined by the jet.Again, if there is significant scattering of the particle in r > r 2 the confinement criterion is then not very strict.
It is important to note that recent simulations of FR II jets by Seo et al. (2023) indicate that shear acceleration appears to dominate shock acceleration, at energies E > 1 EeV.Gradual shear acceleration and non-gradual shear acceleration are increasingly important above about 1 EeV.This seems to justify the present model as a particle acceleration model above about E = 1 EeV, but the model is not so good at lower energies where DSA is dominant.
A further aspect of the gradual shear acceleration problem studied in this paper is the flow line equation solutions of Equation (3) associated with the conservation form of the transport equation (Equation (4)) in (r, p) space.This analysis gives an overall picture of the particle radial transport due to spatial diffusion and the corresponding particle momentum changes over the whole (r, p) phase space domain (see, e.g., Sections 2.3 and 4.3 and Figures 13 and 14 illustrating flow lines for a monoenergetic particle spectrum boundary condition for f 0 as r → ∞ ).b .In the latter case, the flow lines are inward at low momenta at 0 < p < p b and move outward at p > p b .The particles in this case always gain energy (in reality, there must be both gains and losses of the particle momentum because there is momentum space diffusion describing the shear acceleration).The flow lines show that on average the particles always gain energy in a statistical sense in the Figure 15 example.In principle, there should be energy losses at high energies due to synchrotron losses, but these are not included in the present analysis.The energy losses are clearly masked in Figure 15, by the weighting of the results by the spectrum ¥ ¢ ( ) f p , 0 .Observational constraints and particle acceleration possibilities (e.g., radio-jet emission beaming, microquasars, blazars, stellar-mass black hole jets, supermassive black hole jets following the merger of two galaxies, shock acceleration at the head shock, at helical shocks in jets with helical magnetic fields, re-confinement shocks, stellar ablation, injection, and acceleration of particles in shear flows) were discussed in Section 5.These observational constraints are important in developing more general models and theories of particle acceleration in radio jets than the simple shear acceleration model used in the present paper.However, the present model does have some heuristic value in elucidating particle acceleration by cosmic-ray viscosity in radio-jet shear flows.For example, it is clearly easier for particles to enter the shear flow and cross back and forth across the flow, than for particles to swim against the flow as envisaged for particle acceleration by the first-order Fermi mechanism in relativistic shocks.In early models of particle acceleration by the first-order Fermi mechanism at relativistic shocks based on the focused transport equation for cosmic rays (e.g., Kirk & Schneider 1987a, 1987b;Achterberg et al. 2001) the pitch angle diffusion coefficient D μμ used in the models was assumed to be of a given form suggested by quasi-linear theory.For particles to be accelerated by the first-order Fermi mechanism in these models the particles need to be highly beamed along the magnetic field and shock normal in order to swim upstream of the shock.Such highly beamed distributions are inherently unstable, giving rise to waves that isotropize the distribution function leading to quenching of the shock acceleration process.This suggests that a more realistic model of wave-particle interactions and momentum space diffusion would take into account the dynamical coupling of the waves or turbulence scattering the particles and the particle transport coefficients.One might expect both energy losses and energy gains of particles due to momentum space diffusion and perhaps quasi-isotropic particle distribution functions due to enhanced turbulence, but with a net particle energization in such models.The PIC simulations of Sironi et al. (2013) show that particle acceleration in relativistic shocks is not very efficient in environments with magnetization of σ = B 2 /(μ 0 n i m i c 2 ) > 1. Kirk et al. (2023) investigate particle acceleration at ultrarelativistic perpendicular shock fronts.Put another way, a more realistic model for particle acceleration at relativistic shocks should take into account momentum space diffusion due to either cosmic-ray viscosity or due to turbulence.Momentum space diffusion would then lead to both energy gains and energy losses at relativistic shocks.
Particle acceleration in radio jets should in principle take into account all other possible acceleration mechanisms, e.g., by first-order and second-order Fermi mechanisms, shear acceleration, electric fields at shocks, the effects of accelerating reference frames (pseudo-gravity), and particle energy losses.However, these issues are beyond the scope of the present paper.These more general models are necessary to properly evaluate the significance of each acceleration mechanism.171)-( 172) of Webb et al. (2020) or Equations ( 124) and (125), for the case of Bohm diffusion (α = 1 and τ ∝ p) for electrons for case (i) including synchrotron losses (blue curve) and for case (ii) with no synchrotron losses (black curve with χ = 0).β 2 = 0, s = 0.1, k = 1, ò = k/s = 10, r 2 /(100 pc) = 1, B/10 μG = 0.1, and P 0 = (δB/B) 2 = 0.01.supported in part by NASA New Horizons funding through contract NAS5-97271/Task order 30.

Appendix A
In this appendix, we present a version of the analysis of Webb et al. (2019Webb et al. ( , 2020)), which shows how the solution (Equation ( 37)) for f 0 (r, p) arises from the Fourier transform solution, for general cases where ò ≠ 0 or ò = 0.This analysis was not presented in Webb et al. (2020).Webb et al. (2019) only considered the case of ò = 0.The Fourier transform of the solution is useful in evaluating the Green ' as .A1 There are two distinct Green's function solutions of the transport equation (Equation ( 5)).The solution for ¢ ( ) G r p p , ; satisfying the boundary condition (Equation (A1)) is referred to as the monoenergetic spectrum Green's function solution.Another fundamental Green's function solution corresponds to steady injection of particles into the flow at radius r = r 1 with momentum p = p 0 , which are energized in the shear flow region 0 < r < r 2 , where f 0 (r, p) satisfies the transport equation (Equation ( 5)) in the region 0 < r < r 2 and f 0 satisfies the diffusion equation (Equation ( 9)) in the shear-free region r > r 2 .The solution in 0 < r < r 2 is matched at the boundary at r = r 2 by the mixed Dirichlet-Von Neumann boundary condition (Equation ( 18)) in which we set f 0 ( ∞ , p) = 0. We also require that the solution corresponds to no particle sources at r = 0, which is equivalent to the requirement that rS d → 0 as r → 0, where S d is the radial particle diffusive flux (see Equation ( 19)).
The transport equation (Equation ( 5)) in the region 0 < r < r 2 can be expressed in the form (e.g., Webb et al. 2020, Equation (28)) The Green's function of Equation (A2) can be obtained by using the Fourier transform (note there is a misprint in Equation (125) of Webb et al. 2020 where dν should be replaced by dT).Taking the Fourier transform of Equation (A2) with respect to T gives the ordinary differential equation The homogeneous equation (Equation (A5)) has solutions of the form where η = ξ 0 − ξ, and I 0 (z) and K 0 (z) are modified Bessel functions of the first and second kind of order zero (Abramowitz & Stegun 1965, 9.6.1 p. 374).
Using standard results from Morse & Feshbach (1953) we obtain the solution for f0 in the form A 1 0 The results (Equations (A9) and (A10)) generalize the solution for n¢ ¯( ) f r, 0 given by Webb et al. (2020), Equations (131) and (132), which were obtained for ò = 0. Webb et al. (2019) only considered solutions with ò = 0, in which there was a free escape boundary at r = r 2 .
The Green's function solution of Equation (A2) can now be obtained by inverting the Fourier transform solution (Equation (A9)), i.e., ò The solution for f 0 (r, p) can be obtained from Equation (A11) using complex integration, accounting for the branch cuts for the multivalued function n L ¢ ( ) and by taking into account the residue theorem and the poles of the integrand (Equation (A9)) that occur at the zeros of D(Λη 2 ) = 0.This process gives the standard eigenfunction expansion (Equation( 46)) of Webb et al. (2020), namely, The parameters η 1 = η(r 1 ) and η 2 = η(r 2 ).J 0 (z) and J 1 (z) are ordinary Bessel functions of the first kind.The eigenvalues j n satisfy the eigenvalue equation (Equation ( 43)), and the parameters λ n and χ n are related to the j n via the equations  5)).In this approach, Equation (5) is written in the form where Γ s is the shear acceleration coefficient (Equation ( 6)).In the case of monoenergetic injection of particles with momentum p = p 0 at radius r = r 1 , is the source term in Equation (A15).The monoenergetic source Green's function of Equation (5), is given by Equation (46) of Webb et al. (2020).Using integration by parts, one obtains Green's theorem for the transport equation (Equation (A15)) in the differential form   where we assume (r, p) lies inside the region  in Equation (A21).The solution for f 0 (r, p) taking into account the mixed Dirichlet boundary conditions at r = r 2 in the region  reduces to To determine G ¢ ¢ ( ) r p r p , ; , notice that Equation (A15) for f 0 with monoenergetic source term Equation (A16) has the same form as the adjoint Green's function (Equation (A19)).The map of Equation (A15) onto Equation (A19) is obtained by using the transformations in which the Green's function solution is described by Equations (A9)-(A12).Using Equation (A25) in Equations (A9)-(A11), we obtain the Fourier space version of G ¢ ¢ ( ) r p r p , ; .as is the normalization constant.In addition the map (Equation (A25)) implies the transformation A 2 8 0 0 0 Using Equations (A26), (A27), and (A11), we obtain in Equation (A24) we obtain for 0 < r < r 2 the solution form Equation (A30) gives the solution for ¢ ( ) G r p p , ; only for the region 0 < r < r 2 .The solution for ¢ ( ) G r p p , ; for r > r 2 is given by Equation (39).
Because Λ is an even function of n¢, Equation (A30) can be written as a Fourier cosine transform A.2. ò ≠ 0 Asymptotics as ¢  p p In this section, we give a sketch of the asymptotics of the solution for ¢ ( ) G r p p , , as ¢  p p for the case of ò ≠ 0. Similar, but different results apply for ò = 0 case.In the case of ò ≠ 0, Φ(w, w 2 ) from Equation (A32) for large w and w 2 has the approximation The simplest approximation in Equation (A34) is to drop the higher-order terms that are of order O(1/w) in the numerator and denominator leading to the crude approximation .From Equation (B3), the A n (r) in the series (Equation (B2)) can be approximated by On the other hand, in the large n limit, we can immediately see that and therefore, , where M − N + 1 is the number of terms summed over in the range of N n M.This is similar to the harmonic series, which diverges like ~N ln .
Next, we seek the asymptotic form of A n (r = 0).In this case, . Hence, the use of the approximation (Equation (B4)) on J 0 ( j n ) yields Note that G N (0) now asymptotically becomes an alternating series and the Leibniz criterion ensures its convergence (Arfkin 1985, Section 5.3).Finally, we consider the general r-dependent case (i.e., h < < ¯( ) r 0 1).Again, with Equation (B3) and the approximation (Equation (B4)) for both h ( ¯) J j n 0 and J 0 ( j n ), we find suggesting that G N (r) grows like ∼N 2 for r close to r 2 .On the other hand, when we look at the h = ¯0 case (i.e., at r = 0), it turns out that .Similarly, for h < < 0 1, we can find the G N (r) in Equation (B2) does not converge like that in the ò ≠ 0 are the Fourier transforms in Equation (C7), where t = ¢ - t t.From Equation (C4), we obtain Thus, the general solution (Equation( 2)) for f 0 (r, p) reduces to the Fourier space (Equation (C7)).From Equation (C7), we obtain The systematic shear acceleration rate is given by g g g g t and the corresponding acceleration time is given by g The cooling rate due to synchrotron losses and inverse Compton scattering is The expression (Equation (D13)) for g max may also be obtained by setting the the cooling time for the cutoff in the electron energy spectrum.Further discussion of g max and E max is given by Wang et al. (2021).The value of g max in general will depend on the radial profile of γ j (r).The g  ( ˜) F q , are defined in terms of Kummer's confluent hypergeometric function as ∞ .Previous work by Webb et al. studied only the Green's function solution for ¢ ( ) G r p p , ;

. 53 0 Figure 1 .
Figure 1.The blue curve shows the Green's function in Equation (42) (or equivalently, Equation (46) or (47) for ò = 1.Other parameters used here are listed in Equation (48) and p is chosen to be 1 as an example.The black dotted and dashed lines represent the asymptotic power-law spectra given by Equations (49) and (51), respectively.And the vertical black solid line at ¢ = p p separates particles into the accelerated (left) and decelerated (right) groups.We see good agreement between the exact and asymptotic solutions.

Figure 3 .
Figure 3. Panel (a) Green's functions (Equation (46)) for different ò for a fixed spatial dependence on τ outside of the jet (i.e., s = 1, see Equation (12)); Panel (b) the same as panel (a) but with different s and a fixed spatial dependence on the jet profile (i.e., k = 10, see Equation (36)).
p b , and γ s are constants.Here, ¢ = ¢ p p p b and = p p p b and the subscript b refers to the bend in the spectrum at p = p b .

Figure 4 .
Figure 4. Green's functions in Equation (46) with different central jet speed β 0 .See the text for more details.

Figure
Figure5.All particle spectrum as a function of E (energy per nucleus) from air-shower measurements: the high-energy cosmic-ray spectrum scaled by E 2.6 , i.e., E 2.6 × j T vs. E, where E is the total particle energy for the energy range of 10 13 eV < E < 5 × 10 20 eV from a variety of experiments (fromPatrignani et al. 2016).The figure shows the spectrum for j T consists of three basic power laws, (a) j T ∼ E −2.6 for 10 13 eV < E < 3 × 10 15 eV, (b) j T ∼ E −3.0 for 3 × 10 15 eV < E < 10 17 eV and (c) j T ∼ E −3.4 for energies 10 17 eV < E < 5 × 10 18 eV followed by (d) the ankle in the spectrum for E > 5 × 10 18 eV.
), (65), and (66)) in which the W ( ) Equations (82) and (83) are displayed in Figure6.The calculations pertain to a relativistic jet with parameters b b a GeV c 1 , for a radio jet with the same parameters as given in Figure 7 for k = s = 1, ò = k/s = 1, α = 1/3, β 0 = 0.9, β 2 = 0, γ s = 2.65, and p b c = 1 GeV, at observation radius r/r b = 0.5 and a range of observation momenta p (here p b is the break momentum in the spectrum as r → ∞ ).The figure shows six panels with different momenta p at the observation radius r/r 2 = 0.5, which is located inside the jet shear flow region.The different panels correspond to the momenta (a) p = 0.1, (b) p = 1, (c) p = 10, (d) p = 100 (e) p = 10 3 , and (f) p = 10 6 GeV c −1 .The plots show the numerical value of the mean momenta á ¢ ñ ( ) p r p , in each panel.The mean momenta are also shown by the dashed vertical lines.It turns out that á ¢ñ < p p in panels (b)-(f), but á ¢ñ > p pin panel (a).Thus. the particles in general are energized on average by the shear flow in traveling from their source at r → ∞ to the observation point at r = 0.5r 2 .
shows a plot of the momentum difference -á ¢ñ | | p p versus p for the solution.Note that -á ¢ñ | | p p is an increasing function of p for pc > p b c = 1 GeV.Using Equations (92)-(94), one can see that -á ¢ñ > p p 0 for p > p b , i.e., the particles on average gain energy from the shear flow for

Figure 7 .
Figure7.Spectra for the distribution function f 0 (r, p) vs. p/p b for the solution, Equation (37).The particles are accelerated (decelerated) in the shear flow in the region 0 < r < r 2 about the jet axis.The boundary spectrum of particles as r → ∞ , namely ¥ ¢ ( ) f p , 0 is given by Equation(81).Inside the shear flow region, 0 < r < r 2 the solution consists of accelerated and decelerated particles given by Equation (58) et seq.Outside the shear flow region in r > r 2 , f 0 (r, p) is given by Equations (37)-(39).Note that at low momenta, ∂f 0 /∂r is positive, but after the crossover point in the spectra in p > p b the gradient ∂f 0 /∂r is negative.
p) plane.The particles diffuse inward at small p, where ¶

Figure 9 .
Figure 9. Left panel: a plot of á ¢ñ p vs. p for the model of shear acceleration described in Figure 7. Shown is the exact numerical value of á ¢ñ p .The dashed curve represents the asymptotic approximate á ¢ñ p vs. p curve (Equation (92)) as p → ∞. á ¢ñ ~p p 0.78 0.55 .The right panel shows the same information using a plot of á ¢ñ ( ¯) p p 0.78 0.55 vs. p.Here, r/r 2 = 0.5.

Figure 10 .
Figure 10.Plots of -á ¢ñ | | p p vs. p for the model for particle acceleration described in Figures 7-9.Note that the particles are accelerated from momentum ¢ p to momentum p in the case of á ¢ñ -< p p 0 (color red).The particles are decelerated for á ¢ñ -> p p 0 (color blue).

pp
coefficient for particle acceleration in the radio-jet shear flow.From Equations (114) and (115), Comparing this rate with the rate of change of the continuity equation momentum (Equation (29))

Figure
Figure 12.Contour plot of á ¢ ñ ( ( ) ) p r p p log , 10 in the (r, p) plane for the model in Figure 7.Note á ¢ñ < p pat large p, but á ¢ñ > p pat small p.The parameters are the same as in Figure 7.

Figure 13 .
Figure13.Flow lines for the Green's function (Equation (42)) for a monoenergetic spectrum of particles with momentum ¢ p as r → ∞.The arrows indicate the tangent vectors to the flow lines.The dark gray line at r = r 2 is the boundary of the shear flow.The black dashed curve represents á ñ =  r 0 inside the jet.The black solid line starting at r → ∞ with = ¢ p p represents the incoming monoenergetic particles with = ¢ p p .This curve extends into the region 0 < r < r 2 .On this curve á ñ =  p 0. Note that á ñ =  p 0 in the region r > r 2 as there are no changes of energy in this region.
of an acceleration (i.e., [g] = [L][T] −2 ).It should be noted that g is not the acceleration vector of the fluid.In Equation (119),

Figure 14 .
Figure14.More details of the flow lines shown in Figure13.Note that the green curves cut across the á ñ =  p 0 curve horizontally and join continuously with the blue curves.These flow lines correspond to particles that first lose energy, but then gain energy as they cross the á ñ =  p 0 curve to become the red curves.There are also curves on which the particles continuously lose energy (the blue curves).
Figure 15 illustrates the form of the flow lines for the case of a more general spectrum ¥
concomitant of ψ and f 0 .It turns out that the differential operator  is self-adjoint, i.e., is equivalent to Equation (90) ofWebb et al. (2020).Once the adjoint Green's function G¢ ¢ ¢ ( ) r p r p , ; , has been determined, Equation (A23) gives the solution for f 0 (r, p) in the form of Equation (37), where is a large number), according to the Dirichlet's test (Dirichlet's test is described in Wikipedia; however, a more substantial discussion of convergence and divergence test for series and infinite integrals were developed byBromwich 1947 andHardy 1946).B.2.Case of ò = 0Compared with Equation (B5), in this case, the asymptotic form of the eigen equation J 0 (x) = 0 reads as (seeEquation (B4))Consequently, a new large n-limit approximation for j n is obtained, to the approximation (Equation (B7)).For h = ¯1 at r = r 2 , the partial sum G N (r) in Equation (B2) (or in fact Equation (B1) degenerates to G = 0 due to the eigen equation J 0 ( j n ) ≡ 0. However, if we consider the asymptotic behavior of G N (r) as r → r 2 , we can expand h

n
which will result in a diverging G N (0), or equivalently,

whereFrom
r as fixed.At first glance, Equation (C11) looks a bit strange because it seems to imply that the radial dependence of the numerator and denominator must cancel.On the other hand, Equation (C7) makes sense because y n¢ ¯( ) r, is the Fourier transform of the Green's function propagator.The inverse Fourier transform of Equation (C11) in principle gives the required possibility is beyond the scope of the present paper.A formal derivation of Equation (C9) follows from the definition of Equation (C8), namely, delta distribution in Equation (C11).Appendix DIn this appendix, we describe the leaky box model for cosmic-ray acceleration in radio-jet shear flows developed byRieger & Duffy (2019),Webb et al. (2020), andWang et al. (2021).The model includes the acceleration of particles due to radio-jet shear flows, synchrotron, and inverse Compton cooling for electrons and particle escape from the leaky box.The basic transport equation for f 0 has the form g ˜is the relativistic particle Lorentz gamma (for ultrarelativistic particles, the particle speed v ≈ c and the particle , we follow the approach and the notation ofWang et al. (2021)ʼs leaky box model for electron acceleration in radio-jet shear flows.The transport equation (Equation (D1)) can be written in the Fokker˜t c is the synchrotron and inverse Compton cooling rate for g ˜and τ esc is the escape time from the leaky box, with source pg ˜Q 4 2 .τ esc is given by the estimate Δr = L is the radius of the cylindrical box and κ = c 2 τ/3 is the particle diffusion coefficient.The energetic particle momentum diffusion coefficient gg ˜D and the systematic acceleration rate g áD D ñ ˜t sh have the formLiu et al. (2017),Rieger (2019), and Webb et al.  (2018, 2019, 2020)  assume a scattering mean free time of the form = pc/(eB) is the particle gyroradius and Ω B = eB/ (mc) is the gyrofrequency.Here, ℓ b is the outer scale or correlation length, where < < k k 0 b defines the energy containing range of the turbulence, and k b = 2π/ℓ b is the break in the turbulence power spectrum at = k k b .In the inertial range of assume the power spectrum is negligible in the dissipation range of > k k d .
f = U rad /U B is the energy density ratio between the target field (U rad ) from inverse Compton scattering and U B = B 2 /(8π) is the magnetic field energy density.Here σ T is the Thomson scattering cross section.The cooling rate (Equation (D13)) balances the systematic shear acceleration rate at g g acceleration time t acc in Equation (D12).Using Equation (D7) for A 0 and Equation (D13) for A 2 in Equation (D15) gives the equation the particle Lorentz gamma, where ξ = (δB/B) 2 measures the strength of the turbulence.The formula gives the Lorentz particle g ˜at which energy gains due to shear acceleration balance energy losses due to synchrotron and inverse Compton losses.It is equivalent to formula (11) of Rieger & Duffy (2019).with γ j = 3 and with a linear decreasing profile of β j = u j /c as a function of r, namely dβ j /dr = − β 0 /Δr, Taking the outer scale ℓ b = Δr as the lateral width of the shear layer, now gives the equation (12) of Rieger & Duffy (2019) for g max .For B = 30 μG and Δr = 100 pc, g = 3μG and Δr = 0.1 kpc and assumes ξ = (δB/B) 2 = 0.2.For B = 10 μG, Δr = 100 pc, Equation (D20) gives g ~1 Applications of the solutions (Equations (D21)-(D24)) by Wang et al. (2021) for g ( ˜) n to X-ray emission in the FR I jet of Centaurus A (NGC 5128) and to the FR II jet of 3C 273 were studied.Wang et al. (2021) used an r-averaged form of Γ s in which .Webb https:/ /orcid.
Figure 16 from Webb et al. (2019) shows plots of Webb et al. (1984Webb et al. ( , 2020), ),e, t acc , t synch , and t esc refer to the characteristic acceleration time, the synchrotron loss time, and the leaky box escape time.Webb et al. (2020)considered the case of Bohm diffusion for which α = 1.From Equation (170) ofWebb et al. (1984Webb et al. ( , 2020), ), The approximation (Equation (B7)) is extremely important in studying the convergence of G N (r).
n in the limit of large n.