Particle Acceleration by Cosmic Ray Viscosity in Radio-jet Shear Flows

The origin and acceleration of cosmic rays to ultra-high energies has been the subject of many papers in cosmic ray astrophysics (e.g., Hillas 1984; Biermann & Strittmatter 1987; Axford 1992, 1994; Protheroe & Szabo 1992; Biermann 1993; Rachen & Biermann 1993; Stanev et al. 1993; Norman et al. 1995; Blandford et al. 2014).

Acceleration mechanisms for energetic particles in radio jets include second-order Fermi acceleration (e.g., Burn 1975; Achterberg 1979; Eilek 1979; Bicknell & Melrose 1982) and first order Fermi acceleration at shocks (e.g., Axford et al. 1977; Krymsky 1977; Bell 1978a, 1978b; Blandford & Ostriker 1978; Drury 1983; Malkov & Drury 2001).

The physics of the formation of radio jets has been a subject of intense interest since the 1970s (e.g., Blandford & Rees 1974; Blandford & Zjanek 1977; Blandford & Payne 1982; Kommisarov 1990; Kommisarov & Falle 1998; Narayan et al. 2007), involving both MHD models and gravitational physics of black holes and accretion discs.

Caprioli (2015) and Mbarek & Caprioli (2019) argue for a so-called one-shot, "expresso" mechanism involving a single particle interaction with the background plasma flow as a possibility to obtain EeV particles in radio jets, in which the particles are Lorentz boosted by a factor of γ² in energy by a single interaction with the relativistic jet flow, where γ is the Lorentz factor of the flow. This process is distinct from a diffusive energization process.

Matthews et al. (2019) present 2D and 3D simulations of radio jets. They argue that the acceleration of ultra-high-energy cosmic rays (UHECR) can in principle be produced by diffusive shock acceleration (DSA) in the back-flowing material of radio-jet galaxy lobes (e.g., in Faranoff–Riley type II jets). Araudo et al. (2016, 2018) argue that the observations of radio galaxy hotspots suggest that jet termination shocks cannot accelerate cosmic rays to EeV energies.

The acceleration of energetic particles at relativistic shocks due to the first order Fermi mechanism was developed by Kirk & Schneider (1987a, 1987b). To describe particle acceleration at relativistic shocks, it is necessary to take into account the pitch-angle anisotropy of the particle distribution in the vicinity of the shock, where the particle speed and the fluid speed are comparable. This necessitates the use of the pitch-angle focusing transport equation for energetic particle transport. For non-relativistic flows, this equation is related to the drift kinetic transport equation (e.g., Kulsrud 1983; le Roux et al. 2007; le Roux & Webb 2009; Webb et al. 2009), and the focused transport equation for cosmic rays in which the particles are scattered by backward and forward Alfvén waves was derived by Skilling (1975) (see also le Roux et al. 2007; le Roux & Webb 2009). The relativistic version of Skilling's equation was derived by Webb (1985, 1987b; see also Achterberg & Norman 2018a). The pitch-angle focusing equation for relativistic flows for a parallel shock configuration was derived by Kirk et al. (1988). Kirk & Duffy (1999) and Achterberg et al. (2001) give reviews of particle acceleration at relativistic shocks, in which the focused transport equation or the related stochastic differential equations are solved for the particle distribution function.

In addition, it is important to take into account particle energy losses in many applications (e.g., synchrotron losses and inverse Compton losses for electrons; Webb et al. 1984; Blasi 2010) and nuclear spallation and interaction of the particles with the 3 K microwave background near the GZK cut-off (Greisen 1966; Zatsepin & Kuzmin 1966; Abbasi et al. 2008); the original version of the GZK cut-off at about 5 × 10¹⁹ eV assumed that the cosmic ray spectrum at these energies is dominated by protons, which appears not to be the case (Aab et al. 2017a, 2017b).

A further mechanism for the acceleration of energetic particles due to shear in the background flow and cosmic ray viscosity was developed by Berezhko (1981, 1982, 1983), Berezhko & Krymsky (1981), Earl et al. (1988), Jokipii et al. (1989), Jokipii & Morfill (1990), Webb (1989, 1990), Webb et al. (1994), Ostrowski (1990, 1998), Ostrowski (2002), Stawarz & Ostrowski (2002), Rieger (2002), Rieger & Duffy (2004a, 2004b, 2005, 2006, 2016), Ohira (2013), Liu (2015), Liu et al. (2017), Kimura et al. (2018), and Webb et al. (2018a). The analysis is carried out using the particle momentum p in the mean scattering frame or fluid frame, and consequently non-inertial effects must be taken into account, as the mean scattering frame is in general a non-inertial frame.

Other acceleration mechanisms involving particle acceleration in regions containing magnetically reconnecting magnetic fields and flux ropes (e.g., Drake et al. 2013; Zank et al. 2014, 2015; Guo et al. 2015; le Roux et al. 2015, 2016; Li et al.2017) have been developed to explain spacecraft observations of energetic particles in the vicinity of shocks and current sheets in the heliosphere (e.g., Khabarova et al. 2017).

Lemoine (2019) provides a statistical approach to generalized Fermi acceleration, in which the affine connection coefficients and Ricci rotation coefficients are used to describe the particle equation of motion in a local inertial frame (tetrad frame) moving with the background flow, with applications to first order Fermi acceleration at shocks, centrifugal and shear acceleration (cosmic ray viscosity), and unipolar inductive acceleration in the vicinity of black holes. This approach is distinct from (but related to) the moment equation approach and the focused transport equation approach, based on the relativistic Vlasov equation or Boltzmann equation used by Webb (1985, 1989), Kirk et al. (1988), and Achterberg & Norman (2018a, 2018b).

Mostafavi et al. (2017, 2018) and Mostafavi & Zank (2018) have incorporated the effects of viscosity into two-fluid hydrodynamic and magnetohydrodynamic models of pick-up ion modified shocks in the heliosphere and in the very local interstellar medium. These models are analogous to two-fluid cosmic ray modified shock models (e.g., Drury 1983; Malkov & Drury 2001), but with the addition of energetic particle viscosity included in the models. These models lead to a unique shock structure, rather than the earlier two-fluid models of cosmic ray modified shocks, which could result in up to three downstream states for a given upstream state (i.e., including energetic particle viscosity or pick-up ion viscosity leads to a unique shock structure).

A recent paper by Webb et al. (2018a; hereinafter referred to as paper I), on energetic particle acceleration by cosmic ray viscosity in radio-jet shear flows, derived momentum spectra of accelerated and decelerated particles due to scattering back and forth across the flow. The fluid velocity of the jet was assumed to be a cylindrically symmetric flow along the jet axis (the z-axis) of the form: u = u(r) e_z, where u(r) was a monotonic decreasing function of cylindrical radius r about the jet axis. For the sake of analytical simplicity, the distribution function f₀(r, p) of the accelerated particles was assumed to be independent of distance z along the jet axis. The model neglected the effect of energy changes due to adiabatic compression, particle acceleration due to the first order Fermi mechanism at shocks in the flow, and the effect of an accelerating flow on the particle energy changes. In fact, ∇_αu^α = 0 (zero flow divergence) and zero acceleration vector for the flow (${\dot{u}}_{\alpha }={u}^{\beta }{{\rm{\nabla }}}_{\beta }{u}_{\alpha }=0$) are consequences of the assumed shear flow profile used in the model. The model was based on the diffusive cosmic ray transport equation of Webb (1989) for particle transport in relativistic flows. Achterberg & Norman (2018a, 2018b) have recently derived similar diffusive transport equations for cosmic ray transport in relativistic flows. The effect of second-order Fermi acceleration on the particle energization and synchrotron losses were neglected in order to obtain an analytically tractable model. The aim of the paper was to illustrate the energy changes of cosmic rays due to cosmic ray viscosity and shear in radio-jet shear flows.

The aim of the present paper is to further explore this model for particle acceleration in radio-jet shear flows, by using a modified form of the fluid shear and particle scattering time τ (r, p) than that used in paper I. The mean scattering time τ (r, p) used was of the form:

$$\begin{eqnarray}&&\tau (r,p)={\tau }_{0}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }\displaystyle \frac{g(\xi )}{{rd}\xi /{dr}},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&\xi (r)=\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{1+\beta (r)}{1-\beta (r)}\right),\end{eqnarray} \tag{ 2 }$$

is the relativistic rapidity and β = β(r) = u(r)/c is the relativistic β of the flow. In paper I, g(ξ) was taken to be a constant. The net result of this choice was that in most cases studied for monotonic decreasing functions ξ(r) and β(r) of r, the model implied $\tau (r,p)\to \infty$ as r → 0. In this paper, we consider an alternative form for the function g(ξ), namely

$$\begin{eqnarray}&&g(\xi )=k(\xi -{\xi }_{0}),\end{eqnarray} \tag{ 3 }$$

where k is a positive constant. By a judicious choice of ξ(r) (or alternatively β(r)) as a monotonic decreasing function of r, we obtain solutions in which $\tau (r,p)={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$ for the scattering time throughout the jet. Thus τ(r, p) is finite throughout the jet.

In Section 2 we describe the model and diffusive transport equation for cosmic rays in relativistic flows. Section 3 derives Green's function solutions of the transport equation for the choice g(ξ) = k(ξ(r) − ξ₀) where ξ₀ = ξ(0), for the case of monoenergetic injection of particles into the flow with momentum p = p₀ at radius r = r₁, and with free escape of particles from the jet at the outer radius r = r₂ (i.e., 0 < r₁ < r₂ is assumed).

Section 4 uses Green's formula for the transport equation to derive solutions in which a monoenergetic spectrum of particles is specified at the outer boundary at $r={r}_{2}$ and with free escape of particles at r = r₂ and with $p\ne {p}_{0}$. Section 4 also describes more general solutions in which an arbitrary distribution function ${f}_{0}({r}_{2},p)$ is specified on the boundary of the jet at radius r = r₂.

Section 5 describes the effects of the cosmic ray inertia, via a second-order time derivative in a space independent transport equation modeling the particle transport and acceleration, by both cosmic ray viscosity and also by second-order Fermi acceleration by Alfvénic turbulence in the background flow (e.g., Skilling 1975). We include particle momentum losses due to synchrotron losses. The spatial transport and escape of particles from the acceleration region is described by a −f₀/τ_esc source term (this model is somewhat analogous to the leaky box model used to describe secondary and primary cosmic rays in the galaxy). The model emphasizes the role of the characteristics of the second-order telegrapher equation that describes the particle transport. In particular we use the characteristics and bi-characteristics of the transport equation to obtain formulae for the time for particles to be accelerated from momentum p = p₀ up to momentum p based on the bi-characteristics of the transport equation. This gives an upper limit to the particle momentum p = p_max(t) that can be achieved in this model due to momentum space diffusion.

The method to derive the Green's function solutions involving J₀(z) Bessel functions is described in Appendix A. More general Green's function solutions involving J_m(z) Bessel functions (m > 0) are described in Appendix B. In this paper, we only investigate the Green's function solutions and monoenergetic spectrum solutions involving the J₀(z) Bessel functions.

In Appendix C we discuss the modification of the boundary condition for the transport equation at the outer radius of the shear flow at $r={r}_{2}$, for the case where the particle spectrum ${f}_{0}(\infty ,p)$ is specified as $r\to \infty$, and in which the particles have a finite diffusion coefficient in r > r₂, but there is no acceleration of the particles in the region r > r₂. In the most general case, the boundary condition at r = r₂ is a mixed Dirichlet–Von Neumann boundary condition, which also involves the specification of the boundary spectrum as $r\to \infty$. In the simplest case the boundary condition reduces to a Dirichlet boundary condition for f₀(r₂, p) in the limit as $\kappa \to \infty$ in the region r > r₂. In the present paper we only consider solutions that correspond to Dirichlet boundary conditions at r = r₂.

In Appendix D we describe the characteristics and bi-characteristics of the model transport Equation (63), which is a second-order partial differential equation for ${f}_{0}(p,t)$ with independent variables t and p. This analysis is used in Section 5 to describe the maximum momentum p = p_max(t) that the particles can achieve from momentum space diffusion due to shear acceleration (cosmic ray viscosity) and due to second-order Fermi acceleration owing to Alfvén waves.

Appendix E presents an analysis of the steady-state balance between energy gains of particles due to cosmic ray viscosity in radio-jet shear flows, and synchrotron losses, in the context of a leaky box model of particle acceleration. For simplicity we only consider the case where the escape time ${\tau }_{\mathrm{esc}}(p)\to \infty$. The analysis considers the dependence of the solution on the mean particle scattering time $\tau (p)={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$, and in particular the dependence of the solution on α for the cases (a) 0 < α < 1, (b) α = 1 (Bohm diffusion), and (c) α > 1.

Appendix F corrects some formulae from Webb et al. (2018a) Section 5, on constraints on acceleration of 1 EeV protons in radio jets similar to the analysis of Liu et al. (2017).

Section 6 concludes with a summary and discussion.

It is assumed that the particle scattering is sufficiently strong that an isotropic scattering and diffusion model may be used. The model uses the diffusive transport equation for cosmic rays in relativistic flows derived by Webb (1989). The scattering wave frame for the cosmic rays is taken as being coincident with the background plasma flow frame (i.e., the co-moving fluid frame). We neglect the effects of second-order Fermi acceleration. The transport equation for the isotropic distribution function of the particles f₀(x^α, p) (x^α = (ct, x, y, z) is the spacetime position four vector of the particle) has the form Webb (1989):

$$\begin{eqnarray}\begin{array}{rcl} & & {{\rm{\nabla }}}_{\alpha }\left[{{cu}}^{\alpha }{f}_{0}-\kappa \left({g}^{\alpha \beta }+{u}^{\alpha }{u}^{\beta }\right)\left(\displaystyle \frac{\partial {f}_{0}}{\partial {x}^{\beta }}-{\dot{u}}_{\beta }\displaystyle \frac{{\left({p}^{0}\right)}^{2}}{p}\displaystyle \frac{\partial {f}_{0}}{\partial p}\right)\right]\\ & & +\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left[-\displaystyle \frac{{p}^{3}}{3}{{cu}}_{;\beta }^{\beta }{f}_{0}+{p}^{3}{\left(\displaystyle \frac{{p}^{0}}{p}\right)}^{2}\kappa {\dot{u}}^{\beta }\left(\displaystyle \frac{\partial {f}_{0}}{\partial {x}^{\beta }}-{\dot{u}}_{\beta }\displaystyle \frac{{\left({p}^{0}\right)}^{2}}{p}\displaystyle \frac{\partial {f}_{0}}{\partial p}\right)\right.\\ & & -\left.{\rm{\Gamma }}\tau {p}^{4}\displaystyle \frac{\partial {f}^{0}}{\partial p}\right]=Q.\end{array}\end{eqnarray} \tag{ 4 }$$

For the case of strong scattering (ωτ ≪1, where ω = qB/m is the gyro-frequency and τ is the scattering time), the viscous energization coefficient Γ is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{c}^{2}}{30}{\sigma }_{\alpha \beta }{\sigma }^{\alpha \beta },\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{\sigma }_{\alpha \beta }={u}_{\alpha ;\beta }+{u}_{\beta ;\alpha }+{u}_{\beta }{\dot{u}}_{\alpha }+{u}_{\alpha }{\dot{u}}_{\beta }-\displaystyle \frac{2}{3}{u}_{;\gamma }^{\gamma }\left({g}_{\alpha \beta }+{u}_{\alpha }{u}_{\beta }\right)\end{eqnarray} \tag{ 6 }$$

is the shear tensor of the background flow (a more conventional form of the shear tensor is one half of the tensor in (6)). More general forms of the transport equation for both weak and strong scattering are described in Webb (1989).

In the application of the transport Equation (4) to particle acceleration in radio jets, we use the model of Webb et al. (2018a), in which the jet fluid velocity is assumed of the form:

$$\begin{eqnarray}&&{\boldsymbol{u}}=u(r){{\boldsymbol{e}}}_{z},\end{eqnarray} \tag{ 7 }$$

where u(r) is a monotonic decreasing function of r. For the fluid velocity profile (7), the acceleration vector of the fluid ${\dot{u}}_{\alpha }=0$ and the divergence of the fluid velocity four vector ${u}_{;\beta }^{\beta }\equiv {{\rm{\nabla }}}_{\beta }{u}^{\beta }=0$. The net result is a model in which the particles are accelerated solely due to cosmic ray viscosity and shear of the background flow. The scattering time τ(r, p) for particles with momentum p, as measured in the fluid frame, is described by (1) and (2).

Note that the viscous energization term in the transport Equation (4) can be written in the standard Fokker–Planck form:

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{\partial }{\partial p}\left({p}^{4}{\rm{\Gamma }}\tau \displaystyle \frac{\partial {f}_{0}}{\partial p}\right)=\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{f}_{0}\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle \right)\\ \,\,-\displaystyle \frac{{\partial }^{2}}{\partial {p}^{2}}\left[\left\langle \displaystyle \frac{{\left({\rm{\Delta }}p\right)}^{2}}{2{\rm{\Delta }}t}\right\rangle {p}^{2}{f}_{0}\right]\end{array}\end{eqnarray} \tag{ 8 }$$

(Earl et al. 1988; Webb 1989; Achterberg & Norman 2018b), where

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle =\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{4}{\rm{\Gamma }}\tau \right)\end{eqnarray} \tag{ 9 }$$

is the systematic acceleration term and

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\left({\rm{\Delta }}p\right)}^{2}}{2{\rm{\Delta }}t}\right\rangle ={\rm{\Gamma }}{p}^{2}\tau \end{eqnarray} \tag{ 10 }$$

is the average rate of momentum dispersion. From (9) and (10), we note the relationship

$$\begin{eqnarray}&&{p}^{2}\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle =\displaystyle \frac{\partial }{\partial p}\left[\left\langle \displaystyle \frac{{\left({\rm{\Delta }}p\right)}^{2}}{2{\rm{\Delta }}t}\right\rangle {p}^{2}\right],\end{eqnarray} \tag{ 11 }$$

between the systematic acceleration coefficient $\langle {\rm{\Delta }}p/{\rm{\Delta }}t\rangle$ and the momentum dispersion coefficient $\langle {\left({\rm{\Delta }}p\right)}^{2} \rangle /2{\rm{\Delta }}t$.

From Webb et al. (2018a) the viscous shear acceleration coefficient for the model is given by

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{c}^{2}}{15}{\gamma }^{4}{\left(\displaystyle \frac{d\beta }{{dr}}\right)}^{2}=\displaystyle \frac{{\gamma }^{4}}{15}{\left(\displaystyle \frac{{du}}{{dr}}\right)}^{2}.\end{eqnarray} \tag{ 12 }$$

For steady-state solutions of the transport Equation (4), we assume there is no z-dependence to the solution. In this case, (4) reduces to the form,

$$\begin{eqnarray}&&-\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(\kappa r\displaystyle \frac{\partial {f}_{0}}{\partial r}\right)-\displaystyle \frac{{c}^{2}}{15}{\gamma }^{4}{\left(\displaystyle \frac{d\beta }{{dr}}\right)}^{2}\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{4}\tau \displaystyle \frac{\partial {f}_{0}}{\partial p}\right)=Q,\end{eqnarray} \tag{ 13 }$$

where β = u/c, and the source term Q has the form:

$$\begin{eqnarray}&&Q=\displaystyle \frac{{N}_{0}}{8{\pi }^{2}{p}_{0}^{2}{r}_{1}}\delta (p-{p}_{0})\delta (r-{r}_{1}).\end{eqnarray} \tag{ 14 }$$

Equation (9) can be used to estimate the acceleration time for particles of momentum p. Assuming τ = τ₀p^αh(r) and a mean free path λ = vτ (v is the particle speed in the fluid frame or scattering frame), (9) gives the expression

$$\begin{eqnarray}&&{t}_{\mathrm{acc}}^{m}=\displaystyle \frac{p}{\langle {\rm{\Delta }}p/{\rm{\Delta }}{t}_{m}\rangle }=\displaystyle \frac{15}{(4+\alpha ){\gamma }^{4}{\left({du}/{dr}\right)}^{2}\tau },\end{eqnarray} \tag{ 15 }$$

in the scattering frame (i.e., using time as measured in the scattering frame, or fluid frame). In the inertial frame,

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{acc}}^{I} & = & \displaystyle \frac{p}{\langle {\rm{\Delta }}p/{\rm{\Delta }}{t}_{I}\rangle }=\displaystyle \frac{15}{(4+\alpha ){\gamma }^{3}{\left({du}/{dr}\right)}^{2}\tau }\\ & \equiv & \displaystyle \frac{15v}{(4+\alpha ){\gamma }^{3}{\left({du}/{dr}\right)}^{2}\lambda },\end{array}\end{eqnarray} \tag{ 16 }$$

where we have used the result Δt_I/Δt_m = γ for the scattering process taking place in the fluid frame. There are limitations on the use of (16) to estimate the mean particle momentum of particles accelerated in the shear flows (e.g., Webb et al. 2018a).

The steady-state Green's function of (13) describes the particle transport and acceleration in the region 0 < r < r₂, where particles are injected monoenergetically into the flow with momentum p = p₀ at radius r = r₁ (0 < r₁ < r₂). It is assumed that there are no sources of particles on the axis of the jet at r = 0. The boundary conditions on the solution for f₀(r, p) (see paper I) are

$$\begin{eqnarray}&&{f}_{0}\to 0\quad \mathrm{as}\quad r\to {r}_{2},\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{rS}}_{d}=r\left(-4\pi {p}^{2}\kappa \displaystyle \frac{\partial {f}_{0}}{\partial r}\right)\to 0\quad \mathrm{as}\quad r\to 0.\end{eqnarray} \tag{ 18 }$$

The source term (14) describes steady injection of particles with momentum p = p₀ at radius r = r₁. The model can be applied to the case of a "naked" jet, in which there is no back-flowing cocoon (u(r) > 0), or with a back-flowing cocoon with u > 0 near r = 0 and monotonically decreasing with u < 0 at large r (i.e., the case of a back-flowing cocoon; e.g., Webb et al. 2018a).

To obtain analytical solutions of the transport Equation (13), we use the new independent variables,

$$\begin{eqnarray}&&\xi (r)=\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{1+\beta }{1-\beta }\right)={\tanh }^{-1}(\beta ),\quad T=\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right),\end{eqnarray} \tag{ 19 }$$

and assume a scattering mean free time τ(r, p) of the form,

$$\begin{eqnarray}&&\tau ={\tau }_{0}\ {\left(p/{p}_{0}\right)}^{\alpha }g(\xi )/[r(d\xi /{dr})],\end{eqnarray} \tag{ 20 }$$

where g(ξ) is an arbitrary function of ξ and α is constant determining the momentum dependence of mean scattering time τ. Note that ξ is a monotonic increasing function of β. With these choices, the transport Equation (13) with source term (14) reduces to

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{1}{g(\xi )}\displaystyle \frac{\partial }{\partial \xi }\left(g(\xi )\displaystyle \frac{\partial {f}_{0}}{\partial \xi }\right)+\displaystyle \frac{1}{5}\left(\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {T}^{2}}+(3+\alpha )\displaystyle \frac{\partial {f}_{0}}{\partial T}\right)\\ & & \quad =\displaystyle \frac{3{N}_{0}\delta (T)\delta (\xi -{\xi }_{1})}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}g({\xi }_{1})}.\end{array}\end{eqnarray} \tag{ 21 }$$

Here ξ₁ ≡ ξ(r₁) is the value of ξ at r = r₁ (see Webb et al. 2018a).

At this point we choose an appropriate functional form for g(ξ) that determines the scattering time τ(r, p) in (20). As in Webb et al. (2018a), we consider only ultra-relativistic particles for which the particle speed v ∼ c in the fluid frame, where c is the speed of light, and the particle diffusion coefficient κ ∼ c²τ/3. We choose

$$\begin{eqnarray}&&g(\xi )=k(\xi -{\xi }_{0})\end{eqnarray} \tag{ 22 }$$

(other interesting cases with g(ξ₀) = 0 are possible). In the case (22) we choose ξ so that

$$\begin{eqnarray}&&\displaystyle \frac{g(\xi )}{{rd}\xi /{dr}}=\displaystyle \frac{k(\xi -{\xi }_{0})}{{rd}\xi /{dr}}=1,\end{eqnarray} \tag{ 23 }$$

where k is a positive constant. For these choices, we obtain

$$\begin{eqnarray}&&\tau ={\tau }_{0}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }.\end{eqnarray} \tag{ 24 }$$

Integrating (23) leads to the class of shear flows, for which

$$\begin{eqnarray}&&\xi (r)={\xi }_{0}-({\xi }_{0}-{\xi }_{2}){\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k},\quad \beta (r)=\tanh (\xi ).\end{eqnarray} \tag{ 25 }$$

If we choose k > 1, then dξ/dr → 0 as r → 0, but if 0 < k < 1, then $d\xi /{dr}\to \infty$ as r → 0. For the case k = 1, dξ/dr → −ξ₀₂/r₂ as r → 0, where ξ₀₂ = ξ₀ − ξ₂. We choose the cases k ≥ 1, which ensures that dξ/dr is bounded as r → 0.

For the above choice of g(ξ) and ξ and τ(r, p), the transport Equation (21) reduces to the form:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {\eta }^{2}}+\displaystyle \frac{1}{\eta }\displaystyle \frac{\partial {f}_{0}}{\partial \eta }+\displaystyle \frac{1}{5}\left(\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {T}^{2}}+(3+\alpha )\displaystyle \frac{\partial {f}_{0}}{\partial T}\right)\\ \quad =\,\displaystyle \frac{3{N}_{0}\delta (T)\delta (\xi -{\xi }_{1})}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}g({\xi }_{1})},\end{array}\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}\eta (r) & = & {\xi }_{0}-\xi (r)={\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k},\quad {\xi }_{02}={\xi }_{0}-{\xi }_{2},\\ g({\xi }_{1}) & = & k({\xi }_{1}-{\xi }_{0}),\quad T=\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right).\end{array}\end{eqnarray} \tag{ 27 }$$

In the next section we obtain the Green's function solutions of (26) and (27), satisfying the boundary conditions (17) and (18) and corresponding to steady injection of particles with momentum p = p₀ at radius r = r₁.

3.1. The Green's Function

The Green's function solution of the steady-state transport Equation (26), for monoenergetic injection of particles with momentum p = p₀ at radius r = r₁, with free escape of particles at the outer radius of the jet at radius r = r₂, and with no sources of particles on the axis of the jet at r = 0, may be expressed in the form

$$\begin{eqnarray}&&{f}_{0}=D{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{-a}\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{J}_{0}({\lambda }_{n}\eta ){J}_{0}({\lambda }_{n}{\eta }_{1})\exp \left(-{\chi }_{n}| T| \right)}{{\left[{J}_{1}({j}_{0,n})\right]}^{2}{\chi }_{n}},\end{eqnarray} \tag{ 28 }$$

where

$$\begin{eqnarray}&&D=\displaystyle \frac{15{N}_{0}}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k{\eta }_{2}^{2}},\quad {\eta }_{1}=\eta ({r}_{1}),\quad {\eta }_{2}=\eta ({r}_{2}).\end{eqnarray} \tag{ 29 }$$

In (28) J₀(z) and J₁(z) are ordinary Bessel functions of the first kind. The quantities j_0,n denote the positive zeros of the J₀(z)—that is,

$$\begin{eqnarray}&&{J}_{0}({j}_{0,n})=0,\quad n=1,2,3,\,\,\ldots .\end{eqnarray} \tag{ 30 }$$

The variables η and T in (28) are

$$\begin{eqnarray}&&\eta ={\xi }_{0}-\xi ={\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k},\quad T=\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right).\end{eqnarray} \tag{ 31 }$$

The other parameters in the solution (28) are

$$\begin{eqnarray}\begin{array}{rcl}{\xi }_{02} & = & {\xi }_{0}-{\xi }_{2}=\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{1+{\beta }_{02}}{1-{\beta }_{02}}\right),\\ {\beta }_{02} & = & \displaystyle \frac{{\beta }_{0}-{\beta }_{2}}{1-{\beta }_{0}{\beta }_{2}},\quad a=\displaystyle \frac{3+\alpha }{2},\\ {\chi }_{n} & = & {\left({a}^{2}+5{\lambda }_{n}^{2}\right)}^{1/2},\quad {\lambda }_{n}=\displaystyle \frac{{j}_{0,n}}{{\xi }_{02}},\quad n=1,2,3,\,\ldots .\end{array}\end{eqnarray} \tag{ 32 }$$

A derivation of the Green's function solution (28) for f₀ is given in Appendix A.

The limit of large r₂

In the limit of ${r}_{2}\gg r,{r}_{1}$ we can write

$$\begin{eqnarray}\begin{array}{rcl}\eta & = & {\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k}={\rm{\Lambda }}{r}^{k},\quad {\rm{\Lambda }}={\xi }_{02}{r}_{2}^{-k}={r}_{1}^{-k},\\ \eta & = & {\left(\displaystyle \frac{r}{{r}_{1}}\right)}^{k},\quad {\xi }_{02}={\left(\displaystyle \frac{{r}_{2}}{{r}_{1}}\right)}^{k},\quad {\xi }_{02}\gg 1.\end{array}\end{eqnarray} \tag{ 33 }$$

In (33) we have set ${\xi }_{02}={\left({r}_{2}/{r}_{1}\right)}^{k}$, which can be taken as a definition of the value of r₁ in this particular instance. This is because we are interested in the limit ${\xi }_{02}\gg 1$ and ${r}_{2}/{r}_{1}\gg 1$, which corresponds to the eigenspectrum for λ_n becoming very closely spaced (i.e., the eigenspectrum for λ_n becomes continuous in this limit). In this limit,

$$\begin{eqnarray}&&{\lambda }_{n}=\displaystyle \frac{{j}_{0,n}}{{\xi }_{02}}\end{eqnarray} \tag{ 34 }$$

only gives a finite λ_n if ${j}_{0,n}\gg 1$ (because we assume ${\xi }_{02}\gg 1$). Effectively, both ξ₀₂ and r₂ are large in this limit. Using the asymptotic form of Bessel functions for large arguments,

$$\begin{eqnarray}&&{J}_{n}(z)\sim \sqrt{\displaystyle \frac{2}{\pi z}}\cos \left(z-\displaystyle \frac{n\pi }{2}-\displaystyle \frac{\pi }{4}\right),\quad | \arg \ (z)| \lt \pi ,\end{eqnarray} \tag{ 35 }$$

(Abramowitz & Stegun 1965, formula 9.2.1, p. 364) we obtain the approximation:

$$\begin{eqnarray}&&{\left[{J}_{1}({\lambda }_{n}{\xi }_{02})\right]}^{2}\sim \displaystyle \frac{2}{\pi {j}_{0n}}=\displaystyle \frac{2}{\pi {\lambda }_{n}{\xi }_{02}}.\end{eqnarray} \tag{ 36 }$$

In the above limit, the solution for f₀ in (28) can be reduced to an infinite integral over λ in which λ_n → λ and

$$\begin{eqnarray}\begin{array}{rcl}{j}_{0,n} & \sim & \left(n-\displaystyle \frac{1}{2}\right)\pi +\displaystyle \frac{\pi }{4},\quad (n\,\mathrm{integer}),\\ d\lambda & \sim & \displaystyle \frac{({j}_{0,n+1}-{j}_{0,n})}{{\xi }_{02}}=\displaystyle \frac{\pi }{{\xi }_{02}}.\end{array}\end{eqnarray} \tag{ 37 }$$

The solution (28) then reduces to the infinite integral:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0} & \sim & \displaystyle \frac{15{N}_{0}}{16{\pi }^{2}{p}_{0}^{3}{\tau }_{0}k}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{-a}{\displaystyle \int }_{0}^{\infty }d\lambda \ \displaystyle \frac{\lambda }{\sqrt{{a}^{2}+5{\lambda }^{2}}}\\ & & \times \,\exp \left(-\sqrt{{a}^{2}+5{\lambda }^{2}}| T| \right){J}_{0}\left(\lambda \eta \right){J}_{0}\left(\lambda {\eta }_{1}\right).\end{array}\end{eqnarray} \tag{ 38 }$$

Note in the approximation (36) and (37) we have only approximated the eigenvalue equation and have not assumed anything about the size of r. Also note that as r → 0, λη → 0 and the solution for f₀ is finite as r → 0. Using (33) we have η₁ = 1 and $\eta ={\left(r/{r}_{1}\right)}^{k}$. We have not been able to reduce the integral (38) to a simpler form. The main point about the above reduction is that the eigenfunction sum (28) for the solution reduces to an infinite integral in this limit.

3.2. Green's Function Properties

Figure 1 shows a plot of the fluid velocity β(r) versus r/r₂ using formula (25), for the fluid velocity of the jet:

$$\begin{eqnarray}&&\beta (r)=\tanh (\xi ),\quad \xi (r)={\xi }_{0}-{\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k},\end{eqnarray} \tag{ 39 }$$

for a range of values of the parameter k (k = 0.5, 1, 2, 5) for a jet with β₀ = 0.5, β₂ = 0. It shows that the velocity profile has a finite, non-zero radial derivative for k = 1 at r = 0 but du/dr = 0 at r = 0 for k > 1, with the profile becoming flatter near r = 0 with increasing k.

Zoom In Zoom Out Reset image size

Figures 2–7 show some of the characteristics of the Green's function solution (28). Figure 2 shows a plot of Log(f₀(r, p)/D) versus Log(p/p₀), where the normalization constant $D\,=15{N}_{0}/(8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k{\eta }_{2}^{2})$ and the injection rate N₀ is defined in (14). The figure is for a jet with β₀ = 0.5, β₂ = 0, and β₁ = 0.25. The solution (28) is a function of r via the rapidity ξ(r), which in turn depends on the parameter k. However, the characteristics of the solution written in terms of ξ(r) or β(r) do not depend explicitly on the parameter k.

Zoom In Zoom Out Reset image size

The locations r = 0, r = r₁ and r = r₂ correspond to the axis of the jet, the injection radius r₁, and the outer radius of the jet (r = r₂), respectively. The figure shows how the spectrum for f₀(r, p) varies for different observational radial positions for β(r). The figure is analogous to Figure 2 of Webb et al. (2018a), except that $\tau (r,p)={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$ is independent of r (a more complicated radial dependence of τ(r, p) was used in Webb et al. (2018a), where $\tau \to \infty$ as r → 0). The parameter α = 1 in Figure 2. The particle momentum spectrum at large momenta is a power law in p with ${f}_{0}\propto {p}^{-{\mu }_{\infty }}$ (in this particular case ${\mu }_{\infty }=12$). At low momenta (0 < p < p₀), ${f}_{0}\propto {p}^{{\mu }_{0}}$ as p → 0 (μ₀ = 8 for the parameters used in above example). At high energies the spectral index is given by the formula

$$\begin{eqnarray}&&{\mu }_{\infty }={\left[\displaystyle \frac{{\left(3+\alpha \right)}^{2}}{4}+5\displaystyle \frac{{j}_{0,1}^{2}}{{\xi }_{02}^{2}}\right]}^{1/2}+\displaystyle \frac{3+\alpha }{2},\end{eqnarray} \tag{ 40 }$$

which corresponds to the lowest eigenfunction (n = 1) in the series (28). Here j_0,1 ∼ 2.4048 is the first zero of the J₀(z) Bessel function (Abramowitz & Stegun 1965, Table 9.5, p. 409). For low momenta (0 < p < p₀), the power-law index μ₀ is given by

$$\begin{eqnarray}&&{\mu }_{0}={\mu }_{\infty }-(3+\alpha ).\end{eqnarray} \tag{ 41 }$$

The index μ₀ > 0. There is a net radial outflux of particles in the region 0 < r < r₁ at all momenta, except near p ≈ p₀, where the flux is inward. In the region r > r₁ there is a net outward diffusive flux of particles at all momenta (see Webb et al. 2018a). Figure 3 shows the dependence of the asymptotic spectral index ${\mu }_{\infty }$ on the fluid velocity β₀ on the jet axis, for a range of momentum dependencies of the scattering time τ (τ ∝ p^α, α = 0, 1, 2). At non-relativistic flow velocities of the jet,

$$\begin{eqnarray}&&{\mu }_{\infty }={\left[\displaystyle \frac{{\left(3+\alpha \right)}^{2}}{4}+5\displaystyle \frac{{j}_{0,1}^{2}}{{\xi }_{02}^{2}}\right]}^{1/2}+\displaystyle \frac{3+\alpha }{2}\approx \displaystyle \frac{\sqrt{5}{j}_{\mathrm{0,1}}}{{\beta }_{02}}.\end{eqnarray} \tag{ 42 }$$

Thus, ${\mu }_{\infty }\to \infty$ as β₀₂ → 0, implying that shear acceleration is not effective for non-relativistic laminar (non-turbulent) flows. In the opposite limit of ultra-relativistic jets (β₀₂ → 1 and ${\xi }_{02}\to \infty$) we obtain

$$\begin{eqnarray}&&{\mu }_{\infty }\to 3+\alpha \quad \mathrm{as}\quad {\beta }_{02}\to 1\end{eqnarray} \tag{ 43 }$$

(see also Webb et al. 2018a and similar results from Rieger & Duffy 2004a, 2004b, 2005). Thus shear acceleration is quite effective in relativistic jet shear flows. ${\mu }_{\infty }\to 4$ as β₀₂ → 1 for α = 1, which is comparable to the spectral index for particles accelerated by the first order Fermi mechanism at strong non-relativistic shocks with compression ratio r_c = 4 (e.g., Axford et al. 1977; Krymsky 1977; Bell 1978a, 1978b; Blandford & Ostriker 1978).

Zoom In Zoom Out Reset image size

Figure 4 shows the hardening of the spectrum of the accelerated particles with increasing jet velocity β₀ at the source radius (r = r₁) for a range of β₀ (β₀ = 0.3, 0.4, 0.6, 0.8, and 0.99). As the central jet velocity β₀ increases, both the low energy particles with p < p₀ and very-high-energy particles with p > p₀ increase due to the increased diffusive spreading and systematic momentum changes.

Zoom In Zoom Out Reset image size

Figure 5 shows Log(f₀(r, p)/D) on the axis of the jet for the case β = β₀ = 0.5 and α = 1 for a range of injection radii as parameterized by β₁ (β₁ = 0.001, 0.01, 0.2, and 0.5) and for β₂ = 0 and α = 1. The smaller the value of β₁, the smaller the value of f₀ in Figure 5. For smaller β₁, the source is located closer to the free escape boundary at r = r₂, where β₂ = 0, resulting in enhanced escape of particle through the free escape boundary at r = r₂.

Zoom In Zoom Out Reset image size

Figure 6 shows the variation of the spectrum f₀ at the source radius (β = β₁ = 0.25) for a range of diffusion coefficients with κ ∝ p^α (α = 0, 0.5, 1 and 2), for the Green's function (28) with β₂ = 0 and β₀ = 0.5. Figure 7 shows the effect of a back-flowing cocoon for the case β₀ = 0.5, β = β₁ = 0.25, for a range of β₂ (β₂ = 0, −0.1, −0.2) and α = 1. The spectrum hardens as the net shear (as measured by β₀₂ or β₀ − β₂) increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. The Monoenergetic Spectrum Solution

In this section we use Green's formula (Webb et al. 2018a) to solve boundary value problems for the steady-state cosmic ray transport equation:

$$\begin{eqnarray}&&{ \mathcal L }({f}_{0})=-\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r\kappa \displaystyle \frac{\partial {f}_{0}}{\partial r}\right)-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({\rm{\Gamma }}{p}^{4}\tau \displaystyle \frac{\partial {f}_{0}}{\partial p}\right)=Q(r,p).\end{eqnarray} \tag{ 44 }$$

The solution method uses the Green's functions (28) for particle transport derived in Section 3, with monoenergetic source term of the form

$$\begin{eqnarray}&&Q(r,p)=\displaystyle \frac{{N}_{0}\delta (p-{p}_{0})\delta (r-{r}_{1})}{8{\pi }^{2}{p}_{0}^{2}{r}_{1}},\end{eqnarray} \tag{ 45 }$$

and f₀ satisfies the boundary conditions (i) f₀(r₂, p) = 0 and (ii) rκ ∂f₀/∂r → 0 as r → 0.

In this section we obtain the solution of (44) with Q = 0, and for a monoenergetic spectrum of particles specified on the boundary surface at r = r₂ by using Green's formula. The boundary conditions for f₀ are

$$\begin{eqnarray}\begin{array}{rcl}r\kappa \displaystyle \frac{\partial {f}_{0}}{\partial r} & \to & 0\quad \mathrm{as}\quad r\to 0,\\ {f}_{0}(r,p) & = & \displaystyle \frac{{N}_{g}\delta (p-{p}_{0})}{4\pi {p}_{0}^{2}}\quad \mathrm{at}\quad r={r}_{2}\end{array}\end{eqnarray} \tag{ 46 }$$

(i.e., a monoenergetic spectrum of particles is specified on the outer boundary r = r₂ of the jet).

The boundary value problem (46) for the transport Equation (44) with Q = 0 corresponds to specifying a monoenergetic spectrum of particles with momentum $p={p}_{0}$ at cylindrical radius r = r₂ from the jet axis. The particles are accelerated and transported in the jet shear flow in the region 0 < r < r₂, and eventually free escape from the acceleration region 0 < r < r₂ through the boundary r = r₂ (note f₀(r₂, p) = 0 for $p\ne {p}_{0}$). From Webb et al. (2018a), Equation (83) et seq., Green's formula for the transport Equation (44) gives the solution of the boundary value problem for f₀(r, p) in the region:

$$\begin{eqnarray}&&{ \mathcal R }=\left\{(r,p):\quad 0\lt r\lt {r}_{2},\quad 0\lt p\lt \infty \right\},\end{eqnarray} \tag{ 47 }$$

by the formula

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(r,p) & = & {\displaystyle \int }_{0}^{{r}_{2}}r^{\prime} {dr}^{\prime} \ {\displaystyle \int }_{0}^{\infty }{p}^{{\prime} 2}{dp}^{\prime} {\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)Q(r^{\prime} ,p^{\prime} )\\ & & +{\displaystyle \int }_{0}^{\infty }{p}^{{\prime} 2}{dp}^{\prime} {\left[r^{\prime} \kappa ^{\prime} \left({\rm{\Gamma }}^{\prime} \displaystyle \frac{\partial {f}_{0}^{{\prime} }}{\partial r^{\prime} }-{f}_{0}^{{\prime} }\displaystyle \frac{\partial {\rm{\Gamma }}^{\prime} }{\partial r^{\prime} }\right)\right]}_{r^{\prime} =0}^{r^{\prime} ={r}_{2}}\\ & & +{\displaystyle \int }_{0}^{{r}_{2}}r^{\prime} {dr}^{\prime} {\left[{\rm{\Gamma }}^{\prime} {p}^{{\prime} 4}\tau ^{\prime} \left({\rm{\Gamma }}^{\prime} \displaystyle \frac{\partial {f}_{0}^{{\prime} }}{\partial p^{\prime} }-{f}_{0}^{{\prime} }\displaystyle \frac{\partial {\rm{\Gamma }}^{\prime} }{\partial p^{\prime} }\right)\right]}_{p^{\prime} =0}^{p^{\prime} =\infty },\end{array}\end{eqnarray} \tag{ 48 }$$

where ${\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)$ is the adjoint equation Green's function satisfying the transport equation:

$$\begin{eqnarray}&&{{ \mathcal L }}^{\dagger ^{\prime} }\left[{\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)\right]=\displaystyle \frac{\delta (r^{\prime} -r)\delta (p^{\prime} -p)}{{{rp}}^{2}},\quad {{ \mathcal L }}^{\dagger }\equiv { \mathcal L },\end{eqnarray} \tag{ 49 }$$

and Γ satisfies the boundary conditions:

$$\begin{eqnarray}\begin{array}{rcl} & & r^{\prime} \kappa (r^{\prime} ,p^{\prime} )\displaystyle \frac{\partial {\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)}{\partial r^{\prime} }\to 0\quad \mathrm{as}\quad r^{\prime} \to 0,\\ & & {\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)=0\quad \mathrm{at}\quad r^{\prime} ={r}_{2}.\end{array}\end{eqnarray} \tag{ 50 }$$

Note that the transport Equation (44) is self-adjoint (i.e., ${{ \mathcal L }}^{\dagger }={ \mathcal L }$). Equation (48) is Green's formula for the transport Equation (44), which gives the solution for f₀(r, p) for (r, p) in ${ \mathcal R }$, in terms of the adjoint Green's function, the source of particles in ${ \mathcal R }$, and in terms of fluxes through the spatial and momentum boundaries.

Following Webb et al. (2018a), the solution of the boundary value problem (46) for ${f}_{0}(r,p)$ for a monoenergetic spectrum of particles with momentum p = p₀ specified at radius r = r₂ can be obtained by inserting the appropriate boundary conditions in the Green's formula (48) to obtain

$$\begin{eqnarray}&&{f}_{0}(r,p)=-\displaystyle \frac{{N}_{g}}{4\pi }{\left[r^{\prime} \kappa ^{\prime} (r^{\prime} ,p^{\prime} )\displaystyle \frac{\partial {\rm{\Gamma }}^{\prime} }{\partial r^{\prime} }\right]}_{r^{\prime} ={r}_{2},p^{\prime} ={p}_{0}}\end{eqnarray} \tag{ 51 }$$

as the required solution for f₀(r, p). Here ${\rm{\Gamma }}^{\prime} \equiv {\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)$ in (51). In the derivation of (51), the contributions to the particle momentum fluxes through the boundaries at p' = 0 and $p^{\prime} \to \infty$ are set equal to zero.

To solve the boundary value problem (46), we use the Green's function (28) to determine the adjoint Green's function ${\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p)$. This is achieved by using the map,

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{0}}{8{\pi }^{2}}\to 1,\quad r\to r^{\prime} ,\quad p\to p^{\prime} ,\quad {r}_{1}\to r,\quad {p}_{0}\to p,\end{eqnarray} \tag{ 52 }$$

in the Green's function solution (28). These substitutions in (28) give the adjoint Green's function:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Gamma }}(r^{\prime} ,p^{\prime} ;r,p) & = & D\exp \left[-\left(\displaystyle \frac{3+\alpha }{2}\right)T(p^{\prime} ,p)\right]\displaystyle \sum _{n=1}^{\infty }\\ & & \times \displaystyle \frac{{J}_{0}({\lambda }_{n}\eta ){J}_{0}({\lambda }_{n}\eta ^{\prime} )}{{[{J}_{1}({j}_{0,n})]}^{2}{\chi }_{n}}\exp \left[-{\chi }_{n}| T(p^{\prime} ,p)| \right],\end{array}\end{eqnarray} \tag{ 53 }$$

where

$$\begin{eqnarray}&&T(p^{\prime} ,p)=\mathrm{ln}\left(\displaystyle \frac{p^{\prime} }{p}\right).\end{eqnarray} \tag{ 54 }$$

In (53)

$$\begin{eqnarray}&&\eta ={\eta }_{2}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k},\quad \eta ^{\prime} ={\eta }_{2}{\left(\displaystyle \frac{r^{\prime} }{{r}_{2}}\right)}^{k},\quad {\eta }_{2}={\xi }_{02},\end{eqnarray} \tag{ 55 }$$

and

$$\begin{eqnarray}&&D=\displaystyle \frac{15}{{p}^{3}{c}^{2}{\tau }_{0}k{\eta }_{2}^{2}}.\end{eqnarray} \tag{ 56 }$$

Note that we set ${N}_{0}/(8{\pi }^{2})\to 1$ and p₀ → p in (29) to obtain the expression for D in (56). The diffusion coefficient $\kappa (r^{\prime} ,p^{\prime} )$ in the map (52) is given by

$$\begin{eqnarray}&&\kappa (r^{\prime} ,p^{\prime} )={\kappa }_{0}{\left(\displaystyle \frac{p^{\prime} }{p}\right)}^{\alpha }\equiv {\kappa }_{0}\exp \left[\alpha T(p^{\prime} ,p)\right].\end{eqnarray} \tag{ 57 }$$

Using the formulae (52)–(57) in (51) we obtain the monoenergetic spectrum solution:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(r,p) & = & {D}_{1}\exp \left[-\left(\displaystyle \frac{3+\alpha }{2}\right)T(p,{p}_{0})\right]\displaystyle \sum _{n=1}^{\infty }\\ & & \times \,\displaystyle \frac{{j}_{0,n}{J}_{0}({\lambda }_{n}\eta )}{{J}_{1}({j}_{0,n})}\displaystyle \frac{\exp \left(-{\chi }_{n}| T(p,{p}_{0})| \right)}{{\chi }_{n}},\end{array}\end{eqnarray} \tag{ 58 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{D}_{1} & = & \displaystyle \frac{5{N}_{g}}{4\pi {\eta }_{2}^{2}{p}_{0}^{3}},\quad {\lambda }_{n}=\displaystyle \frac{{j}_{0,n}}{{\xi }_{02}},\\ {\chi }_{n} & = & {\left({a}^{2}+5{\lambda }_{n}^{2}\right)}^{1/2},\quad a=\displaystyle \frac{3+\alpha }{2},\quad n=1,2,3,\,\ldots ,\\ T(p,{p}_{0}) & = & \mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right),\quad \eta ={\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k}.\end{array}\end{eqnarray} \tag{ 59 }$$

In the derivation of (58) we used ${J}_{0}^{{\prime} }(z)=-{J}_{1}(z)$ (Abramowitz & Stegun 1965, formula 9.1.28, p. 361) in calculating the derivative of Γ with respect to r'.

4.2. Solution Characteristics

Figures 8 and 9 show some of the main characteristics of the monoenergetic spectrum solution (58) in which monoenergetic particles enter the shear flow region 0 < r < r₂ with momentum p = p₀ and eventually escape from the region through the free escape boundary r = r₂ after being either decelerated or accelerated by their interaction with the shear flow.

Zoom In Zoom Out Reset image size

Figure 8 shows plots of the solution (58) for f₀(r, p)/D₁ versus Log(p/p₀) for the case β₀ = 0.99, β₂ = 0, and α = 1. The momentum spectrum of the particles f₀(r, p)/D₁ are shown at different radial locations parameterized by the flow velocity β(r), which is a monotonic decreasing function of r. The spectra are shown at the locations β = 0.01, 0.1, 0.5, and 0.9. Near the edge of the jet β = 0.01 the spectrum for f₀(r, p) is nearly monoenergetic. Harder momentum spectra or particles with momentum p > p₀ are obtained as one moves closer to the jet axis as r decreases and β(r) increases. Both more accelerated particles and more decelerated particles are obtained as r decreases, because these particles have interacted both with the inner and outer parts of the flow. The particle momentum spectra at large p/p₀ has the power-law form ${f}_{0}\propto {p}^{-{\mu }_{\infty }}$, where ${\mu }_{\infty }$ corresponds to the first term in the eigen-series solution (58)—that is, ${\mu }_{\infty }$ is given by formula (40). At low momenta, the distribution function is a positive power law in momentum p with f₀ ∝ p^μ₀ at p ≪ p₀ (the expression for μ₀ is given in (43)).

Figure 9 shows a similar plot of the distribution function f₀(r, p) to that in Figure 8, but for the case of a jet with central velocity β₀ = 0.8 and β₂ = 0 at the outer boundary of the jet at radius r = r₂. The parameter α = 1. The figure shows the variation in the spectra for f₀(r, p) with ${f}_{0}\propto {p}^{-{\mu }_{\infty }}$ at p ≫ p₀, and f₀ ∝ p^μ₀ at momenta p ≪ p₀. The parameter β = 0.01, 0.1, and 0.5 in Figure 9. The spectrum approaches a monoenergetic spectrum as r → r₂, where β₂ = 0 at the outer boundary of the acceleration region.

Zoom In Zoom Out Reset image size

4.3. General Spectrum Solution at Jet Boundary at $r={r}_{2}$

The general solution for f₀(r, p) in which f₀(r₂, p) is specified at r = r₂ is given from (48) by the integral:

$$\begin{eqnarray}&&{f}_{0}(r,p)=-{\int }_{0}^{\infty }{dp}^{\prime} {p}^{{\prime} 2}{f}_{0}({r}_{2},p^{\prime} ){\left(r^{\prime} \kappa ^{\prime} \displaystyle \frac{\partial {\rm{\Gamma }}^{\prime} }{\partial r^{\prime} }\right)}_{r^{\prime} ={r}_{2}}.\end{eqnarray} \tag{ 60 }$$

Using ${\rm{\Gamma }}^{\prime} (r^{\prime} ,p^{\prime} ;r,p)$ from (53) in (60) gives

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(r,p) & = & \displaystyle \frac{5}{{\eta }_{2}^{2}{p}^{3}}{\displaystyle \int }_{0}^{\infty }{dp}^{\prime} {p}^{{\prime} 2}\exp \left[(\alpha -a)T(p^{\prime} ,p)\right]{f}_{0}({r}_{2},p^{\prime} )\\ & & \times \left(\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{J}_{0}({\lambda }_{n}\eta ){j}_{0,n}}{{J}_{1}({j}_{0,n})}\displaystyle \frac{\exp [-{\chi }_{n}| T(p^{\prime} ,p)| ]}{{\chi }_{n}}\right)\end{array}\end{eqnarray} \tag{ 61 }$$

as the solution for f₀(r, p) in the region 0 < r < r₂, in which f₀(r₂, p) is specified on the boundary at r = r₂. The solution (61) can also be written in the form:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(r,p) & = & \displaystyle \frac{5}{{\eta }_{2}^{2}{p}^{a}}\left[{\displaystyle \int }_{0}^{p}{dp}^{\prime} {p}^{{\prime} a-1}{f}_{0}({r}_{2},p^{\prime} )\right.\\ & & \times \,\left(\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{J}_{0}({\lambda }_{n}\eta ){j}_{0,n}}{{J}_{1}({j}_{0,n}){\chi }_{n}}{\left(\displaystyle \frac{p^{\prime} }{p}\right)}^{{\chi }_{n}}\right)\\ & & +{\displaystyle \int }_{p}^{\infty }{dp}^{\prime} {p}^{{\prime} a-1}{f}_{0}({r}_{2},p^{\prime} )\\ & & \times \,\left.\left(\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{J}_{0}({\lambda }_{n}\eta ){j}_{0,n}}{{J}_{1}({j}_{0,n}){\chi }_{n}}{\left(\displaystyle \frac{p^{\prime} }{p}\right)}^{-{\chi }_{n}}\right)\right].\end{array}\end{eqnarray} \tag{ 62 }$$

In general, it is difficult to include particle energy losses and gains and particle escape all together in an analytic model of particle acceleration in relativistic radio jets, without invoking some approximations and simplifications. Simplified diffusion equation models for particle acceleration in radio-jet shear flows were developed by Berezhko (1981) and Rieger & Duffy (2004a, 2004b, 2006) and Ostrowski (2002), Stawarz & Ostrowski (2002). We generalize this model to include the effects of cosmic ray inertia by using the telegrapher equation model governed by the transport equation,

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{\widetilde{\gamma }}^{2}\tau \displaystyle \frac{{{\rm{\partial }}}^{2}{f}_{0}}{{\rm{\partial }}{t}^{2}}+\gamma \displaystyle \frac{{\rm{\partial }}{f}_{0}}{{\rm{\partial }}t}-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}p}\left({p}^{2}{D}_{{pp}}\displaystyle \frac{{\rm{\partial }}{f}_{0}}{{\rm{\partial }}p}\right)\\ \,\,-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}p}\left({D}_{s}{p}^{4}{f}_{0}\right)+\displaystyle \frac{{f}_{0}}{{\tau }_{{\rm{esc}}}}=Q(p,t),\end{array}\end{array}\end{eqnarray} \tag{ 63 }$$

in which γ is the Lorentz factor of the flow, and f₀/τ_esc represents particle losses from the acceleration region (this model is analogous to the leaky box model used to describe galactic cosmic rays). Here Q(p, t) represents particle sources. In (63) we have assumed that the net effect of averaging the space dependent transport equation, including the telegrapher second-order time derivative term results in a transport equation with space averaged transport coefficients. The exact form of the space averaged transport coefficients must depend on the flow velocity profile of the jet about the jet axis. We do not attempt to carry out this averaging in the present paper. This makes our analysis qualitative rather than quantitative, and more work needs to be done to obtain more quantitative results.

5.1. Telegrapher Equation Effects

The first term in (63) represents a telegrapher equation term, which takes into account the finite inertia of the cosmic ray gas (e.g., Webb et al. 2018a), where

$$\begin{eqnarray}&&{\tilde{\gamma }}^{2}={\gamma }^{2}(1-{\beta }_{p}^{2}{\beta }^{2}/3).\end{eqnarray} \tag{ 64 }$$

Here β = u/c and β_p = v/c give the fluid speed u, and the particle speed v normalized to the speed of light. The second time derivative term in (63) limits the upper momentum p_max that takes into account the finite particle speed and inertia. In the diffusion approximation the second time derivative term is neglected. Retention of this term is important at high energies, because the diffusion approximation allows particles to attain an infinite energy in an infinitesimal time. We discuss the implications of this term in Appendix D.

In Appendix D we show that the characteristic manifold S(t, p) = const. for the second-order partial differential Equation (63) satisfies the characteristic equation

$$\begin{eqnarray}&&G\equiv {\tilde{\gamma }}^{2}\tau {S}_{t}^{2}-{D}_{{pp}}{S}_{p}^{2}=0.\end{eqnarray} \tag{ 65 }$$

The momentum space diffusion coefficient D_pp can be split up into two terms:

$$\begin{eqnarray}&&{D}_{{pp}}={D}_{{pp}}^{\mathrm{sh}}+{D}_{{pp}}^{\mathrm{sf}},\end{eqnarray} \tag{ 66 }$$

where

$$\begin{eqnarray}&&{D}_{{pp}}^{\mathrm{sh}}={\rm{\Gamma }}{p}^{2}\tau ,\quad {D}_{{pp}}^{\mathrm{sf}}=\displaystyle \frac{{p}^{2}{V}_{{\rm{A}}}^{2}}{9\kappa },\quad \kappa =\displaystyle \frac{{v}^{2}\tau }{3},\quad (v\approx c)\end{eqnarray} \tag{ 67 }$$

describe stochastic particle acceleration due to cosmic ray viscosity (D_pp^sh) and second-order Fermi acceleration, owing to Alfvénic turbulence (D_pp^sf).

From Appendix D, the bi-characteristics of (63) (i.e., the characteristics of (65)) are given by the ordinary differential equation:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dp}}{{dt}} & = & \sigma \displaystyle \frac{p}{\tilde{\gamma }}{\left[{\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}{\tau }_{0}^{2}}{\left(\bar{p}\right)}^{-2\alpha }\right]}^{1/2}\quad \mathrm{where}\\ \bar{p} & = & \displaystyle \frac{p}{{p}_{0}}\quad \mathrm{and}\quad \sigma =\pm 1.\end{array}\end{eqnarray} \tag{ 68 }$$

The maximum particle momentum p_max(t) allowed in the telegrapher Equation (63) corresponds to the σ = +1 possibility in (68). In (68) β_A = V_A/c is the relativistic Alfvén speed divided by the speed of light. The case σ = −1 in (68) corresponds to the minimum momentum implied by the telegrapher Equation (63).

In the formulae below we take β_p = 1 corresponding to the case of fully relativistic particles only. In the more general case where the particles are not necessarily relativistic, Equation (68) can be integrated numerically if required.

For the special case of no second-order Fermi acceleration due to Alfvénic turbulence (β_A = 0, ${\rm{\Gamma }}\ne 0$), the σ = +1 version of (68) reduces to

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dt}}=\displaystyle \frac{\sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }}p.\end{eqnarray} \tag{ 69 }$$

Integrating (69) with initial condition p = p₀ at time t = 0 gives the solution

$$\begin{eqnarray}&&p={p}_{0}\exp \left(\displaystyle \frac{\sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right),\end{eqnarray} \tag{ 70 }$$

for the maximum particle momentum due to pure shear acceleration alone. On the other hand, for the case of no shear (${\rm{\Gamma }}=0$) and inclusion of second-order Fermi acceleration (${\beta }_{{\rm{A}}}\ne 0$), (68) with σ = +1, reduces to the equation

$$\begin{eqnarray}&&\displaystyle \frac{d\bar{p}}{{dt}}=\displaystyle \frac{{\beta }_{{\rm{A}}}{\left(\bar{p}\right)}^{1-\alpha }}{\sqrt{3}\tilde{\gamma }{\tau }_{0}}\quad \mathrm{where}\quad \bar{p}=\displaystyle \frac{p}{{p}_{0}}.\end{eqnarray} \tag{ 71 }$$

Integration of (71) with $\bar{p}=1$ at time t = 0 (i.e., p = p₀ at t = 0) gives the solution

$$\begin{eqnarray}&&p(t)={p}_{0}{\left[1+\displaystyle \frac{\alpha {\beta }_{{\rm{A}}}}{\sqrt{3}\tilde{\gamma }}\displaystyle \frac{t}{{\tau }_{0}}\right]}^{1/\alpha },\end{eqnarray} \tag{ 72 }$$

for the maximum particle momentum implied by the telegrapher Equation (63) due to second-order Fermi acceleration alone.

The solution of (68) with σ = +1 for the general case (${\beta }_{{\rm{A}}}\ne 0$ and ${\rm{\Gamma }}\ne 0$) is given by (198)—that is,

$$\begin{eqnarray}&&p(t)={p}_{0}{\left\{\cosh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)+{\left({A}^{2}+1\right)}^{1/2}\sinh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)\right\}}^{1/\alpha },\end{eqnarray} \tag{ 73 }$$

where

$$\begin{eqnarray}&&A=\displaystyle \frac{{\beta }_{{\rm{A}}}}{\sqrt{3{\rm{\Gamma }}}{\tau }_{0}}.\end{eqnarray} \tag{ 74 }$$

Note that (72)–(74) apply to the case $\alpha \ne 0$ (see Appendix D for the case α = 0).

5.1.1. Lower Limits to Acceleration Timescales

It is interesting to invert the formulae (70)–(73) to give t as a function of p and p₀ (i.e., t = t(p, p₀)), which intuitively should give a lower time estimate for particles to be accelerated from momentum p₀ up to momentum p (see, e.g., Drury 1983; Malkov & Drury 2001) for similar formulae for particles to be accelerated from momentum p₀ up to momentum p in a mean time t by the first order Fermi mechanism at astrophysical shocks.

Inversion of (70) gives the formula

$$\begin{eqnarray}&&t=\displaystyle \frac{\tilde{\gamma }}{\sqrt{{\rm{\Gamma }}}}\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right)\end{eqnarray} \tag{ 75 }$$

as a lower limit to the time for particles to be accelerated from momentum p₀ up to momentum p by the viscous shear acceleration mechanism. Similarly we obtain the estimate

$$\begin{eqnarray}&&t={\tau }_{0}\displaystyle \frac{\sqrt{3}\gamma }{\alpha {\beta }_{{\rm{A}}}}\left[{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }-1\right]\end{eqnarray} \tag{ 76 }$$

for the lower limit for time for particles to be accelerated from p = p₀ up to momentum p solely due to second-order Fermi acceleration. More generally, inversion of (73) gives the formula

$$\begin{eqnarray}&&t(p,{p}_{0})=\displaystyle \frac{\tilde{\gamma }}{\sqrt{{\rm{\Gamma }}}\alpha }\mathrm{ln}\left[\displaystyle \frac{{\bar{p}}^{\alpha }+{\left({\bar{p}}^{2\alpha }+{A}^{2}\right)}^{1/2}}{1+{\left({A}^{2}+1\right)}^{1/2}}\right]\end{eqnarray} \tag{ 77 }$$

for the least time for particles to be accelerated from $p={p}_{0}$ up to momentum p due to both shear acceleration (${\rm{\Gamma }}\ne 0$) and due to second-order Fermi acceleration (${\beta }_{{\rm{A}}}\ne 0$).

As an application of these ideas, consider the timescale (75) for particles to be accelerated from p₀ up to p in a radio-jet shear flow. Use (64) for $\tilde{\gamma }$ and note that

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{\gamma }^{4}}{15}{\left(\displaystyle \frac{{du}}{{dr}}\right)}^{2}\sim \displaystyle \frac{{\gamma }_{0}^{4}}{15}{\left(\displaystyle \frac{{u}_{0}}{L}\right)}^{2},\end{eqnarray} \tag{ 78 }$$

where u₀ is the jet's central velocity and L is the width of the jet. From (75) we obtain

$$\begin{eqnarray}&&t(p,{p}_{0})\sim \displaystyle \frac{\sqrt{15}\tilde{\gamma }}{{\gamma }^{2}}\displaystyle \frac{L}{\beta c}\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right).\end{eqnarray} \tag{ 79 }$$

Noting that $\tilde{\gamma }=\gamma \sqrt{1-{\beta }^{2}/3}$, we obtain the formula

$$\begin{eqnarray}&&t(p,{p}_{0})\sim 1.2\times {10}^{4}\displaystyle \frac{\sqrt{1-{\beta }^{2}/3}}{\gamma \beta }{L}_{\mathrm{kpc}}\mathrm{ln}\left(\displaystyle \frac{p}{{p}_{0}}\right)\quad \mathrm{yr},\end{eqnarray} \tag{ 80 }$$

where L_kpc is the jet width in kilo-parsecs. Here we have used $1\ \mathrm{pc}\sim 3\times {10}^{18}\,\mathrm{cm}$ and 1 year = 3.154 × 10⁷ sections. If we consider p₀c = 10⁹ eV and pc = 10¹⁸ eV (i.e., p₀ = 1 GeV/c and p = 1 EeV/c, where p₀ corresponds roughly to the maximum in the cosmic ray intensity spectrum j_T = p²f₀, and p corresponds to an energy of the order of the ankle), then the formula (80) gives

$$\begin{eqnarray}&&t(p,{p}_{0})\sim 2.5\times {10}^{5}\displaystyle \frac{\sqrt{1-{\beta }^{2}/3}}{\gamma \beta }{L}_{\mathrm{kpc}}\quad \mathrm{yr},\end{eqnarray} \tag{ 81 }$$

where β = u/c is the jet speed. Note in (78)–(81) that these formulae depend on our estimate that du/dr ∼ u₀/L. It is interesting to compare this timescale with the light-crossing timescale t_lc across the jet, namely

$$\begin{eqnarray}&&{t}_{{lc}}=\displaystyle \frac{L}{c}=3.2\times {10}^{3}{L}_{\mathrm{kpc}}\quad \mathrm{yr}.\end{eqnarray} \tag{ 82 }$$

Thus the acceleration timescale (81) is about 100 light-crossing timescales (we assume $\sqrt{1-{\beta }^{2}/3}/(\gamma \beta )\sim 1$). Clearly for non-relativistic jets (e.g., β = 0.01) the acceleration time is much larger than 100 t_lc. The above estimates of the acceleration timescale in (80) and (81) can clearly be much greater than these estimates if particle energy losses are taken into account (especially synchrotron losses for electrons, which imply that the acceleration time t(p, p₀) will be much larger in that case).

5.2. Particle Momentum Losses and Gains

In this section we discuss the particle energy losses due to synchrotron losses as well as the particle energy gains due to shear acceleration. Losses due to the inverse Compton effect have a similar forms to that for synchrotron losses, but the detailed physics is different for inverse Compton losses. The transport Equation (63) divided by the Lorentz factor of the flow (γ) gives the time evolution equation for f₀ in the fixed observer's frame, in which u is the velocity of the jet along the axis of the jet. In this frame the transport Equation (63) becomes

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{{\tilde{\gamma }}^{2}\tau }{\gamma }\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {t}^{2}}+\displaystyle \frac{\partial {f}_{0}}{\partial t}-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{2}\displaystyle \frac{{D}_{{pp}}}{\gamma }\displaystyle \frac{\partial {f}_{0}}{\partial p}\right)\\ & & \quad -\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left(\displaystyle \frac{{D}_{s}}{\gamma }{p}^{4}{f}_{0}\right)+\displaystyle \frac{{f}_{0}}{\gamma {\tau }_{\mathrm{esc}}}=\displaystyle \frac{Q(p,t)}{\gamma }.\end{array}\end{eqnarray} \tag{ 83 }$$

The synchrotron losses in (83) are described by D_s and by the systematic particle momentum loss term:

$$\begin{eqnarray}&&{\left\langle \displaystyle \frac{{dp}}{{dt}}\right\rangle }_{s}=-\displaystyle \frac{p}{\gamma {t}_{\mathrm{synch}}}=-\displaystyle \frac{{D}_{s}{p}^{2}}{\gamma },\end{eqnarray} \tag{ 84 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{D}_{s} & = & \displaystyle \frac{1}{{m}_{0}c{\tau }_{r}},\quad {\tau }_{r}=\displaystyle \frac{6\pi {m}_{0}c}{{\sigma }_{T}{B}^{2}},\\ {\sigma }_{T} & = & \displaystyle \frac{8\pi {r}_{0}^{2}}{3},\quad {r}_{0}=\displaystyle \frac{{e}^{2}}{{m}_{0}{c}^{2}},\quad {t}_{\mathrm{synch}}=\displaystyle \frac{{\tau }_{r}}{{\gamma }_{p}{\beta }_{p}}.\end{array}\end{eqnarray} \tag{ 85 }$$

Here m₀ denotes the particle rest mass and ${\gamma }_{p}={\left(1-{\beta }_{p}^{2}\right)}^{-1/2}$ is the particle Lorentz factor; β_p = v/c; σ_T is the Thomson cross-section (m₀ = m_e for electrons and m₀ = m_p for protons). Formulas (85) are in cgs units. In (83)–(85) the transport equation was originally obtained using the fluid frame or the scattering frame to describe the particle momentum. Thus the magnetic field B in (83) refers to the value of B in the fluid frame. The form of the field in the jet frame can be obtained by using the Lorentz transformations between frames. If there is no electric field in the fluid frame, then ${B}_{\parallel }={B}_{\parallel }^{{\prime} }$ and ${B}_{\perp }=\gamma {B}_{\perp }^{{\prime} }$, where the subscripts ∥ and ⊥ refer to components parallel and perpendicular to the fluid velocity u and the primed superscript 'refers to the fluid frame, or scattering frame (e.g., Misner et al. 1973).

Berezhko (1981) and Rieger & Duffy (2006) derived the Green's function solution of the transport Equation (63) with source terms Qδ (t)δ(p − p₀) corresponding to impulsive monoenergetic injection of particles with momentum p = p₀ into a jet shear flow, with no synchrotron losses (D_s = 0), and with no particle escape (${\tau }_{\mathrm{esc}}\to \infty$). They used the solution to illustrate particle acceleration by cosmic ray viscosity in radio-jet shear flows. The mean systematic acceleration rate $\langle {\rm{\Delta }}p/{\rm{\Delta }}t\rangle$ is given by

$$\begin{eqnarray}&&\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle =\displaystyle \frac{1}{\gamma {p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{D}_{{pp}}\right)\quad \mathrm{and}\quad \left\langle \displaystyle \frac{{\left({\rm{\Delta }}p\right)}^{2}}{2{\rm{\Delta }}t}\right\rangle =\displaystyle \frac{{D}_{{pp}}}{\gamma },\end{eqnarray} \tag{ 86 }$$

where D_pp is the momentum dispersion coefficient. Using (66) and (67), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle }_{\mathrm{sh}} & = & \displaystyle \frac{{\rm{\Gamma }}p}{\gamma }\left(4\tau +p\displaystyle \frac{\partial \tau }{\partial p}\right)=\displaystyle \frac{{\rm{\Gamma }}{\tau }_{0}{p}_{0}}{\gamma }(4+\alpha ){\left(\bar{p}\right)}^{\alpha +1},\\ {\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle }_{{sf}} & = & \displaystyle \frac{(4-\alpha )}{3\gamma }\left(\displaystyle \frac{{p}_{0}}{{\tau }_{0}}\right)\left(\displaystyle \frac{{V}_{{\rm{A}}}^{2}}{{c}^{2}}\right){\left(\bar{p}\right)}^{1-\alpha }\quad \mathrm{where}\quad \bar{p}=\displaystyle \frac{p}{{p}_{0}},\end{array}\end{eqnarray} \tag{ 87 }$$

as the systematic acceleration rates due to shear acceleration and second-order Fermi acceleration, respectively.

Consider the case where the dominant acceleration mechanism is due to cosmic ray viscosity and shear in the background plasma flow, and in which synchrotron losses is the main energy loss mechanism. In this case there will be a balance between particle energization and synchrotron losses at the momentum p = p_* if

$$\begin{eqnarray}&&\langle \dot{p}{\rangle }_{\mathrm{sh}}+\langle \dot{p}{\rangle }_{s}=0,\end{eqnarray} \tag{ 88 }$$

at $p={p}_{* }$. Using (84) and (87) in (88), the balance condition (88) reduces to

$$\begin{eqnarray}&&(4+\alpha )\tau p{\rm{\Gamma }}-{D}_{s}{p}^{2}=0.\end{eqnarray} \tag{ 89 }$$

Equation (89) can also be written as

$$\begin{eqnarray}&&(4+\alpha ){\tau }_{0}{p}_{0}{\rm{\Gamma }}{\left(\bar{p}\right)}^{\alpha +1}-{D}_{s}{\left(\bar{p}\right)}^{2}{p}_{0}^{2}=0.\end{eqnarray} \tag{ 90 }$$

For $\alpha \ne 1$, (90) has the solution

$$\begin{eqnarray}&&{p}_{\mathrm{sh}}^{* }\,=\,{p}_{0}{\left(\displaystyle \frac{(4+\alpha ){\rm{\Gamma }}{\tau }_{0}}{{D}_{s}{p}_{0}}\right)}^{1/(1-\alpha )}.\end{eqnarray} \tag{ 91 }$$

In the case 0 < α < 1, synchrotron losses dominate over particle acceleration by shear for $p\gt {p}_{\mathrm{sh}}^{* }$, but shear acceleration dominates synchrotron losses for $p\lt {p}_{\mathrm{sh}}^{* }$. The expression for ${p}_{\mathrm{sh}}^{* }$ in (91) can also be written in the form

$$\begin{eqnarray}&&{p}_{\mathrm{sh}}^{* }\,=\,{p}_{0}{\left[\displaystyle \frac{(4+\alpha ){\gamma }^{4}}{15}\left(\displaystyle \frac{{\tau }_{0}{\tau }_{r}}{{\tau }_{\mathrm{sh}}^{2}}\right)\displaystyle \frac{{m}_{0}c}{{p}_{0}}\right]}^{1/(1-\alpha )},\end{eqnarray} \tag{ 92 }$$

where

$$\begin{eqnarray}&&{\tau }_{\mathrm{sh}}=\displaystyle \frac{1}{{du}/{dr}},\end{eqnarray} \tag{ 93 }$$

which emphasizes the role of the scattering timescale τ₀, the radiation timescale τ_r, and the shear flow timescale τ_sh.

However, if α > 1, synchrotron losses are dominant below p = p_* and shear acceleration dominates at p > p_*. This result, at first sight appears peculiar, and may indicate that the mean free path is so long, that the diffusion approximation breaks down. It needs to be investigated in more detail (e.g., by investigating the leaky box transport Equation (63) in which both particle escape losses and particle energy gains and losses are accounted for; see the discussion in Appendix E). If one includes a second-order time derivative term in the transport Equation (63) (or its space dependent version given in Webb et al. 2018a, 2018b), then the effects of a finite cosmic ray inertia will limit the maximum momentum change rate allowed. From Section 5.1 and (69) and (70), the maximum allowed particle momentum derived from the second-order derivatives in the telegrapher equation is given by

$$\begin{eqnarray}&&{p}_{\max }\,=\,{p}_{0}\exp \left(\displaystyle \frac{\sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)\quad \mathrm{and}\quad \langle \dot{p}\rangle =\displaystyle \frac{\sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }}p\end{eqnarray} \tag{ 94 }$$

(here we have assumed for the sake of analytical simplicity that we can take Γ as a constant, in order to capture the basic physical effects). The analog of the balance condition between synchrotron losses and shear acceleration energy gains due to cosmic ray viscosity is replaced by the equation

$$\begin{eqnarray}&&\langle \dot{p}{\rangle }_{\max }+\langle \dot{p}{\rangle }_{s}=0,\end{eqnarray} \tag{ 95 }$$

which implies the balance equation

$$\begin{eqnarray}&&\displaystyle \frac{\sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }}p-\displaystyle \frac{{D}_{s}}{\gamma }{p}^{2}=0.\end{eqnarray} \tag{ 96 }$$

Equation (96) implies that losses balance gains at the momentum:

$$\begin{eqnarray}&&{p}_{\mathrm{tel}}^{* }\,=\,\displaystyle \frac{\gamma \sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }{D}_{s}}=\displaystyle \frac{{\gamma }^{2}}{\sqrt{15(1-{\beta }^{2}/3)}}\left(\displaystyle \frac{{du}}{{dr}}\right){m}_{0}c{\tau }_{r}.\end{eqnarray} \tag{ 97 }$$

The momentum ${p}_{\mathrm{tel}}^{* }$ in (97) can also be written in the form:

$$\begin{eqnarray}&&{p}_{\mathrm{tel}}^{* }\,=\,{m}_{0}c\left[\displaystyle \frac{{\gamma }^{2}}{\sqrt{15(1-{\beta }^{2}/3)}}\displaystyle \frac{{\tau }_{r}}{{\tau }_{\mathrm{sh}}}\right].\end{eqnarray} \tag{ 98 }$$

Because ${p}_{\mathrm{tel}}^{* }$ is the maximum allowed momentum in the telegrapher transport equation model (63), we expect that ${p}_{\mathrm{tel}}^{* }\,\gt \,{p}_{\mathrm{sh}}^{* }$ is required by causality. We take the break momentum p*, to be defined as

$$\begin{eqnarray}{p}^{* }=\left\{\begin{array}{ll}\min \left({p}_{\mathrm{sh}}^{* },\,{p}_{\mathrm{tel}}^{* }\right) & \mathrm{if}\quad 0\lt \alpha \lt 1,\\ {p}_{\mathrm{tel}}^{* } & \mathrm{if}\quad \alpha \gt 1.\end{array}\right.\end{eqnarray} \tag{ 99 }$$

In the above analysis we have not taken into account particle losses through the spatial boundaries of the acceleration region. Also, the dispersion term $\langle {\left({\rm{\Delta }}p\right)}^{2} \rangle /(2{\rm{\Delta }}t)$ term has not been taken into account in the above simple argument.

Below, we give explicit formulae for τ_r and τ_sh, which may be used to determine ${p}_{\mathrm{tel}}^{* }$ in (98) for both electrons (e) and protons (p). From (85), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{{re}} & = & 7.75\times {10}^{20}/{\left({B}_{\mu G}\right)}^{2}\ {\rm{s}},\\ {\tau }_{{rp}} & = & 4.8\times {10}^{30}/{\left({B}_{\mu G}\right)}^{2}\ {\rm{s}},\end{array}\end{eqnarray} \tag{ 100 }$$

for the synchrotron radiation timescales for electrons (e) and protons (p), where B_μG is the magnetic field induction in microGauss (i.e., 10⁻⁶ Gauss). The shear acceleration timescale is taken as

$$\begin{eqnarray}&&{\tau }_{\mathrm{sh}}=\displaystyle \frac{L}{u}\sim {\left({du}/{dr}\right)}^{-1},\end{eqnarray} \tag{ 101 }$$

where L is a typical length scale for the shear flow. Using (100) and (101), we obtain the formulae

$$\begin{eqnarray}&&{p}_{\mathrm{tel}}^{* (e)}=1.0\times {10}^{15}\displaystyle \frac{{\gamma }^{2}}{\sqrt{1-{\beta }^{2}/3}}\displaystyle \frac{\beta }{L({kpc}){B}_{\mu G}^{2}}\quad \mathrm{eV}/c,\end{eqnarray} \tag{ 102 }$$

$$\begin{eqnarray}&&{p}_{\mathrm{tel}}^{* (p)}=1.2\times {10}^{28}\displaystyle \frac{{\gamma }^{2}}{\sqrt{1-{\beta }^{2}/3}}\displaystyle \frac{\beta }{L({kpc}){B}_{\mu G}^{2}}\quad \mathrm{eV}/c,\end{eqnarray} \tag{ 103 }$$

as expressions for ${p}_{\mathrm{tel}}^{* }$ for electrons (e) and protons (p). Note that

$$\begin{eqnarray}&&\begin{array}{l}{\sigma }_{{Tp}}={\sigma }_{{Te}}{\left(\displaystyle \frac{{m}_{e}}{{m}_{p}}\right)}^{2},\quad {\tau }_{{rp}}={\tau }_{{re}}{\left(\displaystyle \frac{{m}_{p}}{{m}_{e}}\right)}^{3},\\ \quad \quad {p}_{\mathrm{tel}}^{* (p)}={p}_{\mathrm{tel}}^{* (e)}{\left(\displaystyle \frac{{m}_{p}}{{m}_{e}}\right)}^{4}.\end{array}\end{eqnarray} \tag{ 104 }$$

In the above equations, L(kpc) is the typical jet width in kilo-parsecs (e.g., Webb et al. 2018a). For a sharp gradient shear layer at the edge of the jet, τ_sh = u/L could be much larger than the above estimate (e.g., Ostrowski 1990, 1998, 2002).

From Webb et al. (2018a), we assume a scattering time

$$\begin{eqnarray}&&\tau ={\tau }_{0}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }=\displaystyle \frac{3}{(2-q)(4-q)\chi }\left(\displaystyle \frac{{{\ell }}_{b}}{c}\right){\left(\displaystyle \frac{{p}_{0}c}{{ZeB}{{\ell }}_{b}}\right)}^{\alpha }{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha },\end{eqnarray} \tag{ 105 }$$

to describe the particle scattering in the turbulent magnetic field, based on quasi-linear theory (e.g., Jokipii 1966, 1971). In (105)

$$\begin{eqnarray}&&\alpha =2-q,\end{eqnarray} \tag{ 106 }$$

where q is the power-law spectral index for slab turbulence—that is, it is assumed that the power spectrum of magnetic field fluctuations transverse to the field has the form

$$\begin{eqnarray}&&{P}_{{xx}}(k)={P}_{{xx}}({k}_{b}){\left(\displaystyle \frac{k}{{k}_{b}}\right)}^{-q}\quad \mathrm{for}\quad k\gt {k}_{b},\end{eqnarray} \tag{ 107 }$$

where k_b = 1/ℓ_b is the correlation scale for the turbulence. In (105)

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{k}_{b}{P}_{{xx}}({k}_{b})}{{B}^{2}}=\displaystyle \frac{\delta {B}_{x}^{2}}{{\rm{\Phi }}{B}^{2}},\end{eqnarray} \tag{ 108 }$$

is a parameter of order unity (see Webb et al. 2018a for more detail). From (105) we obtain

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{3}{(2-q)(4-q)\chi }\left(\displaystyle \frac{{{\ell }}_{b}}{c}\right){\left(\displaystyle \frac{{r}_{g0}}{{{\ell }}_{b}}\right)}^{\alpha },\quad {r}_{g0}=\displaystyle \frac{{p}_{0}c}{{ZeB}}.\end{eqnarray} \tag{ 109 }$$

The expression (92) for ${p}_{\mathrm{sh}}^{* }$ is independent of p₀. Using (109) for τ₀ in (92), we obtain the formula

$$\begin{eqnarray}&&{p}_{\mathrm{sh}}^{* }\,=\,{m}_{0}{{cP}}^{\delta },\end{eqnarray} \tag{ 110 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}P & = & \displaystyle \frac{(4+\alpha ){\gamma }^{4}}{15}{C}_{0}{\left(\displaystyle \frac{{m}_{0}{c}^{2}}{{ZeB}{{\ell }}_{b}}\right)}^{\alpha }\displaystyle \frac{{{\ell }}_{b}}{c}\left(\displaystyle \frac{6\pi {m}_{0}c}{{\sigma }_{T}{B}^{2}}\right){\left(\displaystyle \frac{\beta c}{L}\right)}^{2}\\ & \equiv & \displaystyle \frac{(4+\alpha ){\gamma }^{4}{\beta }^{2}}{15}{C}_{0}{\left(\displaystyle \frac{{m}_{0}{c}^{2}}{{eB}{{\ell }}_{b}}\right)}^{(\alpha +2)}\displaystyle \frac{6\pi {{\ell }}_{b}^{3}{r}_{0}}{{L}^{2}{\sigma }_{T}}.\end{array}\end{eqnarray} \tag{ 111 }$$

In the second form for P in (111) we assume Z = 1 and

$$\begin{eqnarray}&&{C}_{0}=\displaystyle \frac{3}{(2-q)(4-q)\chi },\quad \delta =\displaystyle \frac{1}{1-\alpha }.\end{eqnarray} \tag{ 112 }$$

The above expression (111) for P may be written as

$$\begin{eqnarray}&&P=\displaystyle \frac{3}{20}(4+\alpha ){C}_{0}{\gamma }^{4}{\beta }^{2}{\left(\displaystyle \frac{{{\ell }}_{b}}{L}\right)}^{1-\alpha }{\left(\displaystyle \frac{{r}_{0}}{L}\right)}^{-1}{\left(\displaystyle \frac{{m}_{0}{c}^{2}}{{eBL}}\right)}^{\alpha +2}.\end{eqnarray} \tag{ 113 }$$

For the case of electrons, using (113) with α = 1/3 and q = 5/3 (corresponding to Kolmogorov turbulence with q = 5/3), setting χ = 1, and using Equation (110) gives

$$\begin{eqnarray}&&{p}_{\mathrm{sh}}^{* (e)}=3.1\times {10}^{14}{\gamma }^{6}{\beta }^{3}\left(\displaystyle \frac{{{\ell }}_{b}}{L}\right){L}_{\mathrm{kpc}}^{-4/3}{B}_{\mu G}^{-7/2}\ \mathrm{eV}/c.\end{eqnarray} \tag{ 114 }$$

Comparing (114) for ${p}_{\mathrm{sh}}^{* (e)}$ with (102) for ${p}_{\mathrm{tel}}^{* (e)}$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{p}_{\mathrm{tel}}^{* (e)}}{{p}_{\mathrm{sh}}^{* (e)}}=3.3{\gamma }^{-4}{\beta }^{-2}{\left(1-{\beta }^{2}/3\right)}^{-1/2}{L}_{\mathrm{kpc}}^{1/3}{\left(\displaystyle \frac{{{\ell }}_{b}}{L}\right)}^{-1}{B}_{\mu G}^{3/2}.\end{eqnarray} \tag{ 115 }$$

From (115) we deduce

$$\begin{eqnarray}&&{p}_{\mathrm{tel}}^{* (e)}\lt {p}_{\mathrm{sh}}^{* (e)}\ \ \Longleftrightarrow \ \ {\gamma }^{4}{\beta }^{2}\sqrt{1-{\beta }^{2}/3}\gt 3.3{L}_{\mathrm{kpc}}^{1/3}{\left(\displaystyle \frac{{{\ell }}_{b}}{L}\right)}^{-1}{B}_{\mu G}^{3/2}.\end{eqnarray} \tag{ 116 }$$

Thus, the telegrapher equation cut-off is smaller than the Fokker–Planck, systematic acceleration cut-off for sufficiently relativistic radio jets. In other words, the telegrapher equation cut-off ${p}_{\mathrm{tel}}^{* (e)}$ should be used for fast enough jets for which the inequality (116) applies. This constraint is fairly severe if we take for example ℓ_b/L = 0.1, in which case condition (116) reduces to

$$\begin{eqnarray}&&{\gamma }^{4}{\beta }^{2}\sqrt{1-{\beta }^{2}/3}\gt 33{L}_{\mathrm{kpc}}^{1/3}{B}_{\mu G}^{3/2}.\end{eqnarray} \tag{ 117 }$$

In this paper, we have revisited the problem of acceleration of UHECR in relativistic jet shear flows, extending the earlier work of Webb (1990) and Webb et al. (2018a), based on the relativistic diffusive transport equation of Webb (1989). The model assumes a relativistic jet shear flow of the form ${\boldsymbol{u}}=u(r){{\boldsymbol{e}}}_{z}$ along the axis of the jet (the z-axis), in which u(r) is a monotonic decreasing function of cylindrical radius r about the jet axis (i.e., du/dr < 0). The particles are accelerated by elastically scattering in the local fluid frame, and by scattering across the jet shear flow. The particles are assumed to be relativistic for which the particle speed v ≈ c, where c is the speed of light. The mean scattering time τ(r, p) for particles with momentum p at cylindrical radius r is assumed to have the functional form

$$\begin{eqnarray}&&\tau (r,p)={\tau }_{0}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }\displaystyle \frac{g(\xi )}{{rd}\xi /{dr}},\end{eqnarray} \tag{ 118 }$$

where g(ξ) is a given function of the jet rapidity,

$$\begin{eqnarray}&&\xi (r)=\displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{1+\beta (r)}{1-\beta (r)}\right),\quad \mathrm{where}\quad \beta (r)=\displaystyle \frac{u(r)}{c}\end{eqnarray} \tag{ 119 }$$

is relativistic β of the flow and c is the speed of light. The monoenergetic source Green's function of the steady-state transport equation obtained by Webb (1990) and Webb et al. (2018a) assumed that g(ξ) was constant. This lead to solutions of the transport equation that have the property $\tau (r,p)\to \infty$ as r → 0 for most models of the shear flow. In this paper we have explored more general choices of g(ξ) chosen, so that

$$\begin{eqnarray}&&r\displaystyle \frac{d\xi }{{dr}}=g(\xi )=-k{\left({\xi }_{0}-\xi \right)}^{\delta },\quad \tau ={\tau }_{0}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }.\end{eqnarray} \tag{ 120 }$$

We explored in detail Green's function solutions and monoenergetic spectrum solutions for the case δ = 1, in which the eigenfunction solutions for the Green's function involves zeroth order Bessel functions J₀(λη), where η = η(r) depends only on the radial profile of the jet speed u(r). There are in fact other Bessel eigenfunction series Green's functions that can be obtained for cases where δ > 1 (see Appendix B for further details).

The solution examples (Figures 1–8) are very similar to those presented in Webb et al. (2018a). This illustrates that the solution details do not depend strongly on the detailed behavior of the mean scattering time τ(r, p) on the axis of the jet as r → 0 (in Webb et al. 2018a, $\tau (r,p)\to \infty$ as r → 0, but in the present model $\tau (r,p)={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$ does not depend on r). For an ultra-relativistic jet, the spectral index of shear accelerated particles (${f}_{0}\propto {p}^{-{\mu }_{\infty }}$), ${\mu }_{\infty }\to 3+\alpha$, and in the non-relativistic limit β₀ → 0, ${\mu }_{\infty }\to \sqrt{5}{j}_{\mathrm{0,1}}/{\beta }_{02}$, where β₀₂ is the relative velocity of β₀ and β₂ and j_0,1 is the first zero of the ordinary Bessel function of order zero J₀(z) (j_0,1 ∼ 2.4048). The spectral index ${\mu }_{\infty }$ is displayed as a function of β₀ in Figure 3 (see also Rieger & Duffy 2006, for discussion of the spectral index of particles accelerated in shear flows by cosmic ray viscosity). For α = 1, ${\mu }_{\infty }\to 4$ as β₀ → 1, which is similar to the spectral index of cosmic rays accelerated by the first order Fermi mechanism in astrophysical shocks (Axford et al. 1977; Krymsky 1977; Bell 1978a, 1978b; Blandford & Ostriker 1978). Allowing for a back-flowing cocoon with an increased overall shear of the flow leads to harder particle spectra.

In Section 3, steady-state Green's function solutions of the relativistic diffusive transport equation (Webb 1989) were obtained for a cylindrically symmetric jet in which the flow velocity u = u(r) e_z. The Green's function consists of an eigenfunction expansion in terms of J₀(λ_nη) Bessel functions where η = η(r) describes the radial dependence of the solution and the λ_n are eigenvalues. The Green's function corresponds to steady injection of particles with momentum p = p₀ at radius r = r₁ from the jet axis, with free escape boundary conditions at radius r = r₂ (r₂ > r₁). The flow profile u(r) is determined in part by requiring the mean scattering time τ(r, p) to be independent of r (i.e., $\tau (r,p)={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$).

Green's formula for the transport equation was obtained in Section 4. It gives the solution of boundary problems for f₀(r, p) in terms of the Green's functions of Section 3, and in terms of the boundary and initial data. In particular, the Green's formula gives solutions of boundary value problems in which a spectrum of particles is specified at the outer boundary of the jet at r = r₂ and in terms of sources of particles inside the jet at radii 0 < r < r₂. The Green's formula was used to obtain monoenergetic spectrum solutions in which ${f}_{0}(r,p)\,\to {N}_{g}\delta (p-{p}_{0})/(4\pi {p}_{0}^{2})$ as r → r₂, and with no sources inside the jet shear flow region 0 < r < r₂. In these solutions there is a radial influx of particles with momentum p = p₀ and an outflux of particles at higher and lower momenta that have either have been accelerated (p > p₀) or decelerated (p < p₀) by their interaction with the shear flow.

Astrophysical applications of the model to the acceleration of UHECR to ultra-high energies in extragalactic radio jets were discussed in Liu et al. (2017) and Webb et al. (2018a). Constraints on the acceleration mechanism were used to delimit the width ${\rm{\Delta }}L$ of the radio jet and magnetic field strength B needed to constrain and accelerate E = 10¹⁸ eV protons in the jet, taking into account proton synchrotron losses and other constraints (e.g., the width of the jet ΔL needs to be greater than the particle mean free path λ in order that the particles are confined to the jet and can interact with the shear flow; the particle acceleration time needs to be less than the dynamical advection time in order for the particles to be accelerated by the shear before being advected out of the region). Equation (99) of Webb et al. (2018a) gives a formula for the gyro-radius of the particles that is too big by a factor of about 10. Correcting for this error results in estimates for the allowed jet width to be approximately one tenth of those given in Webb et al. (2018a). These corrections to (99) et seq. of Webb et al. (2018a) are discussed in detail in Appendix F. Constraints on ΔL and B for the case of isotropic diffusion from Webb et al. (2018a; adjusted for the above error) for particles with E = 10¹⁸ eV (i.e., 1 EeV particles), are 1 kpc < ΔL < 10⁵ kpc for B =1 μG and 10 pc < ΔL < 10³ kpc for B = 100 μG for a radio jet with a Lorentz factor γ_j = 1.1 (β_j = 0.4166). The length of the jet d = 100 ΔL was assumed. AGN jets of these approximate size and magnetic field constraints (Liu et al. 2017 and Webb et al. 2018a) imply that the sources MKN 501 and MKN 421; Dermer (2007), Sahayanathan (2009), Abbasi et al. (2014), Caprioli (2015) are possible sources of EeV protons accelerated by the cosmic ray viscous shear acceleration mechanism.

In Appendix C we discuss how the analysis is altered by assuming that the particles have a finite diffusion coefficient in the region r > r₂, but there is no particle acceleration by shear in the region r > r₂ outside the jet. Assuming that that there is a spectrum of particles ${f}_{0}(\infty ,p)$ as $r\to \infty$ and that there are no sources of particles, in r > r₂, the boundary condition at r = r₂ reduces to the equation

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0}({r}_{2},p)+{r}_{2}\kappa ({r}_{2},p)\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}\left({\int }_{{r}_{2}}^{\infty }\displaystyle \frac{{dr}^{\prime} }{r^{\prime} \kappa (r^{\prime} ,p)}\right)\\ \quad =\,{f}_{0}(\infty ,p).\end{array}\end{eqnarray} \tag{ 121 }$$

Effectively the solutions of Section 4, in which f₀(r₂, p) are specified on a free escape boundary at r = r₂, assume that $\kappa (r,p)\to \infty$ in r > r₂. More precisely, if κ ∝ p^αr^s in r > r₂, the boundary condition (121) reduces to the equation

$$\begin{eqnarray}&&{f}_{0}({r}_{2},p)+\displaystyle \frac{{r}_{2}}{s}\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}={f}_{0}(\infty ,p).\end{eqnarray} \tag{ 122 }$$

The boundary condition (122) is a mixed Dirichlet–Von Neumann boundary condition. For s ≫ 1 (122) is approximately the Dirichlet boundary condition:

$$\begin{eqnarray}&&{f}_{0}({r}_{2},p)\approx {f}_{0}(\infty ,p),\end{eqnarray} \tag{ 123 }$$

which is essentially the free escape boundary condition used in Section 4. The effect of including a finite κ in r > r₂ will result in a harder spectrum of accelerated particles in r < r₂ because there is a finite probability that the particles leaving the acceleration region (r < r₂) will return to the region r < r₂ to be further accelerated by the shear flow. This is an interesting extension of the solutions presented in Sections 3 and 4, which lies beyond the scope of the present paper.

Section 5 discusses the effect of energy losses and particle escape from the acceleration region by using a leaky box type transport equation, including the effects of cosmic ray inertia, momentum space diffusion energy changes due to shear acceleration, and second-order Fermi acceleration due to Alfvén waves. From an analysis of the characteristic manifold and bi-characteristics of a telegrapher type transport equation we obtain a first order ordinary differential equation for p = p(t), which gives the maximum possible momentum of the particles in the model. It is shown in Appendix D that if one uses a diffusion equation analysis to calculate the mean particle momentum $\langle p(t)\rangle$ for the particle acceleration, one obtains a contradictory result that particles can be accelerated to an infinite momentum in a finite time. This non-physical result can be circumvented by including a second-order time derivative term in a telegrapher type transport equation that incorporates the effect of the cosmic ray inertia. Appendix E discusses a steady-state leaky box model for particle acceleration due to cosmic ray viscosity in the presence of synchrotron losses, and in particular we discuss the dependence of the solution on the momentum dependence of the particle scattering time $\tau ={\tau }_{0}{\left(p/{p}_{0}\right)}^{\alpha }$ through the parameter α, for the separate cases (a) 0 < α < 1, (b) α = 1 (Bohm diffusion), and (c) α > 1.

As an application of the p_max(t) formula implied by the telegrapher equation, we estimate a lower limit to the time for particles to be accelerated in a radio-jet shear flow, with a typical width of L ∼ 1 kpc in which the particles are injected with a momentum p₀ = 1 GeV/c and obtain a momentum of p = 1 EeV/c on a timescale of the order of 10⁵–10⁶ yr for a mildly relativistic jet. The estimate of p_max(t) based on the telegrapher equation bi-characteristics in the example turns out to be about 100 light-crossing times of the width of the jet. The analysis assumes that the transport Equation (63) has space averaged transport coefficients of the same form as the space dependent transport equation (e.g., Webb et al. 2018a, Appendix E). The averaged transport coefficients depend on the details of the model parameters for the jet.

Estimates of acceleration times (based on the telegrapher equation bi-characteristics, and also on the Fokker–Planck mean drift momentum term) are used to investigate if the particle momentum losses balance the momentum gains at some critical momentum p = p_* in which the particle energy losses exceed the momentum gains at momenta p > p_*, in which case there is a break in the momentum spectrum and a fall off of the particle momentum spectrum at p > p*. In fact if there is very strong particle acceleration, then there can occur a momentum pile up in the spectrum due to strong deceleration of particles above p* and strong acceleration below p* (e.g., Kardashev 1962). Kardashev (1962) studies a variety of energy loss mechanisms (synchrotron losses, expansion losses, nuclear losses, Bremsstrahlung losses, ionization losses, and inverse Compton losses) and stochastic acceleration of the particles, with application to a variety of sources (e.g., the Crab Nebula, SNR Cass A, Jet in NGC 4486, Cygnus A). Similar ideas have been used to describe diffusive shock acceleration of electrons subject to synchrotron losses in astrophysical shocks (e.g., Webb et al. 1984; Fritz & Webb 1990; Blasi 2010). We find that for sufficiently relativistic flows, if the balance point in momentum, based on the telegrapher equation approach for the balance or break point in the spectrum ${p}_{\mathrm{tel}}^{* }$, is less than the cut-off based on the diffusive transport equation, ${p}_{\mathrm{sh}}^{* }$, then the solution effectively must cut off at ${p}_{\mathrm{tel}}^{* }$ and not at ${p}_{\mathrm{sh}}^{* }$, as assumed in earlier analyses.

It was pointed out by Webb et al. (2018a) that time-dependent, diffusive transport equations for the acceleration of particles involving second-order momentum diffusion are non-causal, as they predict that particles can attain an infinite energy in a finite time. These defects of diffusive transport equations can be partially overcome by generalizing the diffusive transport equation to include a second-order time derivative term representing cosmic ray inertia (e.g., Fisk & Axford 1969) to obtain a telegrapher equation for cosmic ray transport. Analysis of the telegrapher equation characteristic manifolds and bi-characteristics suggests the formula (Webb et al. 2018a):

$$\begin{eqnarray}&&{p}_{\max }(t)={p}_{0}\exp \left[\displaystyle \frac{\sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }}{\left({t}^{2}-\displaystyle \frac{{\tilde{\gamma }}^{2}{\left(r-{r}_{0}\right)}^{2}}{{V}_{\mathrm{cr}}^{2}}\right)}^{1/2}\right],\end{eqnarray} \tag{ 124 }$$

for the maximum particle momentum for particles accelerated by the cosmic ray viscous shear acceleration mechanism in the relativistic jet shear flow for particles injected into the flow with momentum p = p₀ at radius r = r₀ at time t = 0. Here ${V}_{\mathrm{cr}}=c/\sqrt{3}$ is the cosmic ray sound speed for isotropic relativistic particles (e.g., Webb 1987a; Kirk & Webb 1988); $\tilde{\gamma }=\gamma \sqrt{1-{\beta }^{2}/3}$, where γ is the relativistic Lorentz factor of the flow. For simplicity, β, γ, and Γ were assumed constant when integrating the bi-characteristics (numerical integration of the bi-characteristics could in principle be used to obtain more accurate estimates of p_max). A preliminary account of these ideas was presented in Webb et al. (2018a) and in a poster paper at the fall AGU meeting in 2018 December (Webb et al. 2018b). We hope to address some of these issues in more detail in future publications.

G.M.W. acknowledges discussions with Frank Rieger, who pointed out the numerical error in the expression for the particle gyro-radius in Equation (99) of Webb et al. (2018a). We thank the referee for a very thorough report and for pointing out inconsistencies and errors in the original manuscript. G.M.W. and G.L. were supported in part by NASA grant 80NSSC19K0075. P.M. acknowledges the support of a NASA Earth and Space Science Fellowship Grant 16-HELIO16F-0022. J.A.L.R. is supported in part by NASA grants NNX15A165G, 80NSSC19K0276, and NSF-DOE grant PHY-1707247. G.P.Z. has been supported in part by the NSF EPSCoR RII-Track-1 Cooperative Agreement OIA-1655280 and an NSF/DOE Partnership in Basic Plasma Science and Engineering via NSF grant PHY-1707247. A.F.B. acknowledges valuable insights and discussions with Prof. Peter Biermann on UHECR and candidate sources and radio galaxies.

In this appendix we derive the Green's function solution (28) of the steady-state transport Equation (26), including the effects of particle acceleration due to cosmic ray viscosity and fluid shear. We derive the Green's function by two different methods. The first method uses the Fourier transform of the partial differential Equation (26) with respect to T. The second method uses a Fourier–Bessel expansion consisting of separated solutions of (26), which take into account the Delta function source term in (26), in which the Delta function $\delta (\eta -{\eta }_{0})$ is expanded in terms of the Bessel equation eigenfunctions for the problem.

A.1. Fourier Transform Approach

We introduce the Fourier transform with respect to T in the form:

$$\begin{eqnarray}&&{\bar{f}}_{0}=\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }{f}_{0}(r,p)\exp (i\nu T)\ d\nu .\end{eqnarray} \tag{ 125 }$$

Taking the Fourier transform of (26) with respect to T gives the ordinary differential equation,

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{\bar{f}}_{0}}{d{\eta }^{2}}+\displaystyle \frac{1}{\eta }\displaystyle \frac{d{\bar{f}}_{0}}{d\eta }-{{\rm{\Lambda }}}^{2}{\bar{f}}_{0}=\bar{Q}=\displaystyle \frac{3{N}_{0}\delta (\xi -{\xi }_{1})}{\sqrt{2\pi }8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}g({\xi }_{1})},\end{eqnarray} \tag{ 126 }$$

where

$$\begin{eqnarray}&&{{\rm{\Lambda }}}^{2}=\displaystyle \frac{{\left(\nu +{ia}\right)}^{2}+{a}^{2}}{5},\quad a=\left(\displaystyle \frac{3+\alpha }{2}\right).\end{eqnarray} \tag{ 127 }$$

Notice that

$$\begin{eqnarray}&&{\rm{\Lambda }}=\displaystyle \frac{{\left(\nu {{\prime} }^{2}+{a}^{2}\right)}^{1/2}}{\sqrt{5}},\quad \nu ^{\prime} =\nu +{ia}.\end{eqnarray} \tag{ 128 }$$

The homogeneous Equation (126) has solutions of the form

$$\begin{eqnarray}&&y={a}_{1}{I}_{0}({\rm{\Lambda }}\eta )+{a}_{2}{K}_{0}({\rm{\Lambda }}\eta )\equiv {a}_{1}{y}_{1}+{a}_{2}{y}_{2},\end{eqnarray} \tag{ 129 }$$

where η = ξ₀ − ξ, and I₀(z) and K₀(z) are modified Bessel functions of the first and second kind of order zero (Abramowitz & Stegun 1965, 9.6.1, p. 374). The general solution of the inhomogeneous ordinary differential Equation (126) is given by the general formula

$$\begin{eqnarray}\begin{array}{rcl}{\bar{f}}_{0} & = & {y}_{1}\left({c}_{1}-{\displaystyle \int }^{\eta }\displaystyle \frac{\bar{Q}{y}_{2}(\eta ^{\prime} )}{{W}_{\eta ^{\prime} }({y}_{1},{y}_{2})}\ d\eta ^{\prime} \right)\\ & & +{y}_{2}\left({c}_{2}+{\displaystyle \int }^{\eta }\displaystyle \frac{\bar{Q}{y}_{1}(\eta ^{\prime} )}{{W}_{\eta ^{\prime} }({y}_{1},{y}_{2})}\ d\eta ^{\prime} \right)\end{array}\end{eqnarray} \tag{ 130 }$$

(Morse & Feshbach 1953, Vol. 1, p. 530, 5.2.19), where W_η(y₁, y₂) = y₁dy₂/dη − y₂dy₁/dη is the Wronskian of the two independent solutions y₁ and y₂ of the homogeneous Equation (126), which are listed in (129). In (128) and (129) we use the Re(Λ) > 0 branch of the multi-valued function Λ, and we use the variable $\nu ^{\prime} =\nu +{ia}$ rather than ν in the analysis.

Using the properties of I₀(z) and K₀(z) to fit the boundary conditions as r → 0 and as r → r₂, we obtain the solution for ${\bar{f}}_{0}$ in the form

$$\begin{eqnarray}\begin{array}{rcl}{\bar{f}}_{0} & = & \displaystyle \frac{{ \mathcal A }{\eta }_{1}}{{I}_{0}({\rm{\Lambda }}{\eta }_{2})}\left[{I}_{0}({\rm{\Lambda }}{\eta }_{1})F(\eta ,{\eta }_{2})H(r-{r}_{1})\right.\\ & & +\left.{I}_{0}({\rm{\Lambda }}\eta )F({\eta }_{1},{\eta }_{2})H({r}_{1}-r)\right],\end{array}\end{eqnarray} \tag{ 131 }$$

where

$$\begin{eqnarray}&&F(\eta ,{\eta }_{2})={K}_{0}({\rm{\Lambda }}\eta ){I}_{0}({\rm{\Lambda }}{\eta }_{2})-{I}_{0}({\rm{\Lambda }}\eta ){K}_{0}({\rm{\Lambda }}{\eta }_{2}),\end{eqnarray} \tag{ 132 }$$

$$\begin{eqnarray}&&{ \mathcal A }=\displaystyle \frac{3{N}_{0}}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k\sqrt{2\pi }},\quad \eta ={\xi }_{02}{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{k}.\end{eqnarray} \tag{ 133 }$$

The required solution for f₀(r, p) is obtained by using the inverse Fourier transform formula:

$$\begin{eqnarray}&&{f}_{0}=\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }\exp (-i\nu T){\bar{f}}_{0}(\eta ,\nu )\ d\nu .\end{eqnarray} \tag{ 134 }$$

To carry out the inverse Fourier transform (134), it is useful to note that ${\bar{f}}_{0}$ has simple poles in the complex λ-plane at points where

$$\begin{eqnarray}&&{I}_{0}({\rm{\Lambda }}{\eta }_{2})=0\quad \mathrm{or}\quad {J}_{0}(i{\rm{\Lambda }}{\eta }_{2})=0.\end{eqnarray} \tag{ 135 }$$

These poles occur at the points

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{n} & = & \pm i{\lambda }_{n}=\pm i\displaystyle \frac{{j}_{0,n}}{{\eta }_{2}}\\ & \equiv & \pm i\displaystyle \frac{{j}_{0,n}}{{\xi }_{02}},\quad {\lambda }_{n}\equiv | {{\rm{\Lambda }}}_{n}| ,\quad n=1,2,3,\,\ldots ,\end{array}\end{eqnarray} \tag{ 136 }$$

where ${j}_{0,n}$ is the nth zero of J₀(z)—that is, ${J}_{0}({j}_{0,n})=0$ (n a positive integer). The value of ν' = ν + ia at these points are given by

$$\begin{eqnarray}&&{\nu }_{n}^{{\prime} }=\pm i{\chi }_{n},\quad {\chi }_{n}={\left({a}^{2}+5\displaystyle \frac{{j}_{0,n}^{2}}{{\xi }_{02}^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 137 }$$

The multi-valued complex function Λ = Λ(ν') is chosen to have a branch cut from ν' = ia to $\nu ^{\prime} =i\infty$ along the Im(ν') > 0 axis, and also a cut from ν' = −ia to $\nu ^{\prime} =-i\infty$ along the Im(ν') < 0 axis. It turns out that the transform (131) does not jump in value (i.e., is continuous) across the branch cut.

Using standard complex variable theory (i..e. Cauchy's theorem) to invert the solution (131) for ${\bar{f}}_{0}$ gives the Green's function solution for f₀(r, p) in the form

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0} & = & C{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{-a}\displaystyle \sum _{n=1}^{\infty }\left(\displaystyle \frac{{j}_{0,n}{Y}_{0}({j}_{0,n})}{{J}_{1}({j}_{0,n}){\chi }_{n}}\right)\\ & & \times \,{J}_{0}({\lambda }_{n}\eta ){J}_{0}({\lambda }_{n}{\eta }_{1})\exp (-{\chi }_{n}| T| ),\end{array}\end{eqnarray} \tag{ 138 }$$

where

$$\begin{eqnarray}&&C=\displaystyle \frac{15{N}_{0}}{16{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k{\eta }_{2}^{2}}.\end{eqnarray} \tag{ 139 }$$

The solution (138) is equivalent to the solution (28). To show that (138) is equivalent to (28), note that the Wronskian $W({J}_{0},{Y}_{0})$ of J₀(z) and Y₀(z) is given by

$$\begin{eqnarray}&&W({J}_{0},{Y}_{0})={J}_{0}(z){Y}_{0}^{{\prime} }(z)-{Y}_{0}(z){J}_{0}^{{\prime} }(z)=\displaystyle \frac{2}{\pi z}\end{eqnarray} \tag{ 140 }$$

(Abramowitz & Stegun 1965, formula 9.1.60, p. 360). Also note that

$$\begin{eqnarray}&&{J}_{0}^{{\prime} }(z)=-{J}_{1}(z),\quad {Y}_{0}^{{\prime} }(z)=-{Y}_{1}(z)\end{eqnarray} \tag{ 141 }$$

(Abramowitz & Stegun 1965, formula 9.1.28, p. 361). Using (140) and (141) with z = j_0,n, we obtain

$$\begin{eqnarray}&&{Y}_{0}({j}_{0,n}){J}_{1}({j}_{0,n})-{J}_{0}({j}_{0,n}){Y}_{1}({j}_{0,n})=\displaystyle \frac{2}{\pi {j}_{0,n}},\end{eqnarray} \tag{ 142 }$$

which implies

$$\begin{eqnarray}&&{j}_{0,n}{Y}_{0}({j}_{0,n})=\displaystyle \frac{2}{\pi {J}_{1}({j}_{0,n})}.\end{eqnarray} \tag{ 143 }$$

Using the result (143) in (138), the solution (138) reduces to the Green's function (28).

A.2. Fourier–Bessel Series Approach

An alternative derivation of the Green's function (28) (e.g., Webb & Gleeson 1977) is to first note that the problem has a discrete i.e., countable eigenspectrum {λ_n}, where λ_n = j_0,n/η₂ and the j_0,n are the zeros of J₀(x). This spectrum arises because we impose the boundary condition that f₀(r₂, p) = 0 at r = r₂ and impose a regularity condition on the particle diffusive flux as r → 0. The transport Equation (26) admits separable solutions of the form

$$\begin{eqnarray}&&{f}_{0}=\displaystyle \sum _{n=1}^{\infty }{J}_{0}({\lambda }_{n}\eta ){h}_{n}(T),\end{eqnarray} \tag{ 144 }$$

where h_n(T) is to be determined.

To determine h_n(T), first note that we can write

$$\begin{eqnarray}&&\delta (\eta -{\eta }_{1})=\displaystyle \sum _{n=1}^{\infty }{c}_{n}{J}_{0}({\lambda }_{n}\eta ),\end{eqnarray} \tag{ 145 }$$

which is a Fourier–Bessel expansion for δ(η − η₁). By using Bessel's equation, we obtain the relations

$$\begin{eqnarray}&&({\lambda }_{m}^{2}-{\lambda }_{n}^{2}){\int }_{0}^{{\eta }_{2}}\eta {J}_{0}({\lambda }_{m}\eta ){J}_{0}({\lambda }_{n}\eta )d\eta ={\left[{W}_{{mn}}\eta \right]}_{0}^{{\eta }_{2}},\end{eqnarray} \tag{ 146 }$$

where

$$\begin{eqnarray}&&{W}_{{mn}}={g}_{m}(\eta )\displaystyle \frac{{{dg}}_{n}}{d\eta }-{g}_{n}(\eta )\displaystyle \frac{{{dg}}_{m}}{d\eta }\quad \mathrm{and}\quad {g}_{n}(\eta )={J}_{0}({\lambda }_{n}\eta ).\end{eqnarray} \tag{ 147 }$$

Here W_mn is the Wronskian of g_m(η) and g_n(η). From (146) to (147), we obtain the orthogonality relations for the {g_n(η)} eigenfunctions—namely,

$$\begin{eqnarray}&&{\int }_{0}^{{\eta }_{2}}\eta {J}_{0}({\lambda }_{m}\eta ){J}_{0}({\lambda }_{n}\eta )\ d\eta ={N}_{n}{\delta }_{{mn}},\end{eqnarray} \tag{ 148 }$$

where δ_mn is the Kronecker delta symbol and

$$\begin{eqnarray}\begin{array}{rcl}{N}_{n} & = & \mathop{\mathrm{lim}}\limits_{{\lambda }_{m}\to {\lambda }_{n}}\displaystyle \frac{1}{{\lambda }_{m}^{2}-{\lambda }_{n}^{2}}{[{W}_{{mn}}\eta ]}_{{\eta }_{2}}=\displaystyle \frac{{\eta }_{2}}{2{\lambda }_{n}}{\left(\displaystyle \frac{\partial {g}_{m}(\eta )}{\partial {\lambda }_{m}}\displaystyle \frac{\partial {g}_{n}(\eta )}{\partial \eta }\right)}_{{\eta }_{2}}\\ & = & \displaystyle \frac{1}{2}{\eta }_{2}^{2}{\left[{J}_{1}({j}_{0,n})\right]}^{2}\end{array}\end{eqnarray} \tag{ 149 }$$

are the normalization constants.

To determine the c_n, multiply (145) by ηJ₀(λ_mη) and integrate over η from η = 0 to η = η₂, using (148). This gives

$$\begin{eqnarray}&&{c}_{n}={\eta }_{1}{J}_{0}({\lambda }_{n}{\eta }_{1})/{N}_{n}\end{eqnarray} \tag{ 150 }$$

for the constants {c_n} in (145). Next we substitute the eigenfunction expansion (144) into the transport Equation (26) and use the Delta function representation (145) for δ(η − η₁) in the source term in (26). We obtain the equation

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}{h}_{n}}{{{dT}}^{2}}+2a\displaystyle \frac{{{dh}}_{n}}{{dT}}-5{\lambda }_{n}^{2}{h}_{n}(T)=-5{{Ac}}_{n}\delta (T),\end{eqnarray} \tag{ 151 }$$

for the h_n(T). Here we have written the source term in (26) in the form −Aδ(T)δ(η − η₁), where

$$\begin{eqnarray}&&A=\displaystyle \frac{3{N}_{0}}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k{\eta }_{1}}.\end{eqnarray} \tag{ 152 }$$

The homogeneous Equation (151) has independent solutions,

$$\begin{eqnarray}&&{y}_{1}=\exp \left[-(a+{\chi }_{n})T\right],\quad {y}_{2}=\exp \left[({\chi }_{n}-a)T\right].\end{eqnarray} \tag{ 153 }$$

With Wronskian,

$$\begin{eqnarray}&&W({y}_{1},{y}_{2})={y}_{1}{y}_{2}^{{\prime} }-{y}_{2}{y}_{1}^{{\prime} }=2{\chi }_{n}\exp (-2{aT}).\end{eqnarray} \tag{ 154 }$$

Using the general formula (130) to solve the inhomogeneous Equation (151) and requiring h_n(T) to be bounded as $T\to \pm \infty$, we obtain the required solution of (151) in the form:

$$\begin{eqnarray}&&{h}_{n}(T)=\displaystyle \frac{5{{Ac}}_{n}}{2{\chi }_{n}}\exp [-{aT}-{\chi }_{n}| T| ].\end{eqnarray} \tag{ 155 }$$

Substituting the solution (155) for h_n(T) in (144) and using (152) for A, and (149) and (150) for the c_n, the solution (144) reduces to the Green's function (28).

In this appendix, we indicate how the Bessel function eigenfunction expansion for the steady-state Green's function (28) can be generalized for other choices of the function g(ξ) used to define the mean free time τ(r, p) in (1) and (20). In particular, we consider the family of solutions of the transport Equation (21) that are obtained by choosing

$$\begin{eqnarray}&&g(\xi )=-k{\left({\xi }_{0}-\xi \right)}^{\delta },\quad \delta \geqslant 1),\end{eqnarray} \tag{ 156 }$$

and by choosing

$$\begin{eqnarray}&&\displaystyle \frac{g(\xi )}{{rd}\xi /{dr}}=1.\end{eqnarray} \tag{ 157 }$$

The solutions in Sections 3 and 4 correspond to the choice δ = 1 in (156) and (157) (see also Appendix A).

From (156) and (157), we obtain the differential equation

$$\begin{eqnarray}&&r\displaystyle \frac{d\xi }{{dr}}=-k{\left({\xi }_{0}-\xi \right)}^{\delta },\quad \delta \geqslant 1,k\gt 0.\end{eqnarray} \tag{ 158 }$$

To integrate (158), we use the auxiliary variable

$$\begin{eqnarray}&&\eta ={\xi }_{0}-\xi ,\end{eqnarray} \tag{ 159 }$$

where ξ₀ = ξ(0). We consider solutions with η > 0 so that dξ/dr < 0, and ξ(r) is then a monotonic decreasing function of r. Integration of (157), by setting ξ₂ = ξ(r₂), results in the solution of (158) as

$$\begin{eqnarray}\begin{array}{rcl}\xi & = & {\xi }_{0}-{\left[{\xi }_{02}^{-2m}+2{mk}\mathrm{ln}\left({r}_{2}/r\right)\right]}^{-1/(2m)},\quad m=\displaystyle \frac{(\delta -1)}{2},\\ \eta & = & {\left[{\xi }_{02}^{-2m}+2{mk}\mathrm{ln}\left({r}_{2}/r\right)\right]}^{-1/(2m)},\quad {\xi }_{02}={\xi }_{0}-{\xi }_{2},\end{array}\end{eqnarray} \tag{ 160 }$$

where we assume $m=(\delta -1)/2\gt 0$. For the above solution,

$$\begin{eqnarray}&&\beta (r)=\tanh [\xi (r)]\end{eqnarray} \tag{ 161 }$$

is a monotonic decreasing function of r.

Using the form (156) for g(ξ) in the transport Equation (21) results in the equation

$$\begin{eqnarray}\begin{array}{rcl} & & \displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {\eta }^{2}}+\displaystyle \frac{2m+1}{\eta }\displaystyle \frac{\partial {f}_{0}}{\partial \eta }+\displaystyle \frac{1}{5}\left[\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {T}^{2}}+(3+\alpha )\displaystyle \frac{\partial {f}_{0}}{\partial T}\right]\\ & & \quad =\displaystyle \frac{3{N}_{0}\delta (T)\delta (\xi -{\xi }_{1})}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}g({\xi }_{1})}.\end{array}\end{eqnarray} \tag{ 162 }$$

The homogeneous partial differential Equation (162) admits separated solutions of the form

$$\begin{eqnarray}&&{f}_{0}=\exp ({sT}){\eta }^{-m}\left[{{cJ}}_{m}(\lambda \eta )+{{dY}}_{m}(\lambda \eta )\right],\end{eqnarray} \tag{ 163 }$$

where c and d are arbitrary constants, and J_m(z) and Y_m(z) are Bessel functions of the first and second kind of order m, where

$$\begin{eqnarray}&&s=-a\pm \sqrt{{a}^{2}+5{\lambda }^{2}}\quad \mathrm{where}\quad a=\displaystyle \frac{3+\alpha }{2}.\end{eqnarray} \tag{ 164 }$$

It is clear from the separated solutions (163) that steady-state Green's functions of (162) can be constructed using the separated solutions (163), following the approach of Appendix A. In particular we obtain the Green's function analogous to (28) in the form

$$\begin{eqnarray}&&{f}_{0}=\hat{D}{\eta }^{-m}\exp (-{aT})\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{J}_{m}({\lambda }_{n}\eta ){J}_{m}({\lambda }_{n}{\eta }_{1})}{{[{J}_{m+1}({j}_{m,n})]}^{2}}\displaystyle \frac{\exp (-{\chi }_{n}| T| )}{{\chi }_{n}},\end{eqnarray} \tag{ 165 }$$

where J_m(z) are ordinary Bessel functions of the first kind of order m, and the eigenvalues {λ_n} are the positive zeros of the eigenvalue equation

$$\begin{eqnarray}&&{J}_{m}({\lambda }_{n}{\eta }_{2})=0,\quad {\rm{i}}.{\rm{e}}.,\quad {\lambda }_{n}=\displaystyle \frac{{j}_{m,n}}{{\eta }_{2}},\quad n=1,2,\,\ldots ,\end{eqnarray} \tag{ 166 }$$

where the {j_m,n} are the positive zeros of J_m(z)—that is, J_m(j_m,n) = 0 for n = 1,2,3, .... In (165),

$$\begin{eqnarray}\begin{array}{rcl}m & = & \displaystyle \frac{\delta -1}{2},\quad {\chi }_{n}={\left({a}^{2}+5{\lambda }_{n}^{2}\right)}^{1/2}={\left({a}^{2}+5{j}_{m,n}^{2}/{\eta }_{2}^{2}\right)}^{1/2},\,\\ \hat{D} & = & \displaystyle \frac{15{N}_{0}{\eta }_{1}^{-m}}{8{\pi }^{2}{p}_{0}^{3}{c}^{2}{\tau }_{0}k{\eta }_{2}^{2}},\quad {\eta }_{2}={\xi }_{0}-{\xi }_{2}={\xi }_{02},\\ \eta & = & {\xi }_{0}-\xi ={\left[{\eta }_{2}^{-2m}+2{mk}\mathrm{ln}({r}_{2}/r)\right]}^{-1/(2m)},\\ \beta (r) & = & \tanh (\xi ),\ \xi ={\xi }_{0}-\eta .\end{array}\end{eqnarray} \tag{ 167 }$$

The solution (165) reduces to the solution (28) in the limit as m → 0. Second, as $p\to \infty$, the solution (165) has the asymptotic spectrum ${f}_{0}\propto {p}^{-{\mu }_{(m)}}$, where

$$\begin{eqnarray}&&{\mu }_{(m)}=a+{\left({a}^{2}+5{j}_{m,1}^{2}/{\eta }^{2}\right)}^{1/2}.\end{eqnarray} \tag{ 168 }$$

This spectral index is larger than the spectral index ${\mu }_{\infty }\equiv {\mu }_{(0)}$ given by (40). This indicates that the flow profile (167) for β(r) for $m\ne 0$ is less efficient for accelerating particles than the case for m = 0. The detailed properties of the solution (165) for $m\ne 0$ will not be explored further in the present paper.

In this appendix, we discuss a modified boundary condition to apply to the particle momentum spectrum f₀(r, p) and ∂f₀/∂r at the edge of the jet at cylindrical radius r = r₂. The boundary condition used in (46) was that f₀(r₂, p) is specified at the edge of the jet at r = r₂. This boundary condition does not explicitly address the particle transport in the region r > r₂. It could be interpreted as the solution of the transport Equation (44), with $\kappa \to \infty$ in r > r₂.

However, if one assumes there is no particle acceleration in r > r₂, but the diffusion coefficient is given in r > r₂, then the particle transport Equation (44) in r > r₂ reduces to

$$\begin{eqnarray}&&-\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r\kappa \displaystyle \frac{\partial {f}_{0}}{\partial r}\right)=0,\end{eqnarray} \tag{ 169 }$$

where for simplicity we assume that Q(r, p) = 0 in r > r₂. The solution of (169) in r > r₂ has the form

$$\begin{eqnarray}&&{f}_{0}(r,p)=C(p){\int }_{{r}_{2}}^{r}\displaystyle \frac{{dr}^{\prime} }{r^{\prime} \kappa (r^{\prime} ,p)}+D(p),\end{eqnarray} \tag{ 170 }$$

where C(p) and D(p) are "integration constants," and we assume that the integral in (170) is finite. Note that

$$\begin{eqnarray}&&\kappa r\displaystyle \frac{\partial {f}_{0}}{\partial r}=C(p)\end{eqnarray} \tag{ 171 }$$

is a first integral of (169). From (171), it follows that

$$\begin{eqnarray}&&C(p)={r}_{2}\kappa ({r}_{2},p)\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}\end{eqnarray} \tag{ 172 }$$

may be identified as the diffusive particle flux at r = r₂.

Letting $r\to \infty$ in solution (170), we obtain

$$\begin{eqnarray}&&D(p)={f}_{0}(\infty ,p)-{r}_{2}\kappa ({r}_{2},p)\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}\left({\int }_{{r}_{2}}^{\infty }\displaystyle \frac{{dr}^{\prime} }{r^{\prime} \kappa (r^{\prime} ,p)}\right),\end{eqnarray} \tag{ 173 }$$

for the unknown function D(p). Using (172) and (173) for C(p) and D(p), we obtain the solution for f₀(r, p) in $r\geqslant {r}_{2}$ as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0}(r,p) & = & {f}_{0}(\infty ,p)-{r}_{2}\kappa ({r}_{2},p)\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}\\ & & \times \,\left({\int }_{r}^{\infty }\displaystyle \frac{{dr}^{\prime} }{r^{\prime} \kappa (r^{\prime} ,p)}\right).\end{array}\end{eqnarray} \tag{ 174 }$$

Setting r = r₂ in (174) gives the boundary condition at r = r₂ as

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{f}_{0}({r}_{2},p)+{r}_{2}\kappa ({r}_{2},p)\displaystyle \frac{{\rm{\partial }}{f}_{0}({r}_{2},p)}{{\rm{\partial }}r}\left({\int }_{{r}_{2}}^{\infty }\displaystyle \frac{{{dr}}^{{\rm{{\prime} }}}}{{r}^{{\rm{{\prime} }}}\kappa ({r}^{{\rm{{\prime} }}},p)}\right)\\ \,=\,{f}_{0}(\infty ,p).\end{array}\end{array}\end{eqnarray} \tag{ 175 }$$

In order to get a better feel for the implications of the boundary condition (175) consider the case where

$$\begin{eqnarray}&&\kappa (r,p)=\displaystyle \frac{{c}^{2}{\tau }_{0}}{3}{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{\alpha }{\left(\displaystyle \frac{r}{{r}_{2}}\right)}^{s},\quad (s\gt 0).\end{eqnarray} \tag{ 176 }$$

Using (176) in (174), we obtain

$$\begin{eqnarray}&&{f}_{0}(r,p)={f}_{0}(\infty ,p)-\displaystyle \frac{{r}_{2}}{s}\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}{\left(\displaystyle \frac{{r}_{2}}{r}\right)}^{s},\end{eqnarray} \tag{ 177 }$$

for f₀ in r > r₂.

The boundary condition (175) then becomes

$$\begin{eqnarray}&&{f}_{0}({r}_{2},p)+\displaystyle \frac{{r}_{2}}{s}\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}={f}_{0}(\infty ,p).\end{eqnarray} \tag{ 178 }$$

The boundary condition (178) is a mixed Dirichlet–Von Neumann boundary condition. For s ≫ 1, (178) is approximately the Dirichlet boundary condition:

$$\begin{eqnarray}&&{f}_{0}({r}_{2},p)\approx {f}_{0}(\infty ,p),\end{eqnarray} \tag{ 179 }$$

which is the free escape boundary condition used in Section 4. In the opposite limit (s → 0), the boundary condition (178) reduces to an approximate Neumann boundary condition:

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{2}}{s}\displaystyle \frac{\partial {f}_{0}({r}_{2},p)}{\partial r}\approx {f}_{0}(\infty ,p).\end{eqnarray} \tag{ 180 }$$

The detailed application of the boundary condition (178) to particle acceleration in radio-jet shear flows lies beyond the scope of the present paper. Intuitively we expect that the boundary condition (175) or (178) will lead to harder particle spectra, because particles crossing r = r₂ into the region r > r₂ may now return back into the region r < r₂ to be further accelerated by the shear flow inside r = r₂.

It is clear that one could also include a source term $Q(r,p)$ on the right-hand side of (169). Similarly, one could include sources of particles in the region 0 < r < r₂ in the problem.

In this appendix we derive formulas for the maximum particle momentum, based on the characteristic classification of the telegrapher transport Equation (63). A similar discussion of the telegrapher equation, including spatial diffusion, was given in Appendix E of Webb et al. (2018a). The basic, space averaged transport Equation (63) may be written in the form

$$\begin{eqnarray}\begin{array}{rcl} & & \left({\tilde{\gamma }}^{2}\tau \displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {t}^{2}}-{D}_{{pp}}\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {p}^{2}}\right)+\gamma \displaystyle \frac{\partial {f}_{0}}{\partial t}-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{D}_{{pp}}\right)\displaystyle \frac{\partial {f}_{0}}{\partial p}\\ & & \quad -{D}_{s}{p}^{2}\displaystyle \frac{\partial {f}_{0}}{\partial p}+{f}_{0}\left(-4{{pD}}_{s}+\displaystyle \frac{1}{{\tau }_{\mathrm{esc}}}\right)=Q.\end{array}\end{eqnarray} \tag{ 181 }$$

In (181)

$$\begin{eqnarray}&&\tilde{\gamma }=\gamma {\left(1-\displaystyle \frac{{\beta }^{2}{\beta }_{p}^{2}}{3}\right)}^{1/2},\end{eqnarray} \tag{ 182 }$$

where $\beta =u/c$ and β_p = v/c are the fluid speed (β) and the particle speed (β_p) normalized to the speed of light. The characteristic manifold of the partial differential Equation (181) is the manifold S(t, p) = const. on which the Cauchy initial value problem for (181) does not have a unique solution. The function S(t, p) (sometimes referred to as the wave eikonal function in JWKB theory) satisfies the first order partial differential equation:

$$\begin{eqnarray}&&G\equiv {\tilde{\gamma }}^{2}\tau {S}_{t}^{2}-{D}_{{pp}}{S}_{p}^{2}=0\end{eqnarray} \tag{ 183 }$$

(e.g., Sneddon 1957; Courant & Hilbert 1989). Note that (183) arises when one tries to calculate higher order derivatives from lower order derivatives by using the partial differential equation and new coordinates ξ(t, p) and η(t, p) to replace t and p as independent coordinates. If the determinant of the transformation from the old coordinates (t, p) to the new coordinates (ξ, η) is zero, then it is not possible to obtain a unique solution (ξ and η can be identified with the possible solutions for S, which in the present case satisfy (183)). Note that ω = −S_t and k = S_p in (183) can be thought of as the frequency ω and wave number k that result from solving the wave equation implicit in (181) which consists of the second derivative terms in (181).

The momentum diffusion coefficient D_pp in (183) due to both cosmic ray viscosity and due to second-order Fermi acceleration arising from Alfvénic turbulence is given by

$$\begin{eqnarray}&&{D}_{{pp}}={\rm{\Gamma }}{p}^{2}\tau +\displaystyle \frac{{p}^{2}{V}_{{\rm{A}}}^{2}}{9\kappa }\equiv {\rm{\Gamma }}{p}^{2}\tau +\displaystyle \frac{{p}^{2}{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}\tau },\end{eqnarray} \tag{ 184 }$$

where β_A = V_A/c is the ratio of the Alfvén speed to the speed of light (we require V_A/c < 1, which means V_A is the relativistic Alfvén speed). The term ${D}_{{pp}}^{A}={p}^{2}{V}_{{\rm{A}}}^{2}/(9\kappa )$ is the second-order Fermi acceleration coefficient due to Alfvén waves (e.g., Skilling 1975).

Using (183)–(184), we write (183) in the form

$$\begin{eqnarray}&&G={\tilde{\gamma }}^{2}\tau \left[{S}_{t}^{2}-\displaystyle \frac{{p}^{2}}{{\tilde{\gamma }}^{2}}\left({\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}{\tau }^{2}}\right){S}_{p}^{2}\right]={\tilde{\gamma }}^{2}\tau {G}_{-}{G}_{+},\end{eqnarray} \tag{ 185 }$$

where

$$\begin{eqnarray}&&{G}_{\pm }={S}_{t}\pm \displaystyle \frac{p}{\tilde{\gamma }}{\left({\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}{\tau }^{2}}\right)}^{1/2}{S}_{p}.\end{eqnarray} \tag{ 186 }$$

Consider the characteristic equation

$$\begin{eqnarray}&&{G}_{+}={S}_{t}+\displaystyle \frac{p}{\tilde{\gamma }}{\left({\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}{\tau }^{2}}\right)}^{1/2}{S}_{p}=0.\end{eqnarray} \tag{ 187 }$$

The characteristics for (187) (the so-called bi-characteristics of (181)) are given by

$$\begin{eqnarray}&&\displaystyle \frac{{dS}}{d\chi }=0,\quad \displaystyle \frac{{dt}}{d\chi }=1,\quad \displaystyle \frac{{dp}}{d\chi }=\displaystyle \frac{p{[{\rm{\Gamma }}+{\beta }_{{\rm{A}}}^{2}/(3{\beta }_{p}^{2}{\tau }^{2})]}^{1/2}}{\tilde{\gamma }},\end{eqnarray} \tag{ 188 }$$

where χ is the affine parameter (or distance parameter) along the characteristic (e.g., Sneddon 1957). From (188) we obtain the differential equation,

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dt}}=\displaystyle \frac{p}{\tilde{\gamma }}{\left[{\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\beta }_{p}^{2}{\tau }_{0}^{2}}{\left(\bar{p}\right)}^{-2\alpha }\right]}^{1/2},\end{eqnarray} \tag{ 189 }$$

for the maximum particle momentum in the telegrapher Equation (181), where we have assumed the mean free time τ(p) has the form

$$\begin{eqnarray}&&\tau ={\tau }_{0}{\bar{p}}^{\alpha },\quad \bar{p}=\displaystyle \frac{p}{{p}_{0}}.\end{eqnarray} \tag{ 190 }$$

At this point, we could integrate (189) directly to obtain the maximum particle momentum p = p_max(t) for particles initially with momentum p = p₀, as predicted by the telegrapher Equation (181). For analytical simplicity, we restrict the analysis below to relativistic particles for which β_p ≈ 1. It is clear that (189) can be integrated for general β_p if necessary by numerical integration. In the analysis below we consider only the case of relativistic particles in which we set β_p = 1. We prefer here to illustrate the solutions of (189) for relativistic particles (β_p = 1) for the special cases (i) β_A = 0 and ${\rm{\Gamma }}\ne 0$ (pure viscous acceleration case) and (ii) ${\beta }_{{\rm{A}}}\ne 0$ and Γ = 0 (pure second-order Fermi case). In case (i)

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dt}}=\displaystyle \frac{\sqrt{{\rm{\Gamma }}}}{\tilde{\gamma }}p.\end{eqnarray} \tag{ 191 }$$

Integration of (191) with respect to t, with initial condition p = p₀ at time t = 0, gives the solution

$$\begin{eqnarray}&&p={p}_{0}\exp \left(\displaystyle \frac{\sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right),\end{eqnarray} \tag{ 192 }$$

for the maximum particle momentum due to pure shear acceleration alone. Similarly, in case (ii), (189) reduces to the equation

$$\begin{eqnarray}&&\displaystyle \frac{d\bar{p}}{{dt}}=\displaystyle \frac{{\beta }_{{\rm{A}}}{\left(\bar{p}\right)}^{1-\alpha }}{\sqrt{3}\tilde{\gamma }{\tau }_{0}}.\end{eqnarray} \tag{ 193 }$$

Integration of (193) with $\bar{p}=1$ at time t = 0 (i.e., p = p₀ at t = 0) gives the solution

$$\begin{eqnarray}&&p(t)={p}_{0}{\left[1+\displaystyle \frac{\alpha {\beta }_{{\rm{A}}}}{\sqrt{3}\tilde{\gamma }}\displaystyle \frac{t}{{\tau }_{0}}\right]}^{1/\alpha },\end{eqnarray} \tag{ 194 }$$

for the maximum particle momentum, where we assume $\alpha \ne 0$. For the case α = 0, we obtain the solution

$$\begin{eqnarray}&&p(t)={p}_{0}\exp \left(\displaystyle \frac{{\beta }_{{\rm{A}}}t}{\sqrt{3}\tilde{\gamma }{\tau }_{0}}\right).\end{eqnarray} \tag{ 195 }$$

Equations (194) and (195) give the maximum particle momentum for the case of pure second-order Fermi acceleration (i.e., case (ii)).

Comparing the solutions (192)–(195), we see that the solutions for second-order Fermi acceleration and no shear (Equations (194) and (195)) depend on the scattering time τ₀, but the pure shear acceleration solution (192) does not depend on τ₀.

Finally, we note that integration of (189) for the general case where ${\beta }_{{\rm{A}}}\ne 0$ and ${\rm{\Gamma }}\ne 0$ and β_p = 1, satisfying p = p₀ at t = 0 may be written in the form

$$\begin{eqnarray}&&{\bar{p}}^{\alpha }={\left({A}^{2}+1\right)}^{1/2}\sinh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)+\cosh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right),\end{eqnarray} \tag{ 196 }$$

where

$$\begin{eqnarray}&&A=\displaystyle \frac{{\beta }_{{\rm{A}}}}{\sqrt{3{\rm{\Gamma }}}{\tau }_{0}}.\end{eqnarray} \tag{ 197 }$$

From (196) we obtain

$$\begin{eqnarray}\begin{array}{rcl}p(t) & = & {p}_{0}\left\{\cosh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)+{\left({A}^{2}+1\right)}^{1/2}\right.\\ & & \times {\left.\sinh \left(\displaystyle \frac{\alpha \sqrt{{\rm{\Gamma }}}t}{\tilde{\gamma }}\right)\right\}}^{1/\alpha }\end{array}\end{eqnarray} \tag{ 198 }$$

as the general solution for the maximum particle momentum p(t) allowed by the telegrapher Equation (181), including both second-order Fermi acceleration (${\beta }_{{\rm{A}}}\ne 0$) and viscous shear acceleration (${\rm{\Gamma }}\ne 0$) where $\alpha \ne 0$. For the case α = 0, (189) simplifies to

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dt}}=\displaystyle \frac{p\sqrt{{{\rm{\Gamma }}}_{* }}}{\tilde{\gamma }}\quad \mathrm{where}\quad {{\rm{\Gamma }}}_{* }=\left({\rm{\Gamma }}+\displaystyle \frac{{\beta }_{{\rm{A}}}^{2}}{3{\tau }_{0}^{2}}\right).\end{eqnarray} \tag{ 199 }$$

Integration of (199) with p = p₀ at t = 0 gives the solution

$$\begin{eqnarray}&&p={p}_{0}\exp \left(\displaystyle \frac{\sqrt{{{\rm{\Gamma }}}_{* }}t}{\tilde{\gamma }}\right),\end{eqnarray} \tag{ 200 }$$

where Γ_* in (199) contains both shear acceleration (Γ) and second-order Fermi (${\beta }_{{\rm{A}}}^{2}/(3{\tau }_{0}^{2})$) terms.

D.1. An Anomaly for $\langle {\rm{\Delta }}p/{\rm{\Delta }}t{\rangle }_{\mathrm{sh}}$

If one uses the mean acceleration rate for shear acceleration in (87) and integrates the differential equation

$$\begin{eqnarray}&&{\left\langle \displaystyle \frac{{\rm{\Delta }}p}{{\rm{\Delta }}t}\right\rangle }_{\mathrm{sh}}=\displaystyle \frac{{dp}}{{dt}}=(4+\alpha ){\rm{\Gamma }}p\tau \equiv (4+\alpha ){\rm{\Gamma }}{p}_{0}{\tau }_{0}{\bar{p}}^{\alpha +1},\end{eqnarray} \tag{ 201 }$$

with p = p₀ at t = 0, one obtains the solution for p as

$$\begin{eqnarray}&&p(t)={p}_{0}{\left[1-\alpha (4+\alpha ){\rm{\Gamma }}{\tau }_{0}t\right]}^{-1/\alpha },\end{eqnarray} \tag{ 202 }$$

Alternatively, the solution (202) can be expressed in the form

$$\begin{eqnarray}&&t=\displaystyle \frac{1}{\alpha (4+\alpha ){\rm{\Gamma }}{\tau }_{0}}\left[1-{\left(\displaystyle \frac{p}{{p}_{0}}\right)}^{-\alpha }\right].\end{eqnarray} \tag{ 203 }$$

From (202) or (203), we find

$$\begin{eqnarray}&&p\to \infty \quad \mathrm{as}\quad t\to \displaystyle \frac{1}{\alpha (4+\alpha ){\rm{\Gamma }}{\tau }_{0}}.\end{eqnarray} \tag{ 204 }$$

The result (204) is clearly non-physical, as it implies that particles can be accelerated to an infinite momentum in a finite time. This non-causal behavior is a result of neglecting the finite inertia of the cosmic rays (i.e., the telegrapher second-order time derivative term in (181)).

One can obtain solutions of the characteristic Equation (185) corresponding to the solution branch G₋ = 0. This in general gives a lower momentum p_min(t) implicit in the solution of (181). However, it turns out in the case of second-order Fermi acceleration that the minimum momentum p_min can be negative, which of course is not physically acceptable. Also note that the coefficient of f_0,tt in (181) tends to zero as p → 0, suggesting that the f_0,tt term should be neglected at small p (assuming τ ∝ p^α with α > 0).

In this appendix we investigate the steady-state leaky box model (63) for particle acceleration due to cosmic ray viscosity in the presence of synchrotron losses and particle escape. We concentrate on the case where ${\tau }_{\mathrm{esc}}\to \infty$ and Q = 0.

The steady-state transport Equation (63) with Q = 0 has the form

$$\begin{eqnarray}&&-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({p}^{2}{D}_{{pp}}\displaystyle \frac{\partial {f}_{0}}{\partial p}\right)-\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left({D}_{s}{p}^{4}{f}_{0}\right)+\displaystyle \frac{{f}_{0}}{{\tau }_{\mathrm{esc}}}=0.\end{eqnarray} \tag{ 205 }$$

We consider the simplest case where the escape time ${\tau }_{\mathrm{esc}}\to \infty$ (i.e., there is no particle escape from the box). Here

$$\begin{eqnarray}&&{D}_{{pp}}={p}^{2}{\rm{\Gamma }}\tau ={\rm{\Gamma }}{p}^{2}{\tau }_{0}{\bar{p}}^{\alpha }\quad \mathrm{where}\quad \bar{p}=\displaystyle \frac{p}{{p}_{0}}.\end{eqnarray} \tag{ 206 }$$

Equation (205) with ${\tau }_{\mathrm{esc}}\to \infty$ integrates to give the first integral,

$$\begin{eqnarray}&&{\rm{\Gamma }}{p}^{4}\tau \displaystyle \frac{\partial {f}_{0}}{\partial p}+{D}_{s}{p}^{4}{f}_{0}={C}_{0},\end{eqnarray} \tag{ 207 }$$

where C₀ is the integration constant, and where for the sake of simplicity we take Γ to be a constant. Equation (207) can be written in the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\bar{f}}_{0}}{\partial \bar{p}}+\displaystyle \frac{{D}_{s}{p}_{0}}{{\rm{\Gamma }}{\tau }_{0}}{\left(\bar{p}\right)}^{-\alpha }{\bar{f}}_{0}(\bar{p})=\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}{\left(\bar{p}\right)}^{-\alpha -4},\end{eqnarray} \tag{ 208 }$$

where we use the notation ${\bar{f}}_{0}(\bar{p})={f}_{0}(p)$ to denote the function f₀(p) in terms of the normalized momentum $\bar{p}=p/{p}_{0}$.

The integrating factor for Equation (208) depends on whether α = 1 or $\alpha \ne 1$. For $\alpha \ne 1$, the integrating factor $I(\bar{p})$ is

$$\begin{eqnarray}&&I(\bar{p})=\exp \left(\int \chi {\left(\bar{p}\right)}^{-\alpha }d\bar{p}\right)=\exp \left(\chi \displaystyle \frac{{\left(\bar{p}\right)}^{(1-\alpha )}}{(1-\alpha )}\right),\end{eqnarray} \tag{ 209 }$$

where

$$\begin{eqnarray}&&\chi =\displaystyle \frac{{D}_{s}{p}_{0}}{{\rm{\Gamma }}{\tau }_{0}}\equiv \displaystyle \frac{{D}_{s}{p}_{0}^{4}}{{\rm{\Gamma }}{p}_{0}^{3}{\tau }_{0}}\end{eqnarray} \tag{ 210 }$$

is a dimensionless parameter that measures the relative importance of synchrotron losses to shear acceleration due to cosmic ray viscosity. For α = 1 (Bohm diffusion case), we obtain

$$\begin{eqnarray}&&I(\bar{p})=\exp \left(\chi \mathrm{ln}(\bar{p})\right)={\bar{p}}^{\chi },\end{eqnarray} \tag{ 211 }$$

for the integrating factor $I(\bar{p})$.

The general solution of (208) using the integrating factor (209) is

$$\begin{eqnarray}&&{\bar{f}}_{0}(\bar{p})=\displaystyle \frac{{C}_{1}}{I(\bar{p})}+\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}I(\bar{p})}{\int }_{1}^{\bar{p}}I(\bar{p}^{\prime} ){\left(\bar{p}^{\prime} \right)}^{-\alpha -4}\ d\bar{p}^{\prime} ,\end{eqnarray} \tag{ 212 }$$

where C₁ is a second integration constant. In (208)–(212) we identify the constants C₁ and C₀ for $\alpha \ne 1$ as

$$\begin{eqnarray}\begin{array}{rcl}{C}_{1} & = & \exp \left(\displaystyle \frac{\chi }{1-\alpha }\right){\bar{f}}_{0}({\bar{p}}_{0}),\\ {C}_{0} & = & {\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}\left(\displaystyle \frac{\partial {\bar{f}}_{0}({\bar{p}}_{0})}{\partial {\bar{p}}_{0}}+\chi {\bar{f}}_{0}({\bar{p}}_{0})\right).\end{array}\end{eqnarray} \tag{ 213 }$$

For the case α = 1, C₁ = f₀(p₀) and C₀ are given by (213). Thus, in general, ${C}_{0}\ne 0$ and ${C}_{1}\ne 0$.

For the case of Bohm diffusion (α = 1), the solution (212) reduces to the simple form

$$\begin{eqnarray}&&{\bar{f}}_{0}=\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}\displaystyle \frac{{\left(\bar{p}\right)}^{-4}}{(\chi -4)}+{\tilde{C}}_{1}{\left(\bar{p}\right)}^{-\chi },\end{eqnarray} \tag{ 214 }$$

where

$$\begin{eqnarray}&&{\tilde{C}}_{1}={C}_{1}-\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{p}_{0}^{3}{\tau }_{0}(\chi -4)}.\end{eqnarray} \tag{ 215 }$$

The first term in (214) represents pure synchrotron losses. It dominates the solution for p ≫ p₀ if χ > 4. The second term in (214) represents the competition between synchrotron losses and particle acceleration due to cosmic ray viscosity. It dominates the solution at large momenta p ≫ p₀ in the case where 0 < χ < 4 (this is the case of dominant particle acceleration due to cosmic ray viscosity), resulting in ${f}_{0}\propto {\bar{p}}^{-\chi }$ as $p\to \infty$.

For the case 0 < α < 1, the solution (212) for p ≫ p₀ can be approximated by

$$\begin{eqnarray}&&{\bar{f}}_{0}\sim {C}_{1}\exp \left(-\displaystyle \frac{\chi {\left(\bar{p}\right)}^{1-\alpha }}{1-\alpha }\right)+\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}\chi }{\left(\bar{p}\right)}^{-4}\left[1+O\left(\displaystyle \frac{1}{w\chi }\right)\right],\end{eqnarray} \tag{ 216 }$$

where

$$\begin{eqnarray}&&w=\displaystyle \frac{{\left(\bar{p}\right)}^{1-\alpha }}{1-\alpha }.\end{eqnarray} \tag{ 217 }$$

Equation (216) shows that for $0\lt \alpha \lt 1$ the pure synchrotron spectrum dominates the shear viscosity acceleration process as $\bar{p}\to \infty$. This result agrees with the results of (90), where synchrotron losses beat shear acceleration gains at $p\gt {p}_{\mathrm{sh}}^{* }$. The result (216) follows from an analysis of the second term in (212), namely

$$\begin{eqnarray}&&{T}_{2}=\displaystyle \frac{1}{I(\bar{p})}{\int }_{1}^{\bar{p}}\ I(\bar{p}^{\prime} ){\left(\bar{p}^{\prime} \right)}^{-\alpha -4}\ {dp}^{\prime} \end{eqnarray} \tag{ 218 }$$

and by using Laplace's asymptotic method to evaluate the integral (218) for large values of the parameter w. To evaluate (218) for large w, note that $I(\bar{p})=\exp (\chi w)$, so that

$$\begin{eqnarray}\begin{array}{rcl}{T}_{2}(\bar{p}) & = & {\int }_{w^{\prime} =1/(1-\alpha )}^{w^{\prime} =w}\exp \left[-\chi (w-w^{\prime} )\right]\\ & & \times {\left[(1-\alpha )w^{\prime} \right]}^{-4/(1-\alpha )}\ {dw}^{\prime} ,\end{array}\end{eqnarray} \tag{ 219 }$$

and then consecutively use the integration variables

$$\begin{eqnarray}&&\zeta =w-w^{\prime} ,\quad \hat{\zeta }=w\zeta ,\end{eqnarray} \tag{ 220 }$$

followed by the use of Laplaces's method using integration by parts to asymptotically evaluate ${T}_{2}(\bar{p})$ for the case where χw is a large parameter (we omit the details). The net result is the asymptotic solution (216) for f₀ for large p.

Similarly, for the case α > 1, the solution (212) for f₀ may be written in the form

$$\begin{eqnarray}&&{f}_{0}=\displaystyle \frac{1}{I(\bar{p})}\left\{\left[{C}_{1}+\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}{I}_{2}(0,1)\right]-\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}{I}_{2}(0,\tilde{w})\right\},\end{eqnarray} \tag{ 221 }$$

where

$$\begin{eqnarray}&&\tilde{w}={\left(\bar{p}\right)}^{1-\alpha }\equiv (1-\alpha )w.\end{eqnarray} \tag{ 222 }$$

Here

$$\begin{eqnarray}\begin{array}{rcl}{I}_{2}(0,\tilde{w}) & = & {\displaystyle \int }_{0}^{\tilde{w}}\displaystyle \frac{{dw}^{\prime} }{\alpha -1}\exp \left(-\displaystyle \frac{\chi w^{\prime} }{\alpha -1}\right){w}^{{\prime} 4/(\alpha -1)}\ {dw}^{\prime} ,\\ {I}_{2}(0,1) & = & {\displaystyle \int }_{0}^{1}\displaystyle \frac{{dw}^{\prime} }{\alpha -1}\exp \left(-\displaystyle \frac{\chi w^{\prime} }{\alpha -1}\right){w}^{{\prime} 4/(\alpha -1)}{dw}^{\prime} .\end{array}\end{eqnarray} \tag{ 223 }$$

The integrals (223) can be expressed in terms of the incomplete Gamma function γ (a, x):

$$\begin{eqnarray}&&\gamma (a,x)={\int }_{0}^{x}\exp (-t){t}^{a-1}\ {dt}\end{eqnarray} \tag{ 224 }$$

(Abramowitz & Stegun 1965, p. 260, formula 6.5.2). We find

$$\begin{eqnarray}\begin{array}{rcl}{I}_{2}(0,\tilde{w}) & = & {\left(\displaystyle \frac{\alpha -1}{\chi }\right)}^{(\alpha +3)/(\alpha -1)}\displaystyle \frac{1}{\alpha -1}\gamma \left(\displaystyle \frac{\alpha +3}{\alpha -1},\displaystyle \frac{\chi \tilde{w}}{\alpha -1}\right),\\ {I}_{2}(0,1) & = & {\left(\displaystyle \frac{\alpha -1}{\chi }\right)}^{(\alpha +3)/(\alpha -1)}\displaystyle \frac{1}{\alpha -1}\gamma \left(\displaystyle \frac{\alpha +3}{\alpha -1},\displaystyle \frac{\chi }{\alpha -1}\right).\end{array}\end{eqnarray} \tag{ 225 }$$

Note from (209), the integrating factor $I(\bar{p})\to 1$ as $\bar{p}\to \infty$. Because we require ${f}_{0}\to 0$ as $\bar{p}\to \infty$, we set

$$\begin{eqnarray}&&{C}_{1}+\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}{I}_{2}(0,1)=0\end{eqnarray} \tag{ 226 }$$

as a relationship between C₁ and C₀ required for the physical solution. Using (226), the solution for ${f}_{0}(\bar{p})$ in (221) reduces to

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0} & = & -\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}\exp \left(\displaystyle \frac{\chi {\bar{p}}^{1-\alpha }}{\alpha -1}\right)\left[{\left(\displaystyle \frac{\alpha -1}{\chi }\right)}^{(\alpha +3)/(\alpha -1)}\right.\\ & & \times \left.\displaystyle \frac{1}{\alpha -1}\gamma \left(\displaystyle \frac{\alpha +3}{\alpha -1},\displaystyle \frac{\chi \tilde{w}}{\alpha -1}\right)\right].\end{array}\end{eqnarray} \tag{ 227 }$$

Using the expansion

$$\begin{eqnarray}\begin{array}{rcl}\gamma (a,x) & = & {\rm{\Gamma }}(a){x}^{a}{\gamma }^{* }(a,x)={\rm{\Gamma }}(a){x}^{a}\exp (-x)\\ & & \times \,\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{x}^{n}}{{\rm{\Gamma }}(a+n+1)}\end{array}\end{eqnarray} \tag{ 228 }$$

(Abramowitz & Stegun 1965, p. 260, formula 6.5.29), the solution (227) for $\bar{p}\to \infty$ reduces to the form

$$\begin{eqnarray}\begin{array}{rcl}{f}_{0} & = & -\displaystyle \frac{{C}_{0}}{{\rm{\Gamma }}{\tau }_{0}{p}_{0}^{3}}\displaystyle \frac{1}{\alpha +3}{\left(\bar{p}\right)}^{-\alpha -3}\\ & & \times \left(1+\displaystyle \frac{y}{\delta +1}+\displaystyle \frac{{y}^{2}}{(\delta +1)(\delta +2)}+\ldots \right),\end{array}\end{eqnarray} \tag{ 229 }$$

where

$$\begin{eqnarray}&&y=\displaystyle \frac{\chi {\left(\bar{p}\right)}^{1-\alpha }}{\alpha -1},\quad \delta =\displaystyle \frac{\alpha +3}{\alpha -1}.\end{eqnarray} \tag{ 230 }$$

The asymptotic formula for f₀ as $\bar{p}\to \infty$ in (229) is that expected from particle acceleration in shear flows in the absence of losses (e.g., Berezhko & Krymsky 1981; Rieger & Duffy 2006). This agrees with the analysis of (90), in which particle acceleration by shear dominates the synchrotron losses at $p\gt {p}_{\mathrm{sh}}^{* }$ for the case α > 1.

If we assume a finite escape time τ_esc(p) in (205), then we expect f₀ → 0 as $p\to \infty$ for τ_esc, a sufficiently fast increasing function of p. If one assumes an escape time

$$\begin{eqnarray}&&{\tau }_{\mathrm{esc}}={\tau }_{e0}{\bar{p}}^{b}\end{eqnarray} \tag{ 231 }$$

in the leaky box transport Equation (205), and assuming Γ can be taken as a constant, (205) reduces to the second-order ordinary differential equation

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{d}^{2}{\bar{f}}_{0}}{d{\bar{p}}^{2}}+\left[\displaystyle \frac{4+\alpha }{\bar{p}}+\chi {\left(\bar{p}\right)}^{-\alpha }\right]\displaystyle \frac{d{\bar{f}}_{0}}{d\bar{p}}+{\left(\bar{p}\right)}^{-1-\alpha }\\ \quad \times \left[4\chi -\displaystyle \frac{1}{{\rm{\Gamma }}{\tau }_{0}{\tau }_{e0}}{\left(\bar{p}\right)}^{-b-1}\right]{\bar{f}}_{0}=0.\end{array}\end{eqnarray} \tag{ 232 }$$

For the case ${\tau }_{e0}\to \infty$, the transformation of variables,

$$\begin{eqnarray}&&{f}_{0}(\bar{p})=F(y),\quad y=\displaystyle \frac{\chi {\left(\bar{p}\right)}^{1-\alpha }}{\alpha -1},\end{eqnarray} \tag{ 233 }$$

transforms the ordinary differential Equation (232) to the confluent hypergeometric equation

$$\begin{eqnarray}&&y\displaystyle \frac{{d}^{2}F}{{{dy}}^{2}}+\left(\displaystyle \frac{4}{1-\alpha }-y\right)\displaystyle \frac{{dF}}{{dy}}-\displaystyle \frac{4}{1-\alpha }F=0,\end{eqnarray} \tag{ 234 }$$

which has the general solution

$$\begin{eqnarray}\begin{array}{rcl}F & = & {a}_{1}M(A,A,y)+{a}_{2}U(A,A,y)\equiv {a}_{1}\exp (y)\\ & & +{a}_{2}U(A,A,y)\quad \mathrm{where}\quad A=\displaystyle \frac{4}{1-\alpha }.\end{array}\end{eqnarray} \tag{ 235 }$$

Here a₁ and a₂ are arbitrary constants and M(a, b, y) and U(a, b, y) are independent solutions of the Kummer's confluent hypergeometric equation (e.g., Abramowitz & Stegun 1965, Ch. 13, p.504). The solutions (235) could have been used to carry out the asymptotics of the solutions for $\bar{p}\to \infty$ for ${\tau }_{e0}\to \infty$. Further analysis of (232) for a finite escape time τ_e0 is an interesting problem, which lies beyond the scope of the present paper.

Here we correct the formula for r_gp of (99) of Webb et al. (2018a), and correct other Equations in (99) et seq. of that paper. The corrected r_gp is

$$\begin{eqnarray}&&{r}_{\mathrm{gp}}=3.3\times {10}^{20}\left(\displaystyle \frac{E}{{10}^{18}\,\mathrm{eV}}\right){\left(\displaystyle \frac{B}{10\mu {\rm{G}}}\right)}^{-1}\mathrm{cm}.\end{eqnarray} \tag{ 236 }$$

Condition ΔL > r_gp/η, for particle confinement by the jet in Webb et al. (2018a) (103), reduces to

$$\begin{eqnarray}&&{\rm{\Delta }}L\gt 3.3\times {10}^{20}{\eta }^{-1}\left(\displaystyle \frac{E}{{10}^{18}\,\mathrm{eV}}\right){\left(\displaystyle \frac{B}{10\mu {\rm{G}}}\right)}^{-1}\mathrm{cm}.\end{eqnarray} \tag{ 237 }$$

Webb et al. (2018a) (107) for the acceleration time t_acc becomes

$$\begin{eqnarray}\begin{array}{rcl}{t}_{\mathrm{acc}} & = & 1.67\times {10}^{-17}{\left[{\epsilon }_{\perp }N(5/3)\right]}^{-1}\chi {\eta }^{-2/3}{\left({\rm{\Delta }}L\right)}^{4/3}\\ & & \times {\left(\displaystyle \frac{E}{{10}^{18}\mathrm{eV}}\right)}^{-1/3}{\left(\displaystyle \frac{B}{10\mu {\rm{G}}}\right)}^{1/3}{\left({\beta }_{j}-{\beta }_{2}\right)}^{-2}{\gamma }_{j}^{-3}\ {\rm{s}}.\end{array}\end{eqnarray} \tag{ 238 }$$

Formula (106) for the synchrotron loss time remains the same as in Webb et al. (2018a). Condition t_acc < t_synch for particles to be accelerated despite synchrotron losses reduces to the condition

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}L & \lt & 2.1\times {10}^{27}{\left[{\epsilon }_{\perp }N(5/3)\right]}^{3/4}{\eta }^{1/2}{\chi }^{-3/4}{\left(\displaystyle \frac{E}{{10}^{18}\mathrm{eV}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{B}{10\mu {\rm{G}}}\right)}^{-7/4}{\gamma }_{j}^{9/4}{\left({\beta }_{j}-{\beta }_{2}\right)}^{3/2}\ \mathrm{cm}.\end{array}\end{eqnarray} \tag{ 239 }$$

Assuming that ${\rm{\Delta }}L={\rho }_{w}d$ (ρ_w = 0.01), the condition t_acc < t_dyn = d/(cβ_j) (i.e., the corrected (109) of Webb et al. 2018a) gives the constraint

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}L & \lt & 7.95\times {10}^{24}{\left[{\epsilon }_{\perp }N(5/3)\right]}^{3}{\chi }^{-3}{\eta }^{2}\left(\displaystyle \frac{E}{{10}^{18}\,\mathrm{eV}}\right){\left(\displaystyle \frac{B}{10\mu {\rm{G}}}\right)}^{-1}\\ & & \times \,{\beta }_{j}^{3}{\gamma }_{j}^{9}{\left(\displaystyle \frac{{\rho }_{w}}{0.01}\right)}^{-3}\ \mathrm{cm}.\end{array}\end{eqnarray} \tag{ 240 }$$

These equations delimit the ${\rm{\Delta }}L$ and B regions in which the E = 1EeV particles are confined during acceleration in radio jet and correct (99) et seq. of Webb et al. (2018a).

Figures 10 and 11 (the corrected Figure 10 and 11 of Webb et al. 2018a) illustrate the three constraints (i) ΔL > r_gp (237) with η = 1, (ii) the acceleration time t_acc < t_synch (239), and (iii) the acceleration time t_acc < t_dyn (240). The allowed ΔL and B(μG) regions are shown by the gray-shaded regions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size