Cosmic Ray Acceleration in Hydromagnetic Flux Tubes

Anthony Bell; James Matthews; Katherine Blundell; Anabella Araudo

arXiv:1906.02508·astro-ph.HE·June 26, 2019

Cosmic Ray Acceleration in Hydromagnetic Flux Tubes

Anthony Bell, James Matthews, Katherine Blundell, Anabella Araudo

PDF

TL;DR

This paper proposes that hydromagnetic flux tubes in radio galaxy lobes can accelerate cosmic rays to ultra-high energies through Fermi acceleration mechanisms, even without strong turbulence.

Contribution

It introduces a novel environment for cosmic ray acceleration and demonstrates that weak shocks can effectively accelerate particles to ultra-high energies.

Findings

01

CR can reach Hillas energies in flux tubes

02

First-order Fermi acceleration is equivalent to second-order in this context

03

Acceleration occurs even with weak turbulence

Abstract

We find that hydromagnetic flux tubes in back-flows in the lobes of radio galaxies offer a suitable environment for the acceleration of cosmic rays (CR) to ultra-high energies. We show that CR can reach the Hillas (1984) energy even if the magnetised turbulence in the flux tube is not sufficiently strong for Bohm diffusion to apply. First-order Fermi acceleration by successive weak shocks in a hydromagnetic flux tube is shown to be equivalent to second-order Fermi acceleration by strong turbulence.

Equations66

\frac{\partial f}{\partial t} + \frac{\partial ( f u )}{\partial z} - \frac{\partial}{\partial z} (D_{∥} \frac{\partial f}{\partial z}) - \frac{\partial}{\partial x} (D_{⊥} \frac{\partial f}{\partial x}) - \frac{1}{3} \frac{\partial u}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f )}{\partial p} = 0

\frac{\partial f}{\partial t} + \frac{\partial ( f u )}{\partial z} - \frac{\partial}{\partial z} (D_{∥} \frac{\partial f}{\partial z}) - \frac{\partial}{\partial x} (D_{⊥} \frac{\partial f}{\partial x}) - \frac{1}{3} \frac{\partial u}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f )}{\partial p} = 0

f_{\pm} (z, x, p) = m = 0 \sum \infty f_{s m} (p) exp (k_{m \pm} z) cos (\frac{( 2 m + 1 ) π x}{l})

f_{\pm} (z, x, p) = m = 0 \sum \infty f_{s m} (p) exp (k_{m \pm} z) cos (\frac{( 2 m + 1 ) π x}{l})

f_{s m} (p) = F_{m} p^{- \frac{3 u _{-}}{u _{-} - u _{+}}} [\frac{1}{2} + \frac{1}{4} + \frac{k _{m}^{2} D _{∥ -} D _{⊥ -}}{u _{-}^{2}}]^{3 u_{-} / 2 (u_{-} - u_{+})} [\frac{1}{2} + \frac{1}{4} + \frac{k _{m}^{2} D _{∥ +} D _{⊥ +}}{u _{+}^{2}}]^{3 u_{+} / 2 (u_{-} - u_{+})} \times

f_{s m} (p) = F_{m} p^{- \frac{3 u _{-}}{u _{-} - u _{+}}} [\frac{1}{2} + \frac{1}{4} + \frac{k _{m}^{2} D _{∥ -} D _{⊥ -}}{u _{-}^{2}}]^{3 u_{-} / 2 (u_{-} - u_{+})} [\frac{1}{2} + \frac{1}{4} + \frac{k _{m}^{2} D _{∥ +} D _{⊥ +}}{u _{+}^{2}}]^{3 u_{+} / 2 (u_{-} - u_{+})} \times

exp [\frac{3}{2} \frac{u _{-} + u _{+}}{u _{-} - u _{+}}] exp - \frac{3}{2} \frac{u _{-}^{2} + 4 k _{m}^{2} D _{∥ -} D _{⊥ -} + u _{+}^{2} + 4 k _{m}^{2} D _{∥ +} D _{⊥ +}}{u _{-} - u _{+}}

exp [\frac{3}{2} \frac{u _{-} + u _{+}}{u _{-} - u _{+}}] exp - \frac{3}{2} \frac{u _{-}^{2} + 4 k _{m}^{2} D _{∥ -} D _{⊥ -} + u _{+}^{2} + 4 k _{m}^{2} D _{∥ +} D _{⊥ +}}{u _{-} - u _{+}}

\frac{\partial N}{\partial t} = \int_{- \infty}^{\infty} d z \int_{- l /2}^{l /2} d x \int_{p}^{\infty} 4 π p^{2} d p [\frac{1}{3} \frac{\partial u}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f )}{\partial p}] = \frac{4 π}{3} (u_{-} - u_{+}) p^{3} \int_{- l /2}^{l /2} f_{s} d x

\frac{\partial N}{\partial t} = \int_{- \infty}^{\infty} d z \int_{- l /2}^{l /2} d x \int_{p}^{\infty} 4 π p^{2} d p [\frac{1}{3} \frac{\partial u}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f )}{\partial p}] = \frac{4 π}{3} (u_{-} - u_{+}) p^{3} \int_{- l /2}^{l /2} f_{s} d x

giving \frac{\partial N}{\partial t} = \frac{8 ( u _{-} - u _{+} ) l}{3} m = 0 \sum \infty \frac{( - 1 ) ^{m}}{2 m + 1} p^{3} f_{s m} .

giving \frac{\partial N}{\partial t} = \frac{8 ( u _{-} - u _{+} ) l}{3} m = 0 \sum \infty \frac{( - 1 ) ^{m}}{2 m + 1} p^{3} f_{s m} .

\frac{\partial n}{\partial t} = m = 0 \sum \infty \frac{4 l ( - 1 ) ^{m}}{2 m + 1} [u_{-}^{2} + 4 k_{m}^{2} D_{∥ -} D_{⊥ -} - u_{-} + u_{+}^{2} + 4 k_{m}^{2} D_{∥ +} D_{⊥ +} + u_{+}] p^{2} f_{s m} .

\frac{\partial n}{\partial t} = m = 0 \sum \infty \frac{4 l ( - 1 ) ^{m}}{2 m + 1} [u_{-}^{2} + 4 k_{m}^{2} D_{∥ -} D_{⊥ -} - u_{-} + u_{+}^{2} + 4 k_{m}^{2} D_{∥ +} D_{⊥ +} + u_{+}] p^{2} f_{s m} .

D_{B o hm} = \frac{r _{g} c}{3} = \frac{p c}{3 e B _{∥}}

D_{B o hm} = \frac{r _{g} c}{3} = \frac{p c}{3 e B _{∥}}

f_{s m} = (- 1)^{m} \frac{4}{π} \frac{F _{s} p ^{- 4}}{2 m + 1} (\frac{1}{2} + \frac{1}{4} + \frac{χ _{m}^{2}}{4})^{2} (\frac{1}{2} + \frac{1}{4} + 4 χ_{m}^{2})^{1/2} exp (\frac{5}{2} - 4 + 4 χ_{m}^{2} - 1/4 + 4 χ_{m}^{2})

f_{s m} = (- 1)^{m} \frac{4}{π} \frac{F _{s} p ^{- 4}}{2 m + 1} (\frac{1}{2} + \frac{1}{4} + \frac{χ _{m}^{2}}{4})^{2} (\frac{1}{2} + \frac{1}{4} + 4 χ_{m}^{2})^{1/2} exp (\frac{5}{2} - 4 + 4 χ_{m}^{2} - 1/4 + 4 χ_{m}^{2})

{\rm and}\hskip 14.22636pt\frac{\partial n}{\partial t}=u_{s}l\sum_{m=0}^{\infty}\frac{(-1)^{m}}{2m+1}\big{(}\sqrt{16+16\chi_{m}^{2}}+\sqrt{1+16\chi_{m}^{2}}-3\big{)}p^{2}f_{sm}\hskip 19.91684pt{\rm where}\ \ \chi_{m}=(2m+1)\frac{2\pi D_{s}}{u_{s}l}\ .

{\rm and}\hskip 14.22636pt\frac{\partial n}{\partial t}=u_{s}l\sum_{m=0}^{\infty}\frac{(-1)^{m}}{2m+1}\big{(}\sqrt{16+16\chi_{m}^{2}}+\sqrt{1+16\chi_{m}^{2}}-3\big{)}p^{2}f_{sm}\hskip 19.91684pt{\rm where}\ \ \chi_{m}=(2m+1)\frac{2\pi D_{s}}{u_{s}l}\ .

T_{ma x} \approx 0.6 u_{s} B l .

T_{ma x} \approx 0.6 u_{s} B l .

u (z) = u_{0} + u_{1} (z) .

u (z) = u_{0} + u_{1} (z) .

f (z, x, p) = [f_{0} (p) + f_{1} (z, p)] cos (k_{⊥} x)

f (z, x, p) = [f_{0} (p) + f_{1} (z, p)] cos (k_{⊥} x)

u_{0} \frac{\partial f _{1}}{\partial z} - D_{∥} \frac{\partial ^{2} f _{1}}{\partial z ^{2}} + k_{⊥}^{2} D_{⊥} f_{1} = \frac{1}{3} \frac{\partial u _{1}}{\partial z} p \frac{\partial f _{0}}{\partial p}

u_{0} \frac{\partial f _{1}}{\partial z} - D_{∥} \frac{\partial ^{2} f _{1}}{\partial z ^{2}} + k_{⊥}^{2} D_{⊥} f_{1} = \frac{1}{3} \frac{\partial u _{1}}{\partial z} p \frac{\partial f _{0}}{\partial p}

f_{1} (z, p) = \frac{1}{3 u _{0} ξ} p \frac{\partial f _{0}}{\partial p} [\int_{- \infty}^{z} \frac{\partial u _{1}}{\partial z ^{'}} e^{b (z^{'} - z)} d z^{'} + \int_{z}^{\infty} \frac{\partial u _{1}}{\partial z ^{'}} e^{a (z^{'} - z)} d z^{'}]

f_{1} (z, p) = \frac{1}{3 u _{0} ξ} p \frac{\partial f _{0}}{\partial p} [\int_{- \infty}^{z} \frac{\partial u _{1}}{\partial z ^{'}} e^{b (z^{'} - z)} d z^{'} + \int_{z}^{\infty} \frac{\partial u _{1}}{\partial z ^{'}} e^{a (z^{'} - z)} d z^{'}]

where a = \frac{u _{0} ( - 1 - ξ )}{2 D _{∥}}, b = \frac{u _{0} ( - 1 + ξ )}{2 D _{∥}} and ξ = 1 + \frac{4 k _{⊥}^{2} D _{⊥} D _{∥}}{u _{0}^{2}} .

where a = \frac{u _{0} ( - 1 - ξ )}{2 D _{∥}}, b = \frac{u _{0} ( - 1 + ξ )}{2 D _{∥}} and ξ = 1 + \frac{4 k _{⊥}^{2} D _{⊥} D _{∥}}{u _{0}^{2}} .

\frac{\partial N}{\partial t} = \frac{2 l}{π} \int_{0}^{L} d z \int_{p}^{\infty} [\frac{1}{3} \frac{\partial u _{1}}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f _{1} )}{\partial p}] 4 π p^{2} d p = - \frac{8 l}{3} \int_{0}^{L} \frac{\partial u _{1}}{\partial z} p^{3} f_{1} d z .

\frac{\partial N}{\partial t} = \frac{2 l}{π} \int_{0}^{L} d z \int_{p}^{\infty} [\frac{1}{3} \frac{\partial u _{1}}{\partial z} \frac{1}{p ^{2}} \frac{\partial ( p ^{3} f _{1} )}{\partial p}] 4 π p^{2} d p = - \frac{8 l}{3} \int_{0}^{L} \frac{\partial u _{1}}{\partial z} p^{3} f_{1} d z .

f_{k} (k_{∥}) = \frac{i k _{∥} u _{k} ( k _{∥} )}{3} \frac{( D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2} ) - i k _{∥} u _{0}}{( D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2} ) ^{2} + k _{∥}^{2} u _{0}^{2}} p \frac{\partial f _{0}}{\partial p}

f_{k} (k_{∥}) = \frac{i k _{∥} u _{k} ( k _{∥} )}{3} \frac{( D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2} ) - i k _{∥} u _{0}}{( D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2} ) ^{2} + k _{∥}^{2} u _{0}^{2}} p \frac{\partial f _{0}}{\partial p}

and \frac{\partial N}{\partial t} = - \frac{8 L l}{9 u _{0}^{2}} p^{4} \frac{\partial f _{0}}{\partial p} \int_{- \infty}^{\infty} Γ (k_{∥}) ⟨ u^{2} ⟩_{k} d k_{∥} where Γ (k_{∥}) = k_{∥} u_{0} (\frac{D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2}}{k _{∥} u _{0}} + \frac{k _{∥} u _{0}}{D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2}})^{- 1} .

and \frac{\partial N}{\partial t} = - \frac{8 L l}{9 u _{0}^{2}} p^{4} \frac{\partial f _{0}}{\partial p} \int_{- \infty}^{\infty} Γ (k_{∥}) ⟨ u^{2} ⟩_{k} d k_{∥} where Γ (k_{∥}) = k_{∥} u_{0} (\frac{D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2}}{k _{∥} u _{0}} + \frac{k _{∥} u _{0}}{D _{∥} k _{∥}^{2} + D _{⊥} k _{⊥}^{2}})^{- 1} .

Γ (k_{∥}) = (\frac{u _{0}^{2}}{D _{∥}}) \frac{k _{∥}}{k _{c ∥}} (\frac{k _{∥}^{2} + k _{c ⊥}^{2}}{k _{c ∥} k _{∥}} + \frac{k _{c ∥} k _{∥}}{k _{∥}^{2} + k _{c ⊥}^{2}})^{- 1} .

Γ (k_{∥}) = (\frac{u _{0}^{2}}{D _{∥}}) \frac{k _{∥}}{k _{c ∥}} (\frac{k _{∥}^{2} + k _{c ⊥}^{2}}{k _{c ∥} k _{∥}} + \frac{k _{c ∥} k _{∥}}{k _{∥}^{2} + k _{c ⊥}^{2}})^{- 1} .

k_{c u t} = \frac{u _{0}}{D _{∥}} [\frac{1}{2} + \frac{1}{4} + \frac{k _{⊥}^{4} D _{B o hm}^{2}}{u _{0}^{2}}]^{1/2}

k_{c u t} = \frac{u _{0}}{D _{∥}} [\frac{1}{2} + \frac{1}{4} + \frac{k _{⊥}^{4} D _{B o hm}^{2}}{u _{0}^{2}}]^{1/2}

\frac{\partial N}{\partial t} = - \frac{8 L l}{9 D _{∥}} p^{4} \frac{\partial f _{0}}{\partial p} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥} .

\frac{\partial N}{\partial t} = - \frac{8 L l}{9 D _{∥}} p^{4} \frac{\partial f _{0}}{\partial p} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥} .

τ_{F 2} = \frac{9 π}{2} D_{∥} (- \frac{p}{f _{0}} \frac{\partial f _{0}}{\partial p})^{- 1} (\int_{k_{c u t}}^{\infty} 2 π ⟨ u^{2} ⟩_{k} d k_{∥})^{- 1} .

τ_{F 2} = \frac{9 π}{2} D_{∥} (- \frac{p}{f _{0}} \frac{\partial f _{0}}{\partial p})^{- 1} (\int_{k_{c u t}}^{\infty} 2 π ⟨ u^{2} ⟩_{k} d k_{∥})^{- 1} .

τ_{F 1} = \frac{8 D _{∥}}{u _{s}^{2}}

τ_{F 1} = \frac{8 D _{∥}}{u _{s}^{2}}

\frac{τ _{F 1}}{τ _{F 2}} = \frac{16}{9 π} (- \frac{p}{f _{0}} \frac{\partial f _{0}}{\partial p}) (\frac{1}{u _{s}^{2}} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥}) .

\frac{τ _{F 1}}{τ _{F 2}} = \frac{16}{9 π} (- \frac{p}{f _{0}} \frac{\partial f _{0}}{\partial p}) (\frac{1}{u _{s}^{2}} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥}) .

\frac{p ^{4}}{D _{∥}} \frac{\partial f _{0}}{\partial p} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥} = constant .

\frac{p ^{4}}{D _{∥}} \frac{\partial f _{0}}{\partial p} \int_{k_{c u t}}^{\infty} ⟨ u^{2} ⟩_{k} d k_{∥} = constant .

k_{m \pm} u_{\pm} - k_{m \pm}^{2} D_{∥ \pm} + k_{m}^{2} D_{⊥ \pm} = 0

k_{m \pm} u_{\pm} - k_{m \pm}^{2} D_{∥ \pm} + k_{m}^{2} D_{⊥ \pm} = 0

k_{m \pm} D_{∥ \pm} = \frac{1}{2} (u_{\pm} \mp u_{\pm}^{2} + 4 k_{m}^{2} D_{∥ \pm} D_{⊥ \pm}) .

k_{m \pm} D_{∥ \pm} = \frac{1}{2} (u_{\pm} \mp u_{\pm}^{2} + 4 k_{m}^{2} D_{∥ \pm} D_{⊥ \pm}) .

\frac{1}{3 p ^{2} f _{m s}} \frac{\partial ( p ^{3} f _{m s} )}{\partial p} = 1 - \frac{k _{m +} D _{∥ +} - k _{m -} D _{∥ -}}{u _{+} - u _{-}}

\frac{1}{3 p ^{2} f _{m s}} \frac{\partial ( p ^{3} f _{m s} )}{\partial p} = 1 - \frac{k _{m +} D _{∥ +} - k _{m -} D _{∥ -}}{u _{+} - u _{-}}

\frac{p}{3 f _{s m}} \frac{\partial f _{s m}}{\partial p} = - \frac{1}{2} - \frac{u _{-}^{2} + 4 k _{m}^{2} D _{∥ -} D _{⊥ -} + u _{+}^{2} + 4 k _{m}^{2} D _{∥ +} D _{⊥ +}}{2 ( u _{-} - u _{+} )} .

\frac{p}{3 f _{s m}} \frac{\partial f _{s m}}{\partial p} = - \frac{1}{2} - \frac{u _{-}^{2} + 4 k _{m}^{2} D _{∥ -} D _{⊥ -} + u _{+}^{2} + 4 k _{m}^{2} D _{∥ +} D _{⊥ +}}{2 ( u _{-} - u _{+} )} .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Cosmic Ray Acceleration in Hydromagnetic Flux Tubes

AR Bell1, JH Matthews2, KM Blundell2, and AT Araudo3,4

1University of Oxford, Clarendon Laboratory, Parks Road, Oxford OX1 3PU, UK

2University of Oxford, Astrophysics, Keble Road, Oxford OX1 3RH, UK

3Astronomical Institute of the Czech Academy of Sciences, Bocni II 1401, CZ-14100 Prague, Czech Republic

4ELI Beamlines, Institute of Physics, Czech Academy of Sciences, 25241 Dolní Břežany, Czech Republic E-mail:[email protected]

Abstract

We find that hydromagnetic flux tubes in back-flows in the lobes of radio galaxies offer a suitable environment for the acceleration of cosmic rays (CR) to ultra-high energies. We show that CR can reach the Hillas (1984) energy even if the magnetised turbulence in the flux tube is not sufficiently strong for Bohm diffusion to apply. First-order Fermi acceleration by successive weak shocks in a hydromagnetic flux tube is shown to be equivalent to second-order Fermi acceleration by strong turbulence.

keywords:

cosmic rays, acceleration of particles, shock waves, magnetic field, radio galaxies

††pagerange: Cosmic Ray Acceleration in Hydromagnetic Flux Tubes–LABEL:lastpage††pubyear: 2019

1 INTRODUCTION

The origin of ultra-high energy cosmic rays (UHECR) is uncertain, and many possible sources have been proposed. Radio galaxies have long been considered a likely source of UHECR because of their high power, large size and longevity (eg Hardcastle et al 2009, O’Sullivan et al 2009, Wykes et al 2013, Eichmann et al 2018). In Matthews et al (2018, 2019) we made the case from observations and theory that UHECR may be accelerated by shocks in the lobes of radio galaxies. Our numerical simulations demonstrated the presence of shocks in concentrated, approximately annular, flows centred on the jet axis, which we refer to as hydromagnetic flux tubes, or more simply as ‘flux tubes’, that emerge from the high pressure hot-spots at the end of relativistic jets. Back-flows in radio lobes have been discussed by many authors (for example Norman et al 1982, Falle 1991, Scheuer 1995, Saxton et al 2002, Krause 2005, Keppens et al 2008, Mignone et al 2010, Cielo et al 2014, Tchekhovskoy & Bromberg 2016). The flow expands into the lobes where the pressure is lower. As discussed in Matthews et al (2019), from application of the Bernoulli equation, the flow becomes supersonic before being slowed down in one or more shocks as it progresses deeper into the lobe as depicted in an idealised model in figure 1. Our simulations of radio jets support this picture by demonstrating the presence of a succession of shocks in flux tubes with Mach numbers of a few and flow velocities of the order of $c/4$ . With these velocities, shocks are well suited to UHECR acceleration since fully relativistic shocks pose severe problems for acceleration to the highest energies (Kirk & Reville 2010; Lemoine & Pelletier 2010; Reville & Bell 2014; Bell et al 2018). Further discussion of our model and the observational context can be found in Matthews et al (2018, 2019). Simulations show that flux tubes are usually more contorted than the idealised flux tubes in figure 1. However, the flux tube need only be long enough to contain the CR diffusion scaleheights upstream and downstream of a shock, and curvature should not invalidate the model since magnetic field lines threading the tube should cause CR to diffuse predominantly parallel to the tube.

Cosmic ray (CR) acceleration by diffusive shock acceleration (Krymskii 1977, Axford et al 1977, Bell 1978a,b, Blandford & Ostriker 1978) has been studied at length, both theoretically and observationally, in the context of supernova remnants (SNR). The maximum energy of CR accelerated by SNR is constrained by various factors. The most general constraint is set by the physical size $R$ of the shock, leading to a maximum CR energy of $ZuBR$ (Hillas 1984), where $u$ is the flow velocity, $B$ is the magnetic field, and the CR has charge $Ze$ . $ZuBR$ is the characteristic maximum CR energy in eV if $u$ , $B$ & $R$ are in SI units. A related constraint is set by the need for the shock to persist for a time longer than the time taken for CR to be accelerated (Lagage & Cesarsky 1983a,b). For a quasi-parallel shock (magnetic field parallel to or at an angle significantly smaller than $\pi/2$ to the shock normal) the Lagage & Cesarsky limit is equivalent to the Hillas limit if the shock persists for a flow time $R/u$ . For perpendicular shocks, the acceleration is very rapid and the maximum CR energy is set by the Hillas limit (Jokipii 1982,1987). In this paper we assume that the hydrodynamic flow through a flux tube is in steady state on the scale of the time taken for CR acceleration, in which case the limit on the maximum CR energy is better understood through the spatial limit of Hillas than the timescale analysis of Lagage & Cesarsky.

The characteristics of backflows in relation to Bernoulli’s equation are discussed by Matthews et al (2019) where it is found that relatively narrow flux tubes are surrounded by vortices as idealised in figure 1. Flux tubes have two characteristic scale-lengths, namely the length and width of the flux tube, and it is not immediately clear which of these, if either, should be invoked when applying the Hillas limit.

In the case of SNR, CR acceleration can be limited by the time taken to amplify the magnetic field (Zirakashvili & Ptuskin 2008, Bell et al 2013), but this limitation is less stringent in flux tubes since (i) the backflow passes through a succession of shocks such that amplification at one shock prepares a magnetic field for subsequent shocks, (ii) a large magnetic field is generated as the plasma passes from the jet through the termination hot-spot into the backflow (Araudo et al 2015, 2018). Hence the situation is different from that of a single isolated shock at the outer boundary of a SNR. A different analysis is needed to determine the maximum CR energy attained in a steady long-lived back-flow. Here we construct an idealised model designed to reproduce the dominant features relevant to CR acceleration by shocks in flux tubes in radio galaxies.

In section 2 we consider first-order acceleration at a single stationary shock in a non-relativistic back-flow. Acceleration at shocks with finite spatial extent has been considered previously in the context of the Earth’s bow shock (Eichler 1981, Lee 1982, Webb et al 1985). Eichler’s analysis exhibits many of the features we find to apply to flux tubes, most notably that the maximum CR energy is independent of the strength of the turbulent scattering of CR by fluctuations in the magnetic field.

In section 3, we apply second-order Fermi theory (Fermi 1949, Jones 1994) to acceleration by strong turbulence in a flux tube and demonstrate the equivalence between this and first-order acceleration by successive weak shocks.

Our analysis focusses on the application of our model to the lobes of radio galaxies, but there may be other scenarios in which the model might apply such as particle acceleration at shocks inside jets with a variable speed, shocks in jets encountering stars or clouds, or reconfinement shocks in jets. Depending on the circumstances, the model may need to be extended to oblique shocks or relativistic shocks.

We use SI units throughout with CR energies in eV. We use the symbol $T$ for CR energies for consistency with related papers where we wish to avoid possible confusion with the electric field. UHECR may consist of nuclei with various charges and masses. To avoid unnecessary mathematical complication, we consider only the acceleration of protons. The behaviour of heavier nuclei can be read from our equations by noting that if a proton with charge $e$ is accelerated to ultra-high energy $T$ by the processes considered here, then a nucleus with charge $Ze$ is automatically accelerated to energy $ZT$ by exactly the same processes.

2 ACCELERATION AT A SINGLE STATIONARY SHOCK IN A FLUX TUBE

Our model for UHECR acceleration at a single shock in a back-flow is visualised in figure 2. We assume 2D geometry in which the flux tube, which might be more accurately referred to as a flux slab, is infinite and uniform in the $y$ direction. Slab geometry is closer to the actual geometry in a radio lobe where the flow geometry resembles in cross-section an annulus extending in angle around the jet, as may be imagined from figure 1 and as seen in Matthews et al (2019). The flux tube contains a stationary shock at $z=0$ where $z$ is the distance along the flux tube. The flow is from left to right in $z$ . The flux tube of width $l$ occupies the space $-l/2<x<l/2$ across the tube. The flux tube is assumed to be straight in $z$ on the basis that it bends over distances larger than diffusion lengths characterising the process of particle acceleration. The plasma flows at velocity $u(z)$ with a velocity $u_{-}$ into the shock, and exits the shock with a downstream velocity $u_{+}$ .

We assume that there is a uniform magnetic field aligned along the flow direction, $B=B_{\parallel}$ both upstream and downstream of the shock. This parallel magnetic field is made plausible by the assumption that all the plasma in the tube has recently passed through the termination shock, so it is magnetically disconnected from plasma in the lobe outside the tube. The crucial consequence in our model is that CR propagating along a magnetic field line stay inside the flux tube and can only escape from the sides of the flux tube by diffusing across the magnetic field. Moreover, the plasma in the flux tube is stretched in the $z$ direction as it expands out of the hot spot to become supersonic as it passes through the lobe, thus making the field predominantly locally

parallel to the flow rather than perpendicular across the width of the tube. As discussed below, CR diffusion has different diffusion coefficients $D_{\parallel}$ and $D_{\perp}$ parallel and perpendicular to the magnetic field. Overall, the model depicted in figure 2 is a simplification but it elucidates basic principles that could be expanded upon in further work.

The equation for the evolution of the CR distribution function $f$ is (eg Blandford & Eichler 1987)

[TABLE]

where $p$ is the magnitude of momentum, and $f$ is defined in the local rest frame of a plasma moving at velocity $u(z)$ . In this section we assume that $u$ , $D_{\perp}$ & $D_{\parallel}$ are spatially uniform except that they jump discontinuously in value across the shock. Separation of variables gives a steady state solution ( $\partial f/\partial t=0$ ) of the form

[TABLE]

where $f_{-}$ and $f_{+}$ are the solutions upstream and downstream respectively. $f_{sm}(p)$ is the value of the $m$ th Fourier component of both $f_{-}$ and $f_{+}$ at the shock ( $z=0$ ) since $f$ is continuous across the shock: $f_{s}(x,p)=f_{-}(0,x,p)=f_{+}(0,x,p)$ . The boundary condition at the edges of the flux tube is that $f=0$ on the assumption that CR reaching the boundary escape freely. $\partial u/\partial z=0$ everywhere except at $z=0$ .

As shown in the Appendix (equation A6), the time-independent solution of equation (1) is

[TABLE]

where $F_{m}$ are constants of integration and $k_{m}=(2m+1)\pi/l$ . From equation (1), the rate at which the shock accelerates CR to momenta greater than $p$ is

[TABLE]

The resulting differential energy spectrum is $n(p)=-dN/dp$ . Differentiating equation (4), and using equation (A4) for $\partial(p^{3}f_{sm})/\partial p$ , the rate at which CR are produced with momentum $p$ is

[TABLE]

We can relate the parallel and perpendicular diffusion coefficients to the Bohm diffusion coefficient defined as

[TABLE]

where $B_{\parallel}$ is the magnetic field along the flux tube on a scale much larger than the CR Larmor radius $r_{g}$ . In standard plasma theories for the diffusion of particles with a scattering time $\tau_{scat}$ longer than the gyration time ( $\omega_{g}\tau_{scat}>1$ where $\omega_{g}=ceB_{\parallel}/p$ ), the diffusion coefficient across a magnetic field is $D_{\perp}=D_{Bohm}/(\omega_{g}\tau_{scat})$ , and the diffusion coefficient parallel to the magnetic field is $D_{\parallel}=D_{Bohm}\times(\omega_{g}\tau_{scat})$ to within numerical factors of order one that depend on the nature of the scattering (eg Braginskii 1965, Blandford & Eichler 1987). For a parallel shock, the zeroth order uniform parallel magnetic field, and therefore the Bohm diffusion coefficient, is the same upstream and downstream. Consequently, $\sqrt{D_{\parallel-}D_{\perp-}}=\sqrt{D_{\parallel+}D_{\perp+}}=D_{s}$ , where $D_{s}$ is the Bohm diffusion coefficient near the shock.

As shown in the above equations (3-5), CR acceleration depends on $D_{s}$ alone, and not on $D_{\perp}$ and $D_{\parallel}$ separately. This was previously demonstrated by Eichler (1981) who highlighted the important result that the multiple $D_{s}^{2}=D_{\parallel}D_{\perp}$ does not depend on the strength of the turbulence since the factor $\omega_{g}\tau_{scat}$ cancels out. This contrasts favourably with acceleration by an isolated shock where the maximum energy is reduced from the Hillas energy if the CR mean free path is greater than its Larmor radius and $D_{\parallel}>D_{Bohm}$ , in which case $T_{max}\sim(D_{\parallel}/D_{Bohm})^{-1}u_{s}Bl$ . In young SNR this is not a major concern since observation and theory both indicate that the CR diffusion coefficient is close to $D_{Bohm}$ (Stage et al 2006, Uchiyama et al 2007, Bell 2014), but more generally this need not be so. Hence a significant advantage of acceleration in a flux tube is that the maximum energy robustly depends on $D_{Bohm}$ and the maximum CR proton energy is $T_{max}\sim u_{s}Bl$ .

If CR with low momenta are injected into the acceleration process uniformly across the width of the shock, $F_{m}=(-1)^{m}(4/\pi)F_{s}/(2m+1)$ where the spectrum at low momentum is $F_{s}p^{-4}$ . For a strong shock with velocity $u_{s}$ , $u_{-}=u_{s}$ and $u_{+}=u_{s}/4$ . giving

[TABLE]

In equations (7), $f_{sm}$ relates to the CR momentum spectrum at the shock while $\partial n/\partial t$ relates to the momentum spectrum of CR released from the flux tube into the surrounding lobe. $\partial n/\partial t$ is the crucial quantity if we are interested in the lobes of radio galaxies as a source of UHECR. Note that the spectrum of CR escaping the flux tube is not the same as the spectrum at the shock. $\partial n/\partial t$ is plotted in figure 3. The $m=0$ mode makes the major contribution. The neglect of higher $m$ modes would not make a discernable difference to the curve in figure 3.

At low CR momentum, the spectrum follows the usual $p^{-2}$ power law for acceleration at a strong shock. Figure 3 shows that the momentum spectrum turns over where $D_{s}/u_{s}l\approx 0.2$ , which corresponds to a maximum proton energy in eV of

[TABLE]

Nuclei with charge $Ze$ can be expected to be accelerated to $Z$ times this energy. The maximum energy in equation (8) is close to the Hillas energy defined as $u_{s}Bl$ . $T_{max}$ is independent of the value of $\omega_{g}\tau_{scat}$ . Our analysis assumes that the entrance to the flux tube is further from the shock than the diffusion length $D_{\parallel-}/u_{s}$ upstream of the shock, and that the exit end of the flux tube is further from the shock than the diffusion length $4D_{\parallel+}/u_{s}$ downstream of the shock. If either of these conditions is violated the maximum CR energy is determined by losses from the relevant end of the flux tube.

The reason the maximum CR energy is independent of $\omega_{g}\tau_{scat}$ is that the acceleration rate is proportional to $D_{\parallel}^{-1}\propto(\omega_{g}\tau_{scat})^{-1}$ , and the loss rate from the sides of the flux tube is proportional to $D_{\perp}$ which is also proportional to $(\omega_{g}\tau_{scat})^{-1}$ . The maximum CR energy is the energy at which the loss rate is equal to the acceleration rate. Since both the acceleration rate and the loss rate are proportional to $(\omega_{g}\tau_{scat})^{-1}$ , $T_{max}$ is independent of $\omega_{g}\tau_{scat}$ .

3 Comparison with 2nd-order Fermi acceleration

Simulations by Matthews et al (2019) followed fluid elements (modelled as tracer particles embedded in the flow) as they pass along the jet into the hotspot at the head of the jet, and then emerge through a back-flowing flux tube into the radio lobe. Analysis of the history of tracer particles shows that the flow through the flux tube is generally transonic with supersonic episodes during which shocks occur. About 10 percent of the tracer particles pass through a shock with a Mach number greater than three where strong CR acceleration can be expected. More generally, tracer particles pass through more frequently occuring weak shocks with low Mach numbers. We could repeat the analysis of the previous section for a multiplicity of weak shocks distributed randomly along the length of the flux tube. However, a collection of many weak shocks can equivalently be described as turbulence with a turbulent velocity comparable with the large-scale background flow velocity. Flow though weak shocks can be treated as a small velocity perturbation that can be decomposed into Fourier modes as discussed below. The flow structure in the following analysis is one-dimensional, and therefore does not qualify correctly as turbulence, but the underlying acceleration process is closely similar to second-order Fermi acceleration in a turbulent plasma and for ease of expression we will refer to the disordered velocity field as turbulence.

In this section we analyse CR acceleration by considering second-order Fermi acceleration in a flux tube with a time-independent spatially-varying flow velocity directed along the flux tube. The configuration is similar to that in figure 2 except that changes in flow velocity $u$ are no longer localised at a single shock but distributed along the length of the tube. The flow velocity is assumed to be uniform across the width of the flux tube:

[TABLE]

This formulation can be used to represent flow along a flux tube as it passes through a series of constrictions imposed by the surrounding lobe. For transonic flow containing weak shocks $u_{1}\sim u_{0}$ , but we make the approximation appropriate for second-order Fermi that $u_{1}\ll u_{0}$ and keep terms of order $(u_{1}/u_{0})^{2}$ . The implications of increasing the ratio $u_{1}/u_{0}$ to $\sim 1$ can then be considered. In the rest frame of plasma flowing along the flux tube, this looks like waves with amplitude $u_{1}$ travelling in reverse direction toward the hotspot at velocity $-u_{0}$ . The transfer of energy from the turbulence to the CR can be interpreted as transit-time damping in which CR extract energy from the wave by diffusing out of dense compressed regions on the same timescale as the flow alternates between compression and rarefaction. This damping of hydrodynamic waves is well known in other contexts (eg Bell 1983). It is a manifestation of second-order Fermi acceleration since CR preferentially spend time in parts of the flow undergoing compression where they scatter between converging fluid elements.

We apply the free boundary escape condition ( $f=0$ ) at the edges of the flux tube. The CR distribution can be Fourier expanded, $f_{m}\propto\cos((2m+1)\pi x/l)$ as in equation (2). In the case of acceleration by a single shock, the $m=0$ mode is the main source of CR acceleration. The same applies here, so we develop the analysis by ignoring all modes except $m=0$ . As part of the linearisation procedure, we express the CR distribution function as a sum of zeroth and first order contributions:

[TABLE]

where $k_{\perp}=\pi/l$ . $f_{0}$ and $u_{0}$ are both uniform along the flux tube. The analysis proceeds by linearising equation (1) for the CR distribution:

[TABLE]

where $D_{\perp}$ & $D_{\parallel}$ are uniform both along (in $z$ ) and across (in $x$ ) the flux tube. The solution to equation (11) is

[TABLE]

Equation (12) says that any local gradient in $u_{1}$ , which might be viewed as a weak shock, causes a disturbance in the surrounding CR population that decays exponentially over a distance $|b|^{-1}$ upstream and a distance $|a|^{-1}$ downstream, where $a<0$ and $b>0$ .

As in equation (4) of section (2), the rate at which the number of CR with momenta greater than p increases in a long flux tube of length $L$ , $L\gg\max(|a|^{-1},|b|^{-1})$ , is

[TABLE]

We now express $u_{1}(z)$ and $f_{1}(z)$ as Fourier transforms: $u_{1}(z)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}u_{k}(k_{\parallel})e^{ik_{\parallel}z}dk_{\parallel}$ and $f_{1}(z)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}f_{k}(k_{\parallel})e^{ik_{\parallel}z}dk_{\parallel}$ . The application of standard Fourier procedures to equation (11) leads to

[TABLE]

$\langle u^{2}\rangle_{k}dk_{\parallel}$ is the mean square velocity of turbulence in the range $k_{\parallel}$ to $k_{\parallel}+dk_{\parallel}$ . $\Gamma(k_{\parallel})$ has the dimensions of a rate in time. $\Gamma(k_{\parallel})$ controls the contribution of turbulence at wavenumber $k_{\parallel}$ to the acceleration of CR with momentum $p$ . At very long wavelengths $\Gamma(k_{\parallel})$ is proportional to $k_{\parallel}^{2}$ : $\Gamma(k_{\parallel})\rightarrow(u_{0}^{2}/D_{\perp})(k_{\parallel}^{2}/k_{\perp}^{2})$ . At very short wavelengths $\Gamma(k_{\parallel})$ is constant: $\Gamma(k_{\parallel})\rightarrow u_{0}^{2}/D_{\parallel}$ . At intermediate wavelengths, there are two critical wavenumbers: (a) $k_{c\parallel}=u_{0}/D_{\parallel}$ , which determines whether the inverse wavenumber $k_{\parallel}^{-1}$ is shorter ( $k_{\parallel}>k_{c\parallel}$ ) than the characteristic diffusion length along the flux tube, and (b) $k_{c\perp}=k_{\perp}\sqrt{D_{\perp}/D_{\parallel}}$ which determines whether parallel diffusion dominates (for $k_{\parallel}>k_{c\perp}$ ) over perpendicular diffusion. $\Gamma(k_{\parallel})$ can be written in terms of these two critical wavenumbers:

[TABLE]

For comparison with shock acceleration, apart from a dimensionless numerical factor, ${u_{0}^{2}}/{D_{\parallel}}$ is the acceleration rate of CR with diffusion coefficient $D_{\parallel}$ at a shock with velocity $u_{0}$ . The function $\Gamma(k_{\parallel})$ is plotted in figure 4.

If the cut-off wavenumber $k_{cut}$ is defined as the value of $k_{\parallel}$ at which $\Gamma=0.5u_{0}^{2}/D_{\parallel}$ , then

[TABLE]

where $D_{\perp}D_{\parallel}=D_{Bohm}^{2}$ and $D_{Bohm}$ is the Bohm diffusion coefficient as defined in equation (6). Once again we see that diffusion through the sides of the flux tube is controlled by the multiple of $D_{\perp}$ and $D_{\parallel}$ which is independent of the strength of the turbulence as represented by $\omega_{g}\tau_{scat}$ .

From equation (16), power-law turbulence at long wavelengths, $k_{\parallel}<k_{cut}$ , $\langle u^{2}\rangle_{k}\propto k_{\parallel}^{-\gamma}$ , does not contribute significantly to CR acceleration provided $\gamma<2$ . CR either escape through the side of the flux tube before gaining energy, or else remain adiabatically fixed in the same fluid element without being accelerated. On the other hand, at short wavelengths, $k_{\parallel}>k_{cut}$ , power law turbulence makes little contribution to CR acceleration provided $\gamma>0$ since $\langle u^{2}\rangle_{k}\propto k_{\parallel}^{-\gamma}$ decreases towards large $k_{\parallel}$ . Consequently, power law turbulence with $0<\gamma<2$ contributes most strongly to CR acceleration at wavenumbers close to $k_{cut}$ . To a reasonable approximation, since $\Gamma(k_{\parallel})\approx u_{0}^{2}/D_{\parallel}$ for $k_{\parallel}>k_{cut}$ ,

[TABLE]

The rate $dp/dt$ at which CR gain momentum is related to $\partial N/\partial t$ by the continuity equation: ${\partial N}/{\partial t}=Ll\ 4\pi p^{2}f_{0}(dp/dt)$ . The acceleration time for second-order Fermi acceleration time is $\tau_{F2}=p/(dp/dt)$ . Hence from equation (19)

[TABLE]

For comparison, the acceleration time for first-order Fermi shock acceleration (Lagage & Cesarsky 1983a,b) is

[TABLE]

where the representative numerical factor of 8 depends on the ratio of the diffusion coefficients upstream and downstream of the shock. Hence, the ratio of the first- and second-order Fermi acceleration times is

[TABLE]

The precise value of the factor of $16/9\pi$ has limited significance since it is sensitive to assumptions and approximations built into our derivation. More crucially, equation (22) shows, as might be expected, that the rate of second-order Fermi acceleration is of the same order as that of shock acceleration if the characteristic velocity of the turbulence contributing to CR acceleration is comparable to the shock velocity.

The acceleration rate is important since the maximum energy reached by CR is determined by competition between acceleration and loss. $\tau_{F2}\sim\tau_{F1}$ is required if the second-order Fermi process is to accelerate CR to the same energy that can be achieved by shock acceleration (equation 8). That is, the turbulent flow velocity needs to be comparable with the velocity at which the plasma flows down the flux tube if second-order is to compete with first-order Fermi acceleration. For second-order Fermi acceleration to contribute significantly to the acceleration of UHECR, the turbulence must have an energy density comparable to the kinetic energy density of large scale flows in the lobe as previously indicated by Hardcastle (2010) and consistent with O’Sullivan et al (2009) who find that wave speeds in excess of $0.1c$ are needed for UHECR acceleration by stochastic Alfven waves in the radio galaxy Centaurus A.

The other main point for comparison between first- and second-order Fermi acceleration is the slope of the CR spectrum. CR escape through the walls of the flux tube is responsible for the turnover in the spectrum at high energy, but escape through the walls has negligible effect on the spectrum at lower energies. At energies below the turnover, the spectrum at a single isolated shock is determined by competition between acceleration at the shock and CR escape due to advection downstream away from the shock. Both the acceleration rate and the advective loss rate are well determined in shock acceleration, leading naturally to $n(p)\propto p^{-2}$ below the spectral turnover as in equation (7). However, the situation is more complicated in the case of multiple weak shocks. The first order Fermi spectrum is not simply $n(p)\propto p^{-2}$ since the spectrum is steepened by low compression at weak shocks, yet conversely it is flattened by passage through multiple shocks (Bell 1978b).

In the case of second-order Fermi acceleration both the acceleration rate and the loss rate can take a range of values, thus leading to a range of possible spectra depending on competition between the two rates. The acceleration rate is proportional to the amplitude of the turbulence (equation 20). The only possible escape route for low energy CR is by advection with the fluid flow through the downstream exit from the flux tube. If the plasma transit time $L/u_{0}$ through the flux tube is much longer than the acceleration time, CR loss by advection during acceleration is also negligible, and CR number conservation dictates that $\partial N/\partial t$ is the same at all momenta. From equation (19), this leads to a CR momentum spectrum determined by

[TABLE]

According to equation (23) a range of CR spectra are possible in this steady state, depending on the wavenumber spectrum of the turbulence. Turbulence in a flux tube is a special case, but if the turbulence follows a Kolmogorov power law with energy density proportional to $k^{-5/3}$ , such that $\int_{k_{cut}}^{\infty}\langle u^{2}\rangle_{k}dk_{\parallel}\propto p^{2/3}$ , since $k_{cut}\propto D_{\parallel}^{-1}\propto p^{-1}$ , then $f_{0}(p)\propto p^{-8/3}$ and $n(p)\propto p^{-2/3}$ . This would be an unusually flat CR spectrum with most of the CR energy density at the large momentum end of the spectrum. The cause of the flat spectrum is that CR acceleration is slower at high CR momentum, with the result that a relatively greater number of CR are needed at high momentum to maintain a constant flux of CR from low to high momentum. The momentum spectrum given by equation (23) is unlikely except in special circumstances, and an extended calculation is needed that takes account of time dependence, escape through the ends of the flux tube, and feedback processes across the momentum spectrum.

4 Conclusions

We examine the acceleration of CR in hydromagnetic flux tubes with application to the lobes of radio galaxies where flux tubes contain shocks that are mostly low Mach number (Matthews et al 2019). We consider CR acceleration in the framework of two overlapping formalisms: first-order Fermi acceleration by shocks and second-order Fermi acceleration by strong turbulence. We show that the maximum CR energy at a single shock in a flux tube is close to the Hillas energy $ZuBl$ where $u$ is the flow velocity along the flux tube, $B$ is an ordered field aligned along the flux tube, $l$ is the width of the flux tube, and the CR charge is $Ze$ . In contrast to CR acceleration at the outer shocks of SNR, the attainment of the Hillas energy is not dependent on Bohm diffusion in the shock environment. Consequently, the Hillas energy is a relatively robust estimate of the energy to which CR are accelerated.

We show that second-order Fermi acceleration in turbulence is slower than first-order shock acceleration by the ratio of the mean square turbulent velocity to the square of the shock velocity. Hence, second-order acceleration in weak turbulence is less able to compete with CR losses, and cannot accelerate CR to the highest energies. However, if the turbulence is strong in the sense that the mean square turbulent velocity is comparable with the large scale flow velocity, then second-order acceleration competes well with first-order acceleration since it is broadly equivalent to multiple first-order acceleration episodes as a CR passes through a succession of weak shocks.

The flow within lobes of radio galaxies is more complicated (see references in section 1) than allowed by our simplifed model of one-dimensional flow in a steady uniform straight flux tube, and additional acceleration processes may be active. For example, the strong velocity differential across the edges of flux tubes may be a prime location for shear acceleration (Earl et al 1988, Ostrowski 1998, Rieger & Duffy 2006, Caprioli 2015) and for the growth of strong multi-dimensional turbulence due to the Kelvin-Helmholtz instability. Also, variations in the width of flux tubes along their length may give rise to lateral compressions and rarefactions that respectively accelerate and decelerate CR, although resulting changes in CR energy may be approximately adiabatic with limited overall effect. As shown by Achterberg (1981), second-order Fermi acceleration is an umbrella term for a number of different effects that can be combined in a general quasi-linear theory of wave-particle interactions. The more general theory may provide a fruitful basis for an extended analysis in the complicated environment of a flux tube, and it may be possible to encompass first- and second-order Fermi acceleration and shear acceleration in a single comprehensive analysis (Lemoine 2019).

The model we use here could be extended to include the above additional effects, but within its limitations, our analysis supports the proposal by Matthews et al (2018, 2019) that UHECR may be accelerated in the lobes of radio galaxies. Suitably extended, our analysis may be applicable to other steady confined flows such as jets.

5 ACKNOWLEDGEMENTS

We thank Bram Achterberg for helpful comments. This research was supported by the UK Science and Technology Facilities Council under grant No. ST/N000919/1.

6 REFERENCES

Achterberg A., 1981, A&A, 97, 259

Araudo A.T., Bell A.R., Blundell K.M., 2015, ApJ, 806, 304

Araudo A.T., Bell A.R., Blundell K.M., Matthews J.H., 2018, MNRAS, 473, 3500

Axford W.I., Leer E., Skadron G., 1977, Proc 15th Int. Cosmic Ray Conf., 11, 132

Bell A.R., 1978a, MNRAS, 182, 147

Bell A.R., 1978b, MNRAS, 182, 443

Bell A.R., 1983, Phys Fluids, 26, 269

Bell AR, 2014, Braz J Phys, 44, 415

Bell A.R., Araudo A.T., Matthews J.H., Blundell K.M., 2018, MNRAS, 473, 2364

Bell A.R., Schure K.M., Reville B., Giacinti G.,2013, MNRAS, 431, 415

Blandford R.D., Eichler D., 1987, Phys Rep, 154,1

Blandford R.D., Ostriker J.P., 1978, ApJ, 221, L29

Braginskii S.I., 1965, Rev Plasma Phys, 1, 205

Caprioli D., 2015, ApJ, 811, L38

Cielo, S., Antonuccio-Delogu, V., Macciò, A. V., Romeo, A. D., Silk, J., 2014, MNRAS, 439, 2903

Earl J.A., Jokipii J.R., Morfill G., 1988, ApJ, 331, L91

Eichler D., 1981, ApJ, 244, 711

Eichmann, B., Rachen, J. P., Merten, L., van Vliet, A., Becker Tjus, J., 2018, JCAP02, 036

Falle S.A.E.G., 1991, MNRAS, 250, 581

Fermi E., 1949, Phys Rev, 75, 8

Hardcastle M.J., 2010, MNRAS, 405, 2810

Hardcastle M.J., Cheung C.C., Feain I.J., Stawarz L., 2009, MNRAS, 393, 1041

Hillas A.M., 1984, ARA&A, 22, 425

Jones F.C., 1994, ApJSS, 90, 561

Jokipii J.R., 1982, ApJ, 255, 716

Jokipii J.R., 1987, ApJ, 313, 842

Keppens R., Meliani Z., van der Holst B., Casse F., 2008, A&A, 483, 663

Kirk J.G., Reville B., 2010, ApJL, 710, L16

Krause, M., 2005, A&A, 431, 45

Krymskii G.F., 1977, Sov Phys Dokl, 23, 327

Lagage P.O., Cesarsky C.J., 1983a, A&A, 118, 223

Lagage P.O., Cesarsky C.J., 1983b, ApJ, 125, 249

Lee M.A. 1982, J Geophys Res, 87, 5063

Lemoine M., 2019, Phys Rev D, 99, 083006

Lemoine M., Pelettier G., 2010, MNRAS, 402, 321

Matthews J.H., Bell A.R., Blundell K.M., Araudo A.T., 2018, MNRAS, 479, L76

Matthews J.H., Bell A.R., Blundell K.M., Araudo A.T., 2019, MNRAS, 482, 4303

Mignone, A., Rossi, P., Bodo, G., Ferrari, A., Massaglia, S., 2010, MNRAS, 402, 7

Norman M.L. Winkler K.-H., Smarr L., Smith M.D., 1982, A&A, 113, 285

Ostrowski M, 1998, A&A, 335, 134

O’Sullivan S., Reville B., Taylor A.M., 2009, MNRAS, 400, 248

Reiger F.M., Duffy P., 2006, ApJ, 652, 1044

Reville B., Bell A.R., 2014, MNRAS, 439, 2050

Saxton C.J., Sutherland R.S., Bicknell G.V., Blanchet G., Wagner S.J., Metchnik M.V., 2002, A&A, 393, 765

Scheuer P.A.G.S., 1995, MNRAS, 277, 331

Stage M.D., Allen G.E., Houck J.C. Davis J.E., 2006, Nature Physics, 2, 614

Tchekhovskoy A., Bromberg, O., 2016, MNRAS, 461, 46

Uchiyama Y. et al, 2007, Nature, 449, 576

Webb G.M., Bogdan T.J., Lee M. A., Lerche I., 1985, MNRAS, 215, 341

Wykes S., Croston J.H., Hardcastle M.J., Eilek J.A., Biermann P.L., Achterberg A., Bray J.D., Lazarian A., Haverkorn M., Protheroe R.J., Bromberg O., 2013, A&A, 558, A19

Zirakashvili V.N., Ptuskin V.S., 2008, ApJ, 678, 939

APPENDIX A

In this Appendix we find a solution to equation (1) by decomposing the distribution function $f$ into Fourier modes as defined in equation (2). Since $\partial u/\partial z=0$ everywhere except at $z=0$ , the equation is linear in each of the upstream and downstream plasmas separately, with exponential scaleheights $k_{m\pm}$ in $z$ defined by

[TABLE]

where $k_{m}=(2m+1)\pi/l$ , and $u_{\pm}$ , $D_{\parallel\pm}$ & $D_{\perp\pm}$ are the upstream and downstream values of $u$ , $D_{\parallel}$ & $D_{\perp}$ . The solution of equation (A1) is

[TABLE]

$f$ is continuous across the shock $f_{s}(p,x)=f_{-}(0,x,p)=f_{+}(0,x,p)$ . Integrating equation (1) in $z$ across the shock gives

[TABLE]

where $f_{ms}$ is the value of $f_{m}$ at the shock ( $z=0$ ). Substituting equation (A2) into the equation (A3) gives

[TABLE]

In the limit of a strong shock ( $u_{-}=4u_{+}$ ) in a wide flux tube ( $k_{m}=0$ ), the solution is $f_{s}\propto p^{-4}$ , which is the standard power law solution for diffusive shock acceleration.

We assume that the CR scattering length is proportional to the CR Larmor radius, in which case the multiple of $D_{\perp}$ and $D_{\parallel}$ is proportional to $p^{2}$ . Using the standard integral

[TABLE]

equation (A4) can be integrated to give

[TABLE]

where $F_{m}$ is a constant of integration.