Physics of relativistic collisionless shocks: III The suprathermal particles

TL;DR

This paper analyzes the physics of accelerated particles in relativistic collisionless shocks, providing theoretical estimates of their scattering length and distribution, validated by PIC simulations, advancing understanding of shock particle acceleration mechanisms.

Contribution

It offers a new theoretical estimate of particle scattering length in unmagnetized relativistic shocks and compares it with simulation results, enhancing the understanding of particle acceleration.

Findings

The scattering length scales with particle momentum and turbulence parameters.

Good agreement between theoretical estimates and PIC simulation results.

Provides a detailed distribution function for accelerated particles.

Abstract

In this third paper of a series, we discuss the physics of the population of accelerated particles in the precursor of an unmagnetized, relativistic collisionless pair shock. In particular, we provide a theoretical estimate of their scattering length $l_{scatt}(p)$ in the self-generated electromagnetic turbulence, as well as an estimate of their distribution function. We obtain $l_{scatt}(p) \simeq (\gamma_p /\epsilon_B)(p/\gamma_{\infty} mc)^2 (c/\omega_p)$, with p the particle momentum in the rest frame of the shock front, $\epsilon_B$ the strength parameter of the microturbulence, $\gamma_p$ the Lorentz factor of the background plasma relative to the shock front and $\gamma_{\infty}$ its asymptotic value outside the precursor. We compare this scattering length to large-scale PIC simulations and find good agreement for the various dependencies.

Figures

Equations

$$t_{\mathrm{corr}}\equiv \frac{1}{⟨\mathbf{\delta }\mathbf{F}(t)\cdot \mathbf{\delta }\mathbf{F}(t)⟩}{\int }_0^{+\infty }\mathrm{d}\tau ⟨\mathbf{\delta }\mathbf{F}(t+\tau )\cdot \mathbf{\delta }\mathbf{F}(t)⟩,$$

$$⟨\mathbf{\delta }\mathbf{F}(t+\tau )\cdot \mathbf{\delta }\mathbf{F}(t)⟩\propto ⟨\mathbf{\delta }\mathbf{B}(t+\tau )\cdot \mathbf{\delta }\mathbf{B}(t)⟩.$$

$$t_{\mathrm{corr}}=$$

$$\times {\int }_0^{+\infty }\mathrm{d}\tau e^{-2g\left[x(t+\tau )-x(t)\right]}{e}^{+i{k^{\prime }}^{\alpha }{\Delta x_{|\mathrm{w}}}_{\alpha }}{\mathcal{S}}_{|\mathrm{w}}({\mathbf{k}}^{\prime }),$$

$$⟨\mathbf{\delta }{{\mathbf{B}}_{|\mathrm{w}}}_{{\mathbf{k}}_1^{\prime }}\cdot \mathbf{\delta }{{\mathbf{B}}_{|\mathrm{w}}}_{{\mathbf{k}}_2^{\prime }}^{\star }⟩=\delta \left({\mathbf{k}}_1^{\prime }-{\mathbf{k}}_2^{\prime }\right){\mathcal{S}}_{|\mathrm{w}}({\mathbf{k}}_1^{\prime }).$$

$${\Delta x_{|\mathrm{w}}}^t$$

$${\Delta x_{|\mathrm{w}}}^x$$

$${\Delta x_{|\mathrm{w}}}^{\perp }$$

$${\mathcal{R}}_{{\mathbf{k}}^{\prime }}$$

$$=\{\begin{array}{ll}\pi \delta \left({k^{\prime }}^{\alpha }{\Delta x_{|\mathrm{w}}}_{\alpha }/\tau \right)& (g\to 0)\\ & \\ \frac{4g\mu }{4g^2{\mu }^2+{\left({k^{\prime }}^{\alpha }{\Delta x_{|\mathrm{w}}}_{\alpha }/\tau \right)}^2}& (g\ne 0)\end{array}$$

$${\mathcal{S}}_{|\mathrm{w}}({\mathbf{k}}^{\prime })\equiv$$

$$\times \Theta \left[\Delta k_{\perp }^{\prime }-\left(k_{\perp }^{\prime }-\widehat{k_{\perp }^{\prime }}\right)\right]\Theta \left[k_x^{\prime }+\Delta k_x^{\prime }\right]$$

$$\times \Theta \left[\Delta k_x^{\prime }-k_x^{\prime }\right].$$

$$t_{\mathrm{corr}}=\{\begin{array}{ll}0& \mathrm{if}\Delta k_x^{\prime }<\nu \frac{{\widehat{k^{\prime }}}_{\perp }-\Delta k_{\perp }^{\prime }}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})}\\ & \\ \frac{\pi }{2}\frac{1}{\left(1-v_{{\omega }^{\prime }}\right)}\frac{1}{{\gamma }_{\mathrm{w}}\Delta k_x^{\prime }}& \mathrm{if}\Delta k_x^{\prime }>\nu \frac{{\widehat{k^{\prime }}}_{\perp }-\Delta k_{\perp }^{\prime }}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})},\end{array}$$

$$t_{\mathrm{corr}}=\{\begin{array}{ll}{\alpha }_1\frac{g_{|\mathrm{w}}}{{\gamma }_{\mathrm{w}}\Delta k_{\perp }^{\prime }\widehat{k_{\perp }^{\prime }}}& \mathrm{if}\Delta k_x^{\prime }<\nu \frac{{\widehat{k^{\prime }}}_{\perp }-\Delta k_{\perp }^{\prime }}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})}\\ & \\ {\alpha }_2\frac{1}{\left(1-v_{{\omega }^{\prime }}\right)}\frac{1}{{\gamma }_{\mathrm{w}}\Delta k_x^{\prime }}& \mathrm{if}\Delta k_x^{\prime }>\nu \frac{{\widehat{k^{\prime }}}_{\perp }-\Delta k_{\perp }^{\prime }}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})}\end{array}$$

$$\mu \simeq \frac{{\pi }^x}{{\pi }^t}$$

$$\frac{\mathrm{d}\mu }{\mathrm{d}s}\simeq \frac{1}{{{\pi }^t}^2}\left({\pi }^t\frac{\mathrm{d}{\pi }^x}{\mathrm{d}s}-{\pi }^x\frac{\mathrm{d}{\pi }^t}{\mathrm{d}s}\right),$$

$$\frac{\mathrm{d}{\pi }^{\alpha }}{\mathrm{d}s}=q\left({\pi }_{\beta }-q\delta A_{\beta }\right)\frac{\partial A^{\beta }}{\partial x_{\alpha }}.$$

$$\Delta \mu (t)\simeq$$

$${\delta A^{\prime }}^{\alpha }=e^{-gx}\int \frac{{\mathrm{d}}^3{k}^{\prime }}{(2\pi {)}^3}{{e^{\prime }}_{{\mathbf{k}}^{\prime }}}^{\alpha }{e}^{i{k^{\prime }}_{\mu }{x^{\prime }}^{\mu }},$$

$$\Delta \mu (t)\simeq$$

$$\times {p^{\prime }}_{\alpha }{{e^{\prime }}_{{\mathbf{k}}^{\prime }}}^{\alpha }{e}^{-gx}{e}^{i{k}_{\mu }^{\prime }{x}^{\prime \mu }},$$

$$⟨{{e^{\prime }}_{{\mathbf{k}}_1^{\prime }}}^{\alpha }{{{e^{\prime }}_{{\mathbf{k}}_2^{\prime }}}^{\beta }}^{\star }⟩=(2\pi {)}^3\delta \left({\mathbf{k}}_1^{\prime }-{\mathbf{k}}_2^{\prime }\right){{\mathcal{S}}_{{\mathbf{k}}_1^{\prime }}}^{\alpha \beta }.$$

$${{\mathcal{S}}_{{\mathbf{k}}^{\prime }}}^{\alpha \beta }={\mathcal{S}}_{\parallel }{\delta }_{x^{\prime }}^{\alpha }{\delta }_{x^{\prime }}^{\beta }+{\mathcal{S}}_{\perp }\left({\delta }_{y^{\prime }}^{\alpha }{\delta }_{y^{\prime }}^{\beta }+{\delta }_{z^{\prime }}^{\alpha }{\delta }_{z^{\prime }}^{\beta }\right).$$

$$⟨\delta B^{\prime 2}⟩$$

$$\simeq \int \frac{{\mathrm{d}}^2{k}_{\perp }^{\prime }\mathrm{d}{k}_x^{\prime }}{(2\pi {)}^3}{k}_{\perp }^{\prime 2}{\mathcal{S}}_{\parallel }.$$

$$⟨\Delta {\mu }^2⟩$$

$${\nu }_{\mathrm{scatt}}(p)=\frac{e^2⟨\delta B^2⟩}{p^2}\{\begin{array}{ll}{\alpha }_3\frac{g^{\prime }}{{\gamma }_{\mathrm{w}}\Delta k_{\perp }^{\prime }\widehat{k_{\perp }^{\prime }}}& \left[\Delta k_x^{\prime }\lesssim \frac{\widehat{k_{\perp }^{\prime }}}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})}\right]\\ & \\ {\alpha }_4\frac{1}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})\Delta k_x^{\prime }}& \left[\Delta k_x^{\prime }\gtrsim \frac{\widehat{k_{\perp }^{\prime }}}{{\gamma }_{\mathrm{w}}(1-v_{{\omega }^{\prime }})}\right]\end{array}$$

$${\nu }_{\mathrm{scatt}}\simeq \frac{e^2⟨\delta B^2⟩}{p^2}\frac{1}{{\gamma }_{\mathrm{w}}\overline{k^{\prime }}},$$

$$l_{\mathrm{scatt}}(p)\sim {\gamma }_{\mathrm{p}}{ϵ}_B^{-1}{\left(\frac{p}{{\gamma }_{\infty }m}\right)}^2{\omega }_{\mathrm{p}}^{-1},$$

$${\gamma }_{\mathrm{w}}\left({\beta }_{\mathrm{w}}+{\mu }_{|\mathrm{w}}\right)\frac{\partial }{\partial x}{f}_{\mathrm{b}}=\frac{1}{2}\frac{\partial }{\partial {\mu }_{|\mathrm{w}}}\left[\left(1-{\mu }_{|\mathrm{w}}^2\right){\nu }_{scatt|\mathrm{w}}\frac{\partial }{\partial {\mu }_{|\mathrm{w}}}{f}_{\mathrm{b}}\right].$$

$$\left(1-{\gamma }_{\mathrm{w}}^2{\theta }_{|\mathrm{w}}^2\right)\frac{\partial }{\partial {\nu }_{scatt|w}x}{f}_{\mathrm{b}}=$$

$${\gamma }_{\mathrm{w}}^3\frac{1}{{\gamma }_{\mathrm{w}}{\theta }_{|\mathrm{w}}}\frac{\partial }{\partial \left({\gamma }_{\mathrm{w}}{\theta }_{|\mathrm{w}}\right)}\left[{\gamma }_{\mathrm{w}}{\theta }_{|\mathrm{w}}\frac{\partial }{\partial \left({\gamma }_{\mathrm{w}}{\theta }_{|\mathrm{w}}\right)}{f}_{\mathrm{b}}\right].$$

$$f_{\mathrm{b}}\simeq \exp \left[-k{\nu }_{scatt|w}(p_{|\mathrm{w}})x\right]h_k\left({\gamma }_{\mathrm{w}}{\theta }_{|\mathrm{w}}\right),$$

$$f_{\mathrm{b}}(x,p_{|\mathrm{w}},{\mu }_{|\mathrm{w}})\simeq C(p_{|\mathrm{w}})e^{-4{\gamma }_{\mathrm{w}}^3{\nu }_{scatt|w}(p_{|\mathrm{w}})x-{\gamma }_{\mathrm{w}}^2{\theta }_{|\mathrm{w}}^2},$$

Full text

I Introduction

In recent decades, relativistic collisionless shock waves have emerged as outstanding phenomena of high energy astrophysics. As natural consequences of the powerful outflows associated with objects such as gamma-ray bursts, pulsar wind nebulae or active galactic nuclei, these shock waves seemingly dissipate with high efficiency the energy that drives them into powerlaws of relativistic charged particles, up to extreme energies, see e.g. Bykov et al. (2012); Sironi et al. (2015); Pelletier et al. (2017) for reviews. By now, the basic mechanism of particle acceleration at shock waves is well understood in the test particle limit, in which one neglects the backreaction of these suprathermal particles on the shock and its environment, e.g. Bednarz and Ostrowski (1998); Kirk et al. (2000); Achterberg et al. (2001); Lemoine and Pelletier (2003); Keshet and Waxman (2005); Lemoine et al. (2006); Niemiec et al. (2006). Much of the ongoing effort focuses on the nonlinear backreaction of accelerated particles on the magnetized turbulence, and how this turbulence controls in turn the acceleration and radiation processes. The acceleration of particles, the generation of turbulence and the structuring of the shock are indeed recognized as three inseparable aspects of collisionless shocks.

For relativistic shocks, the level of magnetization of the unshocked medium is expected to be so low – $\sigma\,\equiv\,(u_{\rm A}/c)^{2}\,\sim\,10^{-9}$ for the interstellar medium ($u_{\rm A}$ Alfvén four-velocity) – that the magnetized turbulence is believed to be entirely generated by micro-instabilities acting in the shock precursor, which arise from the interpenetration of the beam of suprathermal particles and the unshocked (background) plasma Medvedev and Loeb (1999); *1999ApJ...511..852G. In other words, the suprathermal particles self-generate the magnetized turbulence in which they themselves scatter to be accelerated, and possibly also radiate at late times.

The present paper belongs to a series in which we explore the physics of unmagnetized relativistic pair shock fronts, using detailed analytical developments that we compare to dedicated particle-in-cell (PIC) simulations Lemoine et al. (2019a). In Paper I Pelletier et al. (2019), we have argued that the microturbulence in the shock precursor is of a mostly magnetostatic nature in a frame noted $\mathcal{R}_{\rm w}$, which is found to move sub-relativistically with respect to the background plasma. In Paper II Lemoine et al. (2019b), we have developed a model of the heating and deceleration of the background plasma through its interaction with the microturbulence. Finally, Paper IV Vanthieghem et al. (2019) will address the evolution of the microturbulence in the shock precursor.

In this Paper III, we focus on the physics of the suprathermal particles. In particular, we calculate their scattering length in the microturbulence and derive their distribution function in the precursor. As in the other papers of this series, we confront these theoretical predictions to dedicated PIC simulations. The scattering length $l_{\rm scatt}$ is a crucial quantity because it determines the time it takes to complete a Fermi cycle around the shock, hence the acceleration timescale, and hence the maximum energy of accelerated particles. Estimates of this acceleration timescale have been provided through test-particle simulations in a general model of the microturbulence Plotnikov et al. (2013) or through a direct measurement in long-timescale PIC simulations Sironi et al. (2013), yet a clear analytical determination of $l_{\rm scatt}$ is still missing. This calculation is not straightforward, in particular because of nontrivial dependencies on the anisotropy of the turbulence spectrum Achterberg et al. (2007). This length scale, defined in the shock rest frame, reveals a nontrivial scaling with $\gamma_{\infty}$, which plays an important phenomenological role as we argue in Sect. V.

In order to model the backreaction of the suprathermal particles on the precursor, one needs to determine the evolution of their distribution function. In the following we carry out such a calculation in both the steady-state and time-dependent regimes, the latter being most relevant to PIC simulations. This model also provides us a way to infer $l_{\rm scatt}$ from PIC simulations, which we compare with its theoretical value. Besides, those results help develop the scenario of heating and slowdown of the background plasma, presented in Paper II.

This paper is laid out as follows. Section II defines the main quantities used in this work; Sect. III presents a theoretical calculation of the scattering length scale in a microturbulence in relativistic motion with respect to the shock front; Sect. IV discusses the properties of the suprathermal particle population, how it is distributed in the precursor, while clearly defining the stationary and time-dependent regimes; finally, Sect. V summarizes our results and draws some phenomenological consequences. We use Gaussian cgs units with $k_{\rm B}=c=1$.

II Physical and simulation setups

The general setup is as follows: we consider an unshocked background plasma – quantities indexed with subscript p – inflowing from $+\infty$ towards $-\infty$ along the $x-$axis (shock normal), in the (lab) frame in which the shock lies at rest at $x\,=\,0$, noted $\mathcal{R}_{\rm s}$. The Lorentz factor $\gamma_{\rm p}$ of the background plasma, before it enters the shock precursor, is written $\gamma_{\infty}$, its density $n_{\infty}$, etc. Quantities indexed with ${}_{|\rm w},\,_{|\rm p},\,_{|\rm d}$ are respectively defined in the turbulence frame $\mathcal{R}_{\rm w}$ (see Lemoine et al. (2019a); Pelletier et al. (2019); Lemoine et al. (2019b) and below), in the background plasma rest frame and in the downstream rest frame, which coincides with the reference frame of our PIC simulations. Quantities that are not indexed are by default considered in the shock frame, unless they are proper.

The shock precursor is defined as that region (of finite extent $\ell_{\rm prec}$) at positive values of $x$ in which the background plasma (indexed with p) interpenetrates a gas of suprathermal particles (indexed with b). The latter correspond to particles that have been reflected on the shock, or accelerated through a Fermi-like process by multiple interactions around the shock front. As a result of the shock-crossing energy and momentum conserving equations in the fluid limit Blandford and McKee (1976), the typical momentum of shock-heated particles is $\sim\,\gamma_{\infty}m$, and that of suprathermal particles is larger by a factor of the order of a few to ten. It is expected that these suprathermal particles exhibit a Maxwell-Jüttner thermal distribution that turns over into a power-law of index $-s$ at large values of the momentum.

As discussed in various references, e.g. Medvedev and Loeb (1999); *1999ApJ...511..852G; Wiersma and Achterberg (2004); Lyubarsky and Eichler (2006); Achterberg and Wiersma (2007); Achterberg et al. (2007); Lemoine and Pelletier (2010, 2011); Rabinak et al. (2011) for theoretical considerations, Kato (2007); Spitkovsky (2008a, b); Martins et al. (2009); Keshet et al. (2009); Nishikawa et al. (2009); Sironi et al. (2013) for PIC simulations, the interpenetration of the beam of suprathermal particles and the background plasma gives rise to current filamentation instability (CFI), which generates an electromagnetic turbulence on plasma length scales $c/\omega_{\rm p}$ Weibel (1959); Davidson et al. (1972); Bret et al. (2004, 2008, 2010a, 2010b). In the unmagnetized limit, the dominant mode appears to be the transverse CFI, which forms current density filaments oriented along the shock normal, surrounded by (toroidal) magnetic fields and (radial) electric fields. In the precursor of a relativistic shock, the interaction between the suprathermal beam and the background plasma is strongly asymmetric. In the $\mathcal{R}_{\rm s}$ frame, the incoming background plasma is dense, (initially) cold and fast, moving with Lorentz factor $\gamma_{\rm p}$ and forming a tight beam in momentum space with opening angle $\lesssim\,1/\gamma_{\rm p}$, while the suprathermal particle population is tenuous, relativistically hot, slow and roughly isotropic. As discussed in previous papers of this series Lemoine et al. (2019a); Pelletier et al. (2019); Lemoine et al. (2019b), there exists a particular frame $\mathcal{R}_{\rm w}$, in which the filamentation instability gives rise to essentially magnetostatic turbulence modes, which move at subrelativistic velocities relative to the background plasma. In this frame, the suprathermal particles are highly energetic, forming a beam with typical momentum $\gamma_{\infty}\gamma_{\rm p}$ and opening angle $\sim\,1/\gamma_{\rm p}$, while the background plasma is (at least initially) cold and nearly isotropic, drifting at subrelativistic velocity. Due to their relatively low momenta, the plasma particles are mainly trapped in the magnetic structures in $\mathcal{R}_{\rm w}$. By contrast, the suprathermal particles are scattered on scales much larger that the size of those structures, and so propagate in an essentially ballistic manner.

It is thus important to carefully distinguish these populations. In our numerical simulations, we define the background plasma particles as those particles that propagate with negative $x-$velocity continuously since their injection into the precursor, i.e., their $x-$velocity has never changed sign. Oppositely, the suprathermal (or shock-reflected) particles are those moving with positive $x-$velocity, independently of the number of turn-arounds. Figure 1 shows the phase space $p_{x}/m$ vs $x$ of the total and background plasma populations in one such simulation (with $\gamma_{\rm\infty|d}=100$).

Our 2D3V (2D in configuration space, 3D in momentum space) numerical simulations are performed using the finite-difference time-domain PIC code calder. They describe the self-consistent formation and evolution of an unmagnetized, relativistic collisionless shock wave in an electron-positron pair plasma. The shock is initialized by injecting the background plasma from the right-hand side and having it reflect specularly on a conducting wall, as in Spitkovsky (2008b). The plasma is injected with proper temperature $T_{\infty}/m=10^{-2}$ and Lorentz factor $\gamma_{\infty|\rm d}\,=\,10$ in the simulation frame, which coincides with the downstream rest frame. Hence the Lorentz factor of the background plasma relative to the shock front is $\gamma_{\infty}\,=\,17$. Both simulations initially employ 10 macro-particles per cell and per species, with a cell size $\Delta x\,=\,\Delta y\,=\,0.1\,\omega_{\rm p}^{-1}$. The PIC code makes use of the Godfrey-Vay filtering algorithm Godfrey and Vay (2014) and the Cole-Karkkainen finite difference field solver Cole (1997, 2002); Kärkkäinen et al. (2006) in order to mitigate the Čerenkov instability while preserving a large time step $\Delta t=0.99\Delta x$.

The asymmetry of the beam-plasma interaction in the precursor is clearly illustrated in Fig. 2, which displays the density variations of both populations as a result of the CFI in the $\gamma_{\rm\infty|d}\,=\,100$ case. The background plasma (middle) is seen to develop mildly nonlinear current filaments in response to the magnetic modulations (top), while the suprathermal particles (bottom) show very weak fluctuations only.

The present paper focuses on the characteristics of the suprathermal particle population. We first discuss the scattering length of these particles in the microturbulence – thereby confirming, in particular, that it is much larger than the typical transverse size of a filament – then consider their distribution function in the precursor. As shown in earlier papers of this series Lemoine et al. (2019a); Pelletier et al. (2019); Lemoine et al. (2019b), a key quantity is the kinetic pressure $p_{\rm b}$, parametrized in terms of the incoming momentum flux through $\xi_{\rm b}\,\equiv\,p_{\rm b}/\left(\gamma_{\infty}^{2}n_{\infty}m\right)$.

III Scattering length scale

The physics of scattering of suprathermal particles in the downstream flow of an unmagnetized, relativistic collisionless shock is relatively easy to understand. As PIC simulations have demonstrated Spitkovsky (2008b, a); Chang et al. (2008); Martins et al. (2009); Keshet et al. (2009); Sironi and Spitkovsky (2011); Sironi et al. (2013), the downstream turbulence is essentially magnetostatic in $\mathcal{R}_{\rm d}$ and distributed on short length scales $\lambda_{\delta B}\,\sim\,\mathcal{O}\left(10\,\omega_{\rm p}^{-1}\right)$. A particle thus suffers a typical deflection by an angle $\delta\theta\,\simeq\,\pm\lambda_{\delta B}/r_{\rm g}$ upon crossing a coherence length scale $\lambda_{\delta B}$, $r_{\rm g}$ denoting the gyration radius in the average magnetic field. There follows the pitch-angle scattering frequency $\nu_{\rm scatt}\,\simeq\,\delta\theta^{2}/\lambda_{\delta B}\,\simeq\,\lambda_{\delta B}/r_{g}^{2}$. Hence, the scattering length, i.e. the length beyond which the acquired deflection becomes of the order of unity, is $l_{\rm scatt}\,\simeq\,r_{\rm g}^{2}/\lambda_{\delta B}$. That $l_{\rm scatt}(p)\,\propto\,p^{2}$ is expected because the turbulence wavelength is much smaller than the gyroradius of accelerated particles, see Pelletier et al. (2009); Lemoine and Pelletier (2010); Plotnikov et al. (2011, 2013).

In the shock precursor, this scattering issue needs a careful analysis because: (i) the beam of suprathermal particles is strongly anisotropic in the $\mathcal{R}_{\rm w}$ frame in which the turbulence is of a mostly magnetic nature; (ii) this frame is in relativistic motion with respect to the $\mathcal{R}_{\rm s}$ shock rest frame; (iii) the turbulence itself is strongly anisotropic.

The strong relativistic beaming of the suprathermal particles thus entails that their scattering frequency transforms nontrivially from $\mathcal{R}_{\rm w}$ to $\mathcal{R}_{\rm s}$, namely, $\nu_{\rm scatt}\simeq 4\gamma_{\rm w}^{3}\nu_{\rm scatt|w}$ (see Sec. IV.1). The anisotropy of the turbulence also plays an important role. For instance, describing the Weibel turbulence as a collection of infinitely long filaments oriented along the $x-$axis leads to the conservation of the conjugate canonical momentum associated to $\dot{x}$, thus precluding pitch-angle diffusion of the suprathermal particles Achterberg et al. (2007).

In the following, we first compute the correlation time of the random force experienced by a suprathermal particle, then its scattering length in the framework of an extended quasilinear theory, taking into account the anisotropy of the turbulence, its growth in the precursor and the relative motion between the $\mathcal{R}_{\rm w}$ frame and the shock front.

III.1 The correlation time of the random force

From a statistical mechanics perspective, the correlation time of the random force that a particle suffers along its trajectory is a crucial quantity. It is defined as

[TABLE]

with $\bm{\delta}\mathbf{F}(t)\,=\,q\left(\bm{\delta}\mathbf{E}+\bm{\beta}\times\bm{\delta}\mathbf{B}\right)$ the random Lorentz force exerted at time $t$, at position $\bm{x}(t)$ of the particle. Formally, $t$ and $\tau$ in Eq. (1) above represent the time in the lab frame along the particle trajectory; however, for suprathermal particles, it is a good approximation to use for $\bm{x}(t)$ a straight line trajectory, because the scattering length scale is expected to be much larger than the correlation length, as a result of the small-scale nature of the turbulence.

In the $\mathcal{R}_{\rm s}$ frame, $\bm{\delta}\mathbf{E}\,=\,-\bm{\beta_{\rm w}}\times\bm{\delta}\mathbf{B}$, hence

[TABLE]

The prefactor depends on quantities that vary slowly on $t_{\rm corr}$ timescales, which therefore drop out when taking the ratio in Eq. (1). Using $\bm{\delta}\mathbf{B}\,=\,\gamma_{\rm w}\bm{\delta}\mathbf{B_{|\rm w}}$, with $\bm{\delta}\mathbf{B_{|\rm w}}$ defined in the $\mathcal{R}_{\rm w}$ frame, we decompose the latter in $\mathcal{R}_{\rm w}$ plane waves, to obtain

[TABLE]

with $\mathcal{N}\,\equiv\,\int{\rm d}\bm{k^{\prime}}\,{{\mathcal{S}}_{|\rm w}}(\mathbf{k^{\prime}})$, and $\Delta x_{|\rm w}^{\alpha}$ the displacement of the particle along its world line during the time interval $\Delta t$. We also defined the power spectrum through

[TABLE]

The turbulence is assumed stationary in $\mathcal{R}_{\rm s}$, and we have extracted from $\delta B$ its spatial profile, taken in the form $\propto\exp(-2gx)$. This choice means that the fluctuation spectrum $S_{|\rm w}$ is normalized to a value of $\langle B_{\rm w}^{2}\rangle$ close to the shock front. Note that the spatial growth rate $g$, which is assumed independent of $\mathbf{k^{\prime}}$, is expressed in $\mathcal{R}_{\rm s}$ while the plane-wave expansion is performed in $\mathcal{R}_{\rm w}$.

The above decomposition of the magnetic field into a slowly evolving envelope $\exp(-g\,x)$ times a plane wave decomposition should be understood as a simplified description of the turbulence in the precursor; it assumes, in particular, that the power spectrum is preserved throughout the precursor while the turbulence magnetic energy grows. Note also that the quantity ${k^{\prime}}^{\alpha}\,{\Delta x_{|\rm w}}_{\alpha}$ is a Lorentz scalar, hence it can be expressed in any frame.

The quantity ${\Delta x_{|\rm w}}^{\alpha}$ represents the spacetime difference in $\mathcal{R}_{\rm w}$ between the trajectory at $t+\tau$ and that at $t$. If the particle is initially emitted with a pitch angle cosine $\mu$ along the shock normal, in the $\mathcal{R}_{\rm s}$ frame, then

[TABLE]

Note that, because it is expressed in $\mathcal{R}_{\rm s}$, $\mu$ is in principle determined at random between [math] and $1$. In the following, we approximate $\beta_{\rm w}\,\rightarrow\,-1$ when possible. We also note that the combination $-\omega^{\prime}{\Delta x_{|\rm w}}^{t}+k^{\prime}_{x}{\Delta x_{|\rm w}}^{x}$ can be rewritten $\gamma_{\rm w}(1+\mu)\left(1-v_{\omega^{\prime}}\right)\tau$, introducing the quantity $v_{\omega^{\prime}}\,=\,\omega^{\prime}/k^{\prime}_{x}$; $\omega^{\prime}$ is real by definition, since the growing part of the magnetic turbulence has been extracted previously.

For the purely transverse mode of the CFI, $\omega^{\prime}\,=\,0$. However, at finite $k_{x}^{\prime}$ (yet $k_{x}^{\prime}\,\ll\,k_{\perp}^{\prime}$), the cold fluid dispersion relation of the CFI yields $\omega^{\prime}\,=\,\beta_{\rm b|\rm w}k^{\prime}_{x}\left[1+\mathcal{O}(\xi_{\rm b})\right]$ Plotnikov et al. (2013), with $\beta_{\rm b|\rm w}$ the velocity of the beam in the $\mathcal{R}_{\rm w}$ frame. Then $v_{\omega}^{\prime}\,\approx\,1$, because $\beta_{\rm b|\rm w}\,\simeq\,1$ to an accuracy of order $1/\gamma_{\rm b|\rm w}^{2}\,\sim\,1/\gamma_{\infty}^{4}$. In a realistic shock precursor, where the filamentation instability reaches a nonlinear stage, and where oblique modes may contribute to shaping the turbulence, $\omega^{\prime}$ may not obey the above relation, yet we retain the scaling $\omega^{\prime}\,\propto\,k_{x}^{\prime}$ and discuss the influence of $v_{\omega^{\prime}}$ on the result.

We now define the resonance function

[TABLE]

To keep the integrals analytically tractable, we approximate the spectrum with a constant patch in $\mathbf{k^{\prime}}-$space, centered on $(0,\widehat{k^{\prime}_{\perp}})$ with extension $(\pm\Delta k^{\prime}_{x},\pm\Delta k^{\prime}_{\perp})$, describing a mostly transverse instability as it should:

[TABLE]

We note that $g\,=\,g_{|\rm w}/\gamma_{\rm w}$, since the growth time of the instability $g^{-1}$ in $\mathcal{R}_{\rm s}$ is $\gamma_{\rm w}$ times that experienced by the background plasma in $\mathcal{R}_{\rm w}$ Vanthieghem et al. (2019). Furthermore, $g_{|\rm w}\,\ll\,\widehat{k^{\prime}_{\perp}}$ for the filamentation instability, hence the real part $g\mu$ is a small quantity relative to ${k^{\prime}}^{\alpha}{\Delta x_{|\rm w}}_{\alpha}/\tau$.

Consider first the limit $g\,\rightarrow\,0$. Then the integral in Eq. (3) depends on the quantities $\Delta k^{\prime}_{x}$, $\widehat{k^{\prime}_{\perp}}$ and $\Delta k^{\prime}_{\perp}$ characterizing the spectrum as well as on $(1-v_{\omega^{\prime}})$. More specifically,

[TABLE]

with $\nu=(1-\mu)^{1/2}/(1+\mu)^{1/2}\,\approx\,1$. The importance of $v_{\omega^{\prime}}$ is thus clear. If $v_{\omega^{\prime}}$ is small compared to unity, then the condition specified on the rhs of the above equations amounts to whether $\Delta k_{x}^{\prime}\,<\,\Delta k_{\perp}^{\prime}/\gamma_{\rm w}$ or not, assuming $\widehat{k^{\prime}_{\perp}}\,\sim\,\Delta k_{\perp}^{\prime}$. In the limit $\gamma_{\rm w}\,\gg\,1$, the latter condition $\Delta k_{x}^{\prime}\,>\,\Delta k_{\perp}^{\prime}/\gamma_{\rm w}$ appears likely, hence one should expect $t_{\rm corr}\,\propto\,\left(\gamma_{\rm w}\Delta k^{\prime}_{x}\right)^{-1}$. There follows the expected result that (in the limit of small-angle deflections) a nonvanishing $\Delta k^{\prime}_{x}$ is required to ensure a finite $t_{\rm corr}$, and therefore pitch-angle diffusion.

However, if $v_{\omega^{\prime}}\,\simeq\,1-\epsilon$, as for the linear growth of the CFI, the former condition may hold. The spectrum is then such that it forbids the resonance of a particle with waves, preventing pitch-angle diffusion in this linearized limit. Obviously, any finite width to the dispersion relation, as characterized by the contribution in $g$ for instance, will lead to resonance broadening and permit pitch-angle diffusion. To first order in $g$, one obtains

[TABLE]

with $\alpha_{1}\,=\,\mu\ln\left[(\widehat{k^{\prime}_{\perp}}+\Delta k^{\prime}_{\perp})/(\widehat{k^{\prime}_{\perp}}-\Delta k^{\prime}_{\perp})\right]/(1-\mu^{2})\,\approx\,1$ and $\alpha_{2}\,=\,\pi\left[1+\mathcal{O}\left(g^{\prime}/\widehat{k^{\prime}_{\perp}}\,\Delta k^{\prime}_{x}/\Delta k^{\prime}_{\perp}\right)\right]/2\,\approx\,1$. We are particularly interested in the dependence of $t_{\rm corr}$ on $\gamma_{\rm w}$ and, interestingly, both limits leads to $t_{\rm corr}\,\propto\,(\gamma_{\rm w}\widehat{k^{\prime}_{\perp}})^{-1}$. The prefactor can be smaller or larger than unity, depending on which limit applies. Assuming for instance $\Delta k^{\prime}_{\perp}\,\approx\,\widehat{k^{\prime}_{\perp}}$, the first limit implies $\gamma_{\rm w}\widehat{k^{\prime}_{\perp}}t_{\rm corr}\,\sim\,\mathcal{O}\left(g_{|\rm w}/\widehat{k^{\prime}_{\perp}}\right)$, which is typically an order of magnitude smaller than unity in the linear growth phase of the CFI in the precursor of a relativistic shock. The second limit yields $\gamma_{\rm w}\widehat{k^{\prime}_{\perp}}t_{\rm corr}\,\sim\,\mathcal{O}\left(\widehat{k^{\prime}_{\perp}}/\Delta k^{\prime}_{x}\right)$, which is typically expected to be somewhat larger than unity. We will compare this prediction to measurements made in PIC simulations in the following.

Once the correlation time of the random force is known, one can estimate the scattering length (in $\mathcal{R}_{\rm s}$) by noting that, over $t_{\rm corr}$, the particle suffers a deflection of the order of $\pm ct_{\rm corr}/r_{\rm g}$, so that $l_{\rm scatt}\,\approx\,r_{\rm g}^{2}/(ct_{\rm corr})$.

III.2 Quasilinear estimate of $l_{\rm scatt}$

We now carry out a quasilinear calculation of the pitch angle diffusion coefficient of suprathermal particles (with $\gamma_{|\rm p}\gg 1$) in the shock rest frame. In order to keep track of conserved quantities in the possible limit of time-independent or $x-$independent turbulence (describing infinitely long filaments), we rely on a Hamiltonian formalism for the equations of motion. We first note that, for suprathermal particles, the canonical conjugate momentum $\pi^{\alpha}\,\equiv\,p^{\alpha}\,+\,q\,\delta A^{\alpha}$ coincides with the momentum $p^{\alpha}$ to a small error of order $\lambda/r_{\rm g}\,\ll\,1$, because the four-vector potential $|\bm{\delta}\mathbf{A}|\,\sim\,\lambda\delta B$ in order of magnitude. We thus use the approximation

[TABLE]

to describe the evolution of the pitch angle cosine $\mu\,\equiv\,p^{x}/p$ in the shock frame, as

[TABLE]

with the Hamilton equation

[TABLE]

In the above two equations, the conservation of $\mu$ along the particle trajectory indexed by the affine coordinate $s$ is manifest if the turbulence is both time and $x-$independent. Under standard assumptions, if neither of these conditions holds, the pitch angle may start to diffuse; then, the error associated to Eq. (10) is bound to decrease in time, so that our approximation will become more and more accurate. In effect, this error is bounded by the range of variation of the four-vector potential: $\Delta\mu\,\sim\,|\Delta\delta A|/\pi^{t}$; it thus remains fixed in time while the r.m.s. of the pitch angle cosine distribution increases through diffusion.

To simplify the notations, all unprimed variables are understood to be defined in $\mathcal{R}_{\rm s}$ in this section, while primed variables are defined in $\mathcal{R}_{\rm w}$. Substituting Eq. (12) in (11), and using (10) gives

[TABLE]

Here, an extra factor of $1/p$ has appeared because ${\rm d}s\,\equiv\,{\rm d}\tau/p$, with ${\rm d}\tau$ a time interval defined in $\mathcal{R}_{\rm s}$; furthermore, ${\delta A^{\prime}}^{\alpha}$, ${p^{\prime}}_{\alpha}$ are now primed variables. The partial derivatives are more conveniently expressed in terms of primed partial derivatives, and ${A^{\prime}}^{\alpha}$ can be decomposed in plane waves with polarization four-vectors ${{e^{\prime}}_{\mathbf{k^{\prime}}}}^{\alpha}$:

[TABLE]

where, as in Sec. III.1, the spatial $x-$profile of ${\delta A^{\prime}}^{\alpha}$ has been extracted from the plane wave decomposition. Let us stress again that $g$ represents the growth length scale in the shock rest frame, and that it does not depend on $\mathbf{k^{\prime}}$.

Consequently,

[TABLE]

where it is understood that $\omega^{\prime}$ is real, since the growing part has been extracted previously.

The statistical properties of the microturbulence in the $\mathcal{R}_{\rm w}$ can be approximated through the correlators

[TABLE]

In particular, for an anisotropic axisymmetric configuration with $k^{\prime}_{\perp}\,>\,k^{\prime}_{x}$, we can set

[TABLE]

The case of $S_{\parallel}(\mathbf{k^{\prime}})>S_{\perp}(\mathbf{k^{\prime}})$ corresponds to filaments elongated along $x$, as considered here. The power spectrum $S_{\parallel}$ of $\mathbf{\delta A^{\prime}}$ is related to the power spectrum $S_{|\rm w})$ of $\mathbf{\delta B^{\prime}}$, defined in Eq. (4), through $\mathcal{S}_{|\rm w}(\mathbf{k^{\prime}})\,=\,k_{\perp}^{\prime 2}\mathcal{S}_{\parallel}$, since

[TABLE]

We now approximate the trajectory as rectilinear, as in Eq. (5) and evaluate the evolution of pitch angle cosine over a time interval $\Delta t$ assumed much larger than the coherence time of the electromagnetic force $t_{\rm corr}$:

[TABLE]

where the response function $\mathcal{R}_{\mathbf{k^{\prime}}}$ defined in Eq. (6) appears through the time integration. We have approximated ${p^{\prime}}_{\alpha}\,{p^{\prime}}_{\beta}\langle{{e^{\prime}}_{\mathbf{k_{1}^{\prime}}}}^{\alpha}\,{{{e^{\prime}}_{\mathbf{k_{2}^{\prime}}}}^{\beta}}^{\star}\rangle\,\simeq\,p_{x}^{\prime 2}\mathcal{S}_{\parallel}\,(2\pi)^{3}\,\delta\left(\mathbf{k_{1}^{\prime}}-\mathbf{k_{2}^{\prime}}\right)$ because $\mathcal{S}_{\parallel}\,>\,S_{\perp}$ and $p_{x}^{\prime}/p_{\perp}^{\prime}\,\sim\,\mathcal{O}(\gamma_{\rm w})$.

Everywhere in the integral, we can take the simplifying limit $\beta_{\rm w}\,\simeq\,-1$. We also note that $\langle\delta B^{\prime 2}\rangle\,=\,\langle\delta B^{2}\rangle/\gamma_{\rm w}^{2}$, $\langle\delta B^{2}\rangle$ denoting the rms measured in the shock frame $\mathcal{R}_{\rm s}$. To compute the above integral, we use the same power spectrum as in Sec. III.1 and we pay attention to the lowest order term in $g^{\prime}/\widehat{k^{\prime}_{\perp}}$. From the definition $\nu_{\rm scatt}\,=\,\left\langle\Delta\mu^{2}\right\rangle/2\Delta t$, we eventually obtain

[TABLE]

with $\alpha_{3}\,=\,\mu(1+\mu)^{2}\,\ln\left[(\widehat{k^{\prime}_{\perp}}+\Delta k^{\prime}_{\perp})/(\widehat{k^{\prime}_{\perp}}-\Delta k^{\prime}_{\perp})\right]\,\approx\,1$ and $\alpha_{4}\,=\,\pi(1-\mu)(1+\mu)^{2}/2\,\approx\,1$. In both limits, one can verify that, to within a factor of the order of unity, $\nu_{\rm scatt}\,\simeq\,t_{\rm corr}/r_{\rm g}^{2}$, with $r_{\rm g}\,=\,p/(e\langle\delta B^{2}\rangle^{1/2})$, as expected.

To encompass both limits, we write:

[TABLE]

where $\overline{k^{\prime}}$ is a wavenumber of approximate value $\widehat{k^{\prime}}_{\perp}^{2}/g^{\prime}$ if the $g^{\prime}\rightarrow 0$ resonance with waves cannot be satisfied (corresponding to the first limit) or $\Delta k^{\prime}_{x}$ in the opposite case. Up to a numerical prefactor, which can be as large as an order of magnitude or so, we will assume in the following that $\overline{k^{\prime}}\,\sim\,\omega_{\rm p}$, so that the above leads to our estimate for the scattering length $l_{\rm scatt}\,=\,\nu_{\rm scatt}^{-1}$:

[TABLE]

with $\epsilon_{B}\,=\,\langle\delta B^{2}\rangle/\left(4\pi\gamma_{\infty}^{2}nm\right)$. Here, we have used the background plasma Lorentz factor $\gamma_{\rm p}$ as a proxy for $\gamma_{\rm w}$, which represents a good approximation, see Paper I Pelletier et al. (2019).

The above estimate of the scattering length can be understood in a simpler way if one omits the anisotropy of the turbulence. Consider a magnetostatic turbulence with typical wavenumber $k^{\prime}$ (in $\mathcal{R}_{\rm w}$). For beam particles of Lorentz factor $\gamma^{\prime}$ in $\mathcal{R}_{\rm w}$, we have $k^{\prime}r^{\prime}_{\rm g}\,\simeq\,(k^{\prime}/\omega_{\rm p})\epsilon_{B}^{-1/2}\gamma_{\infty}\gamma^{\prime}/\gamma_{\rm w}\,\gg\,1$. Therefore the particles suffer small-angle interactions each time they cross a coherence length $\sim k^{\prime-1}$ of the microturbulence; the correlation time of the force thus reads $\sim k^{\prime-1}$, and the angular scattering frequency $\nu_{\rm s}^{\prime}\,\sim\,k^{\prime}/(k^{\prime}r^{\prime}_{\rm g})^{2}\,\sim\,\omega_{\rm p}(k^{\prime}/\omega_{\rm p})^{-1}\epsilon_{B}\gamma^{\prime-2}\gamma_{\infty}^{2}/\gamma_{\rm w}^{2}$. This scattering frequency is defined in $\mathcal{R}_{\rm w}$, and to convert it to $\mathcal{R}_{\rm s}$, one needs to multiply it by $\gamma_{\rm w}^{3}$ – see Eq. (28) further below – while expressing $\gamma^{\prime}\,\simeq\,\gamma_{\rm w}\gamma$. This gives $\nu_{\rm scatt}\,\sim\,\omega_{\rm p}(k^{\prime}/\omega_{\rm p})^{-1}\epsilon_{B}\gamma_{\rm w}^{-1}(\gamma_{\infty}m/p)^{2}$, as obtained above. The origin of the $\gamma_{\rm w}$ factor in the scattering length thus results from the motion of $\mathcal{R}_{\rm w}$ relative to $\mathcal{R}_{\rm s}$, not from the anisotropy. The latter rather introduces the various possible values of $\overline{k^{\prime}}$.

IV Distribution in the precursor

In astrophysical sources, it is generally expected that, on the time scale that observations probe, the shock wave and the acceleration process have reached a steady state. This means, in particular, that leptons and hadrons have been accelerated up to their respective maximal energies, which are determined through the competition between the energy loss time scale or the age of the source and the characteristic acceleration time scale. Most PIC simulations, however, do not reach a steady state, because particle acceleration, when effective, appears as an unbounded process on the simulation time scale. To make a proper comparison between theoretical predictions and PIC simulations, or to extrapolate the results of these simulations to astrophysical objects, it thus becomes important to distinguish which regime, steady state or not, applies.

To see this, consider the transport of an accelerated particle of Lorentz factor $\gamma\,>\,\gamma_{\infty}$ in the precursor, as viewed in the shock rest frame. The penetration length scale of this particle is controled by its scattering length scale $l_{\rm scatt}(\gamma)$, so that, on general grounds, one expects the distribution function of those particles ${\rm d}N/{\rm d}\gamma{\rm d}x$ to fall off exponentially on length scales $x\,\gtrsim\,l_{\rm scatt}(\gamma)$. Hence, at a distance $x$, one typically finds particles with a Lorentz factor $\gamma$ such that $l_{\rm scatt}(\gamma)\,\sim\,x$. Since $l_{\rm scatt}(\gamma)$ is a growing function of $\gamma$, the larger its $\gamma$, the further away from the shock front the particle can propagate. Our particle of Lorentz factor $\gamma$ moves in a ballistic manner on timescales $t\,\ll\,l_{\rm scatt}$, but it starts to diffuse in the turbulence on timescales $t\,\gg\,l_{\rm scatt}$. Once the PIC simulation has reached a duration $t_{\rm max|d}$ such that it exceeds the acceleration time scale to produce particles of Lorentz factor $\gamma$, and such that the precursor extends well beyond the maximal distance where such particles can be found, one does not expect further evolution of the distribution function for those particles beyond time $t_{\rm max|d}$, hence the PIC simulation has reached a steady state for that Lorentz factor.

As the precursor extends at velocity $\simeq\,c$, and because the acceleration timescale is of the order of $l_{\rm scatt}$ in the shock frame, those two conditions amount to the same: $t_{\rm max|d}\,\gg\,l_{\rm scatt}(\gamma)$. Then, the region close to the shock front with $x\,<\,l_{\rm scatt}(\gamma)$, populated by particles of Lorentz factor $\lesssim\,\gamma$ can be considered in steady state, while particles of larger Lorentz factor $\gamma$ or in regions further away from the shock front evolve in a time-dependent manner as the simulation runs. In the following, we thus distinguish between these two limits.

IV.1 Stationary state

To derive the distribution function of accelerated particles in the steady state, one can use the kinetic equation given in Paper II Lemoine et al. (2019b), which describes the evolution of the distribution function in a mixed coordinate system, with space variables given in the shock rest frame, and momentum variables expressed in the Weibel frame, in which the turbulence is mostly magnetostatic. We simplify further this equation by neglecting the inertial force terms, corresponding to the limit in which the scattering length of the accelerated particles is much larger than the transition scale of the shock, over which significant slowdown of the Weibel frame occurs. This simplification allows us to express analytically the approximate suprathermal particle distribution function, but it prevents a detailed comparison with PIC simulations in regions in which the background plasma effectively slows down. Unfortunately, a closed-form solution which takes this deceleration into account does not seem at hand. Note that the deceleration enters both through the inertial correction, through the presence of $\gamma_{\rm w}$ in the equation, as well as through the dependence of the scattering frequency over $\gamma_{\rm w}$. In a 3D momentum space, the equation takes the form

[TABLE]

This equation has been previously solved in full generality in Ref. Kirk et al. (2000) through an expansion in eigenfunctions, in the limit $\gamma_{\rm w}\,\rightarrow\,\gamma_{\infty}$ and $\nu_{\rm scatt|\rm w}$ independent of $p$. Here we obtain an approximate solution including the dependence of $\nu_{\rm scatt|\rm w}$ on $p$.

This equation assumes $p_{|\rm w}^{t}\,\simeq\,p_{|\rm w}$, and $\mu_{|\rm w}=p_{|\rm w}^{x}/p_{|\rm w}$ represents the cosine of the angle of the particle momentum with the shock normal. Finally, the right-hand side models the stochastic force felt by the particle in the magnetostatic turbulence, as characterized by a scattering frequency $\nu_{\rm scatt|w}$ (“Weibel frame”). To simplify it further, we neglect the dependence of $\nu_{\rm scatt|w}$ on $\mu_{|\rm w}$ and use the small-angle approximation $\mu_{|\rm w}\,\simeq\,1-\theta_{|\rm w}^{2}/2$. The latter is a good approximation, because $|1-\mu_{|\rm w}|\,\lesssim\,1/\gamma_{\rm w}\,\ll\,1$ for ultrarelativistic particles traveling upstream. The above equation can be then rewritten as

[TABLE]

Assuming separation of variables,

[TABLE]

we obtain the simple solution, with $k\,=\,4\gamma_{\rm w}^{3}$,

[TABLE]

where $C(p_{|\rm w})$ is a function of $p_{|\rm w}$. This function can be determined by matching, on the shock front, the full solution in the upstream half-plane with that in the downstream half-plane, for a proper choise of boundary conditions, see e.g. Kirk et al. (2000). This procedure then determines $C(p_{|\rm w})$ as a power-law of index $s+2$. The relationship between $\nu_{\rm scatt|w}(p_{|\rm w})$ and the scattering length $l_{\rm scatt}(p)$, which we have defined earlier in the shock frame, is not trivial because of the relativistic motion of the microturbulence frame $\mathcal{R}_{\rm w}$ with respect to the shock front. The scattering length $l_{\rm scatt}$ is indeed defined as the length scale over which suprathermal particles are deflected by an angle of order unity in the $\mathcal{R}_{\rm s}$ frame, in a time interval $\Delta t\,=\,l_{\rm scatt}$; however, in the $\mathcal{R}_{\rm w}$ frame, this corresponds to deflection by an angle $\sim 1/\gamma_{\rm w}$, over a time scale $\Delta t_{|\rm w}\,\simeq\,\left(\gamma_{\rm w}^{2}\nu_{\rm scatt|\rm w}\right)^{-1}$. Accounting for time dilation $\Delta t_{|\rm w}\,=\,\gamma_{\rm w}(1+\mu)\Delta t$, one eventually obtains $\nu_{\rm scatt|\rm w}\,\sim\,\gamma_{\rm w}^{-3}l_{\rm scatt}^{-1}$. Alternatively, computing infinitesimal variations of pitch angle $\Delta\theta_{|\rm w}$ and $\Delta\theta$ in the respective $\mathcal{R}_{\rm w}$ and $\mathcal{R}_{\rm s}$ frames, respectively on timescales $\Delta t_{|\rm w}$ and $\Delta t$, one finds

[TABLE]

Then, the average of $\nu_{\rm scatt|\rm w}$ over an isotropic distribution in $\mu_{|\rm w}$ in the interval $[\beta,1]$ gives

[TABLE]

which matches the previous estimate. In the following, we re-absorb for convenience the prefactor times the factor $4$ appearing in front of $\nu_{\rm scatt|\rm w}$ into our expression for $\nu_{\rm scatt}\equiv l_{\rm scatt}$ (determined up to a prefactor) so that $\nu_{\rm scatt}$ and $f_{\rm b}\propto\exp\left(-\nu_{\rm scatt}x\right)$.

In terms of shock frame pitch-angle cosine, one also derives

[TABLE]

hence

[TABLE]

Using $C(p_{|\rm w})\,\propto\,p_{|\rm w}^{-s-2}$, expressing $p_{|\rm w}=\gamma_{\rm w}(1-\beta_{\rm w}\mu)p$ and approximating $\beta_{\rm w}\,\simeq\,-1$, one recast $f_{\rm b}$ as

[TABLE]

The above introduces a pivot momentum $p_{\rm m}$ and its corresponding scattering frequency in the shock frame $\nu_{\rm m}\,\equiv\,\nu(p_{\rm m})$; $p_{\rm m}$ characterizes here the minimum momentum of the beam distribution function at a distance $\nu_{\rm m}^{-1}$ away from the shock front.

One is usually interested in power-law solutions $C(p)\propto p^{-s-2}$; however, in making comparison to PIC simulations, one must keep in mind that the extent of the power-law in the simulation is rather small, of the order of one decade, which modifies the scalings of $\beta_{\rm b}$ and $\xi_{\rm b}$ below. Hence, we will also consider exponentially suppressed powerlaw-like solutions. To be complete, note that one can also derive the above solution directly in the shock frame, by writing the pitch-angle scattering operator in terms of shock frame coordinates $p$ and $\mu$, and taking the asymptotic limit $\beta_{\rm w}\,\rightarrow\,-1$. One then finds that $f_{\rm b}$ obeys the following equation:

[TABLE]

One can verify that the differential operator commutes with $p(1+\mu)=p_{|\rm w}/\gamma$, as expected since this operator is nothing but the pitch angle scattering operator in $\mu_{|\rm w}$ at constant $p_{|\rm w}$, and that $f_{\rm b}$ as given in Eq. (31) is a solution to the above equation.

The beam distribution function can be normalized through the parameter $\xi_{\rm b}$, which we recall is defined as the ratio of the pressure of suprathermal particles to the incoming asymptotic momentum flux at infinity Lemoine et al. (2019a); Pelletier et al. (2019); Lemoine et al. (2019b),

[TABLE]

At distances $x\,\gg\,\nu_{\rm m}^{-1}$, and in the limit in which the distribution function can be considered as an unbounded power-law in momentum of index $-s-2$, one finds a power-law behavior for $\xi_{\rm b}(x)$,

[TABLE]

and

[TABLE]

If $p_{\rm m}\,\sim\,\gamma_{\infty}m$, then an order of magnitude for $\nu_{\rm m}^{-1}$ is $\epsilon_{B}^{-1}\,\gamma_{\rm w}\,c/\omega_{\rm p}$, or about $10-10^{2}\,c/\omega_{\rm p}$ in the shock vicinity, where $\epsilon_{B}\,\sim\,0.01-0.1$ and $\gamma_{\rm w}$ is of the order of a few.

Figure 3 plots the dependence $\xi_{\rm b}(x)$ observed in our two PIC simulations for shock Lorentz factors $\gamma_{\infty|\rm d}\,=\,10$ and $100$. The overlaid law $\xi_{\rm b}(x)\,\propto\,x^{-0.4}$ (dashed line) confirms that $\xi_{\rm b}(x)$ indeed follows a power-law scaling over the first $\sim\,300\,c/\omega_{\rm p}$, before turning over into an exponentially suppressed behavior. This length scale depends directly on the integration time of the simulation, $t_{\rm max|d}=3600\,\omega_{\rm p}^{-1}$ for $\gamma_{\infty|\rm d}=10$ and to $t_{\rm max|d}=6900\,\omega_{\rm p}^{-1}$ for $\gamma_{\infty|\rm d}=100$. The slight difference between the observed spatial power-law and that predicted above, $\xi_{\rm b}\,\propto\,x^{-0.1}$ for $s=2.2$, may be attributed to the difference between the spectrum of accelerated particles in the simulation and a pure powerlaw, and/or possibly to the fact that $\gamma_{\rm p}$ is evolving in this region in the numerical simulation. Repeating the above calculation of $\xi_{\rm b}$ for such a power-law with exponential suppression an order of magnitude above the injection threshold, one indeed finds a steeper power-law for $\xi_{\rm b}$, which transits into an exponentially suppressed dependence further away. Note that the theoretical profiles that we have derived ignore the evolution of $\epsilon_{B}$, which, although slow, may impact further the power-law behavior.

For reference, we plot in Fig. 4 the all-particle spectra ${\rm d}N/{\rm d}\gamma_{|\rm d}$ for our two reference PIC simulations, as integrated over the transverse dimension and over a box of length $300\,c/\omega_{\rm p}$ along $x$, centered at various positions $x_{|\rm d}$ as indicated. This figure shows how the spectrum evolves in the precursor. In particular, the expected (in 3D momentum space) powerlaw ${\rm d}N/{\rm d}\gamma_{|\rm d}\,\propto\,\gamma_{|\rm d}^{-2.2}$ Bednarz and Ostrowski (1998); Kirk et al. (2000); Achterberg et al. (2001); Lemoine and Pelletier (2003); Keshet and Waxman (2005); Lemoine et al. (2006); Niemiec et al. (2006); Sironi et al. (2015) is recovered in the steady-state regime, in the near precursor and downstream.

The 2D distribution function, which is needed for a proper comparison to PIC simulations, can be obtained by replacing the scattering operator according to: $\partial_{\mu_{|\rm w}}(1-\mu_{|\rm w}^{2})\partial_{\mu_{|\rm w}}\,\rightarrow\,\partial_{\theta_{|\rm w}}\partial_{\theta_{|\rm w}}$. One then finds a similar distribution function, up to the substitution $(1-\mu)/(1+\mu)\,\rightarrow\,(1/2)(1-\mu)/(1+\mu)$ in the exponential, and of course $-s-2\,\rightarrow\,-s-1$. The power-law behavior for $\xi_{\rm b}(x)$, however, remains unchanged.

Returning to 3D, the number density of suprathermal particles,

[TABLE]

follows a power-law $n_{\rm b}\,\propto\,(\nu_{\rm m}x)^{(1-s)/2}$ at large distances compared to $\nu_{\rm m}^{-1}$. Finally, the beam bulk velocity,

[TABLE]

evolves from $\beta_{\rm b}\,\simeq\,-0.27$ at the shock front $x\,\rightarrow\,0$, to $\beta_{\rm b}\,\rightarrow\,0$ at large distances $x\,\gg\,\nu_{\rm m}^{-1}$. The suprathermal particles are not fully isotropic in the downstream or shock front frame, but their average velocity in the shock frame remains well subrelativistic.

IV.2 Time-dependent regime

It is of interest to study the distribution function of suprathermal particles in a time-dependent regime in order to carry out a meaningful comparison to PIC simulations. Such a comparison will notably provide a direct estimate of the scattering length of particles as a function of their energy, which can then be compared to our quasilinear estimate Eq. (22). To approximate this distribution function in the time-dependent regime, $t\,\ll\,l_{\rm scatt}$, To this goal, we assume that particles move along straight lines but can at any time interact to suffer a deflection into the opposite half-space of pitch angle – i.e., those with $\mu\,>\,0$ are deflected in the half-space $\mu\,<\,0$ and vice versa – with a mean waiting time of $l_{\rm scatt}$. We suppose that the shock front injects ${\rm d}\dot{N}/{\rm d}p{\rm d}\mu$ particles (per unit transverse area of the shock front) per unit time, momentum interval and pitch angle cosine interval. We also distinguish between forward- and backward-moving beam particles: the backward-moving particles result from the deflection of forward-moving particles. Since about half or more of the particles come back to the shock front after experiencing only one interaction, we neglect the possibility of multiple interactions in this time-dependent regime.

The forward-moving beam particles, with distribution function $f_{\rm b>}(x)$, thus correspond to the injected population that has not experienced any deflection up to distance $x$, and reads

[TABLE]

since $x/\mu$ indicates the time spent since injection at the shock front. Integrating over pitch-angle cosine, and assuming isotropic injection at the shock front, we derive the position dependent density

[TABLE]

keeping in mind that $l_{\rm scatt}$ also depends on $p$. The logarithmic divergence at $x\,\rightarrow\,0$ is an artefact that results from our assumption of an infinitely thin shock front; in the following, we will regularize it as ${\rm d}n_{\rm b>}/{\rm d}p{\rm d}x\,\approx\,{\rm d}\dot{N}/{\rm d}p$ as $x\,\ll\,l_{\rm scatt}$.

The backward-moving particles, although fewer in number, play a special role, as will be shown in the following. Their distribution can also be obtained assuming straight-line trajectories. Consider such a particle at $x$ at $t$, with momentum $p$ and pitch-angle cosine $\mu$. In our approximation, this particle results from the deflection of a forward-moving particle at some point $x_{1}\,\geq\,x$, and some time $t_{1}\,\leq\,t$, which itself was emitted by the shock front at some time $t_{0}\,\leq\,t_{1}$, with pitch-angle cosine $\mu_{0}$ and momentum $p_{0}$. For simplicity, we neglect the order of unity energy gain experienced by the particle suffers during its deflection and set $p_{0}\,=\,p$. If $G(x,t,\mu;x_{1},t_{1},\mu_{1})$ denotes the propagator representing rectilinear propagation without interaction, from coordinates $(x_{1},t_{1},\mu_{1})$ to $(x,t,\mu)$, then the density of backward-moving particles can be written

[TABLE]

Here, $P(\mu_{1<};\mu_{1>})$ denotes the probability of deflecting particle with incoming $\mu_{1>}$ into $\mu_{1<}$ upon interaction; we use $P(\mu_{1<};\mu_{1>})\,=\,\Theta\left[-\mu_{1<}\,\mu_{1>}\right]$. For ballistic transport,

[TABLE]

Assuming isotropic injection at the shock front, we then obtain

[TABLE]

since $t\,\ll\,l_{\rm scatt}$ by assumption. Note that the coordinates are expressed in the shock frame; in terms of downstream frame coordinates, which are more appropriate for a direct comparison with the PIC simulation, the above becomes

[TABLE]

with $\beta_{\rm d}$ the velocity of the downstream relative to the shock front; in a 2D PIC simulation, $\beta_{\rm d}\,\simeq\,-1/2$. The quantity in the numerator represents the difference between the total precursor length and the distance between the shock and the particle, in the simulation rest frame. It thus indicates that the backward-moving particle density vanishes at the tip of the precursor, then increases linearly as one nears the shock front. Quite interestingly, the above dependence of the backward-moving suprathermal particles distribution on distance offers a direct means to infer the scattering length scale of these particles from a PIC simulation. By contrast, this scattering length scale does not affect the distribution of forward-moving suprathermal particles in the time-dependent regime, since they propagate along nearly straight lines. It does control the distribution of the suprathermal particle population in the stationary state, but, as discussed above, this stationary state is limited to the lowest momenta range for a PIC simulation of reasonable duration, and is thus sensitive to the details of the turbulence in the vicinity of the shock transition.

The two panels of Fig. 5 plot the quantity $(p/p_{\rm m})^{-2}\Delta x{\rm d}N_{>}/{\rm d}N_{<}$, as measured in our PIC simulations, with $p$ the momentum of a particle, $\Delta x\,=\,x_{\rm max}-x$ the distance between the tip of the accelerated particle population and $x$, and ${\rm d}N_{>},{\rm d}N_{<}$ the ratio between the populations of forward- to backward-moving particles at $p$ and $x$. According to Eq. (43), this provides a direct estimate of the scattering length divided by $(p/p_{\rm m})^{2}$. In these two panels, the color code indicates the position $x$ where the estimate is taken. At large distances to the shock (corresponding to yellow/red colors), and at large momenta, where the assumed time-dependent regime should hold, all curves nearly lie on top of each other, indicating a coherent value of the scattering length. Furthermore, these curves do not depend on $p$, indicating that the scattering length does indeed scale as $p^{2}$. The gray band shows the value corresponding to the theoretical estimate given in Eq. (22), with a width corresponding to an uncertainty of a factor $3$ on either side. The satisfactory agreement obtained for all (large) values of $p$, (large) values of $x$ and the different values of $\gamma$ suggests that the simple formula (22) captures the leading dependencies of the scattering length in the shock precursor.

V Discussion – Conclusions

In this third paper of a series dedicated to the physics of unmagnetized, relativistic collisionless pair shocks, we have discussed the characteristics of the suprathermal population in the shock precursor. In particular, we have provided a theoretical estimate of the scattering length of these particles on the microscopic “Weibel” like turbulence, taking proper account of the anisotropy of this turbulence, its growth in the precursor and its relativistic motion with respect to the shock front, see Lemoine et al. (2019a); Pelletier et al. (2019); Lemoine et al. (2019b). We have obtained the formula: $l_{\rm scatt}(p)\,\approx\,\epsilon_{B}^{-1}\gamma_{\rm w}(p/p_{\rm m})^{2}\,c/\omega_{\rm p}$, with $p_{\rm m}\,=\,\gamma_{\infty}m$ the typical injection momentum of the accelerated particle distribution and $\gamma_{\rm w}$ the Lorentz factor of the turbulence relative to the shock front. We have compared this value to values extracted from dedicated large-scale PIC simulations and found satisfactory agreeement, both in terms of the prefactor and of the exponent of $p$. Furthermore, we have characterized the distribution of the suprathermal particles in the shock precursor and deduced the spatial profile of their momentum flux. Here again, reasonable agreement is found when making comparison with PIC simlulations.

A comment is worthwhile on the $\gamma_{\rm w}\simeq\gamma_{\rm p}$ prefactor in $l_{\rm scatt}(p)$. This prefactor, which derives from the relativistic motion between the “Weibel frame” and the shock frame (in which the suprathermal particle beam is roughly isotropic), bears important consequences for the phenomenology of such shock waves. This is so because $l_{\rm scatt}(p)$ sets the typical value of the residence time in the $\mathcal{R}_{\rm s}$ frame, hence the typical value of the acceleration timescale. The residence time in the downstream plasma is expected to be $\gamma_{\infty}$ times shorter, because the downstream flow is not relativistic relative to the shock front. The acceleration timescale has been directly measured in long-timescale PIC simulations Sironi et al. (2013). At present times, however, such simulations can only probe the early development of the power-law of accelerated particles, because $t_{\rm acc}(p)\,\propto\,p^{2}$ implies that the maximum energy $p_{\rm max}\,\propto\,t^{1/2}$. Furthermore, most of the accelerated particles in such simulations have gained energy in the shock vicinity. This is clearly seen in Fig. 1: the largest extent of the spectrum in $p_{x}$ of the accelerated particles is reached in the near precursor. In our model, this behavior can be ascribed to the scattering length increasing with $\gamma_{\rm w}$, and hence with the penetration depth into the precursor (since $\gamma_{\rm w}$ drops from $\gamma_{\infty}$ at the tip of the precursor down to $\sim 1$ at $x\,\sim\,0$), which renders grazing Fermi orbits more likely in a given simulation time. This also implies that current PIC simulations fail to access the acceleration timescale in the far precursor that controls the generation of the highest-energy particles in a realistic setting.

Most likely, the transport of these highest-energy particles in the far upstream would then become dominated by an external magnetic field, even if as weak as that of the interstellar medium, because the microturbulence implies an acceleration timescale $t_{\rm acc}\,\propto\,p^{2}$, while a regular magnetic field guarantees $t_{\rm acc}\,\propto\,p$. Hence, in spite of the inefficient scattering implied by the relativistically moving “Weibel” turbulence, the acceleration timescale in the interstellar magnetic field remains short enough to ensure that synchrotron GeV photons can be produced in the early phases of a highly energetic gamma-ray burst afterglow Plotnikov et al. (2013).

Acknowledgements.

We acknowledge financial support from the Programme National Hautes Énergies (PNHE) of the C.N.R.S., the ANR-14-CE33-0019 MACH project and the ILP LABEX (under reference ANR-10-LABX-63) as part of the Idex SUPER, and which is financed by French state funds managed by the ANR within the ”Investissements d’Avenir” program under reference ANR-11-IDEX-0004-02. This work was granted access to the HPC resources of TGCC/CCRT under the allocation 2018-A0030407666 made by GENCI. We acknowledge PRACE for awarding us access to resource Joliot Curie-SKL based in France at TGCC Center.