Towards understanding the physics of collisionless relativistic shocks

TL;DR

This paper explores the physical mechanisms behind collisionless relativistic shocks, emphasizing the role of filamentation instabilities and presenting a Monte Carlo model for particle acceleration that extends beyond PIC simulation capabilities.

Contribution

It provides a theoretical analysis of shock physics and introduces a non-linear Monte Carlo model for particle acceleration in relativistic shocks.

Findings

Filamentation instabilities seed particle scattering and preheat the plasma.

The Monte Carlo model captures particle back-reaction effects beyond PIC simulation limits.

Insights into shock dissipation and particle acceleration mechanisms.

Abstract

Relativistic astrophysical collisionless shocks represent outstanding dissipation agents of the huge power of relativistic outflows produced by accreting black holes, core collapsed supernovae and other objects into multi-messenger radiation (cosmic rays, neutrinos, electromagnetic radiation). This article provides a theoretical discussion of the fundamental physical ingredients of these extreme phenomena. In the context of weakly magnetized shocks, in particular, it is shown how the filamentation type instabilities, which develop in the precursor of pair dominated or electron-ion shocks, provide the seeds for the scattering of high energy particles as well as the agent which preheats and slows down the incoming precursor plasma. This analytical discussion is completed with a mesoscopic, non-linear model of particle acceleration in relativistic shocks based on Monte Carlo techniques. This Monte Carlo model uses a semi-phenomenological description of particle scattering which allows it to calculate the back-reaction of accelerated particles on the shock structure on length and momentum scales which are currently beyond the range of microscopic particle-in-cell (PIC) simulations.

Figures

Equations

$$(e_{\mathrm{d}}+p_{\mathrm{d}})u_{\mathrm{d}}^2+p_{\mathrm{d}}=n_{\mathrm{u}}m{c}^2{u}_{\mathrm{u}}^2,$$

$$(e_{\mathrm{d}}+p_{\mathrm{d}}){\Gamma }_{\mathrm{d}}{u}_{\mathrm{d}}=n_{\mathrm{u}}m{c}^2{\Gamma }_{\mathrm{u}}{u}_{\mathrm{u}}$$

$$w_{\mathrm{d}}=\frac{1}{3}\frac{1+3{\sigma }_{\mathrm{d}}}{1+{\sigma }_{\mathrm{d}}}.$$

$${\beta }_{\mathrm{d}}^2-(1+w_{\mathrm{d}}){\beta }_{\mathrm{u}}{\beta }_{\mathrm{d}}+w_{\mathrm{d}}=0.$$

$$p_{\mathrm{d}}=(1-w_{\mathrm{d}}){\gamma }_{\mathrm{sh}}^2{n}_{\mathrm{u}}m{c}^2.$$

$${\beta }_{\mathrm{dHT}}=-{\beta }_{\mathrm{u}}\frac{B_{\mathrm{u},\perp }}{B_{\mathrm{u},\parallel }}={\beta }_{\mathrm{sh}}{\gamma }_{\mathrm{sh}}\tan {\theta }_{B,|\mathrm{u}}.$$

$${\gamma }_{\mathrm{sh}}\sin {\theta }_{B,|\mathrm{u}}<1.$$

$${\ell }_{\mathrm{prec}|u}=\left(1-{\beta }_{\mathrm{sh}}\right)t_{u|u}\simeq \frac{{\omega }_{\mathrm{L},0}^{-1}}{2{\gamma }_{\mathrm{sh}}^3}$$

$$P_{\mathrm{cr}}={\xi }_{\mathrm{cr}}{\gamma }_{\mathrm{sh}}^2{n}_{\mathrm{u}}m{c}^2.$$

$$W_{\mathrm{em}}={\xi }_B{\gamma }_{\mathrm{sh}}^2{n}_{\mathrm{u}}m{c}^2.$$

$$\sigma \equiv \frac{B_{\perp }^2/4\pi }{{\gamma }_{\mathrm{sh}}^2{n}_{\mathrm{u}}m{c}^2}.$$

$$g_{inst|u}>{\gamma }_{\mathrm{sh}}{\omega }_{\mathrm{c}},$$

$$\sigma <{\gamma }_{\mathrm{sh}}^{-2}\frac{g_{inst|u}^2}{{\omega }_{\mathrm{pi}}^2}$$

$${\omega }_{\mathrm{pb}}={\left(\frac{4\pi n_{cr|u}{e}^2}{{\gamma }_{\mathrm{b}}m}\right)}^{1/2}$$

$$G=\frac{1-{\beta }_{\mathrm{rel}}{\mu }_1}{1-{\beta }_{\mathrm{rel}}{\mu }_2}.$$

$$⟨G⟩={\left(1+\frac{2}{3}|{\beta }_{\mathrm{rel}}|\right)}^2.$$

$$s=1-\frac{\ln \left(1-P_{\mathrm{esc}}\right)}{\ln ⟨G⟩}\simeq 1+3\left|\frac{{\beta }_{\mathrm{d}}}{{\beta }_{\mathrm{rel}}}\right|,$$

$$⟨1-P_{\mathrm{esc}}({\mu }_2)⟩⟨G({\mu }_1,{\mu }_2{)}^{s-1}⟩=1,$$

$$s\simeq 1-\frac{\ln \left(1-P_{\mathrm{esc}}\right)}{\ln ⟨G⟩}.$$

$${\gamma }_{\mathrm{cm}}=\frac{{\gamma }_{\mathrm{sh}}}{(1+{\xi }_{\mathrm{cr}}^2{\gamma }_{\mathrm{sh}}^2{)}^{1/2}}.$$

$$f_y=-{\rho }_{\mathrm{cr}}{E}_y\sim {\beta }_{\mathrm{sh}}{\xi }_{\mathrm{cr}}{n}_{\mathrm{u}}{\gamma }_{\mathrm{sh}}eB.$$

$$(-{\omega }^2+{\omega }_{\mathrm{p}}^2+k_{\perp }^2{c}^2)\delta A_{x|\mathrm{w}}$$

$${\omega }^2\delta {\rho }_{bg|w}$$

$$\delta j_x=\sum _{\alpha }{\rho }_{\alpha }\delta {\beta }_{x\alpha }+{\beta }_{\alpha }\delta {\rho }_{\alpha }$$

$${\rho }_{\alpha }\delta {\beta }_{x\alpha }=-\frac{{\omega }_{p\alpha }^2}{4\pi {\gamma }_{\alpha }^2}\delta A_x.$$

$${\omega }^2=-{\xi }_{\mathrm{cr}}{\beta }_{b|w}{\omega }_{\mathrm{p}}^2\frac{k_{\perp }^2{c}^2}{{\omega }_{\mathrm{p}}^2+k_{\perp }^2{c}^2}.$$

$$(\omega +{\xi }_{\mathrm{cr}}{k}_{x}c{)}^2=-{\xi }_{\mathrm{cr}}{\beta }_{b|w}{\omega }_{\mathrm{p}}^2\frac{k_{\perp }^2{c}^2}{{\omega }_{\mathrm{p}}^2+k^2{c}^2}.$$

$$\sum _s{\alpha }_s^{01}=0.$$

$$\delta E_{x|\mathrm{w}}\equiv -\frac{1}{c}{\partial }_t\delta A_{x|\mathrm{w}}\sim \frac{g_{inst|w}}{{\omega }_{\mathrm{pi}}}{B}_{|\mathrm{w}}.$$

$$\Delta {ϵ}_{|\mathrm{u}}={\beta }_{w|u}c\Delta p_{x|\mathrm{u}}.$$

$$k\dot{T}=\frac{1}{3}{\beta }_{w|u}c{F}_{d|u}.$$

Full text

1 Introduction

Relativistic collisionless shocks are extreme phenomena of plasma physics which form at the interface of the relativistic outflows of powerful astrophysical sources (e.g. gamma-ray bursts, pulsar wind nebulae and active galactic nuclei) and their environment. These shock fronts are not only outstanding dissipation agents, which convert into relativistic heat and disorder a well ordered relativistic flow, they also appear to produce quasi-power law spectra of accelerated particles up to very high energies through a relativistic variant of the Fermi process. These particles either escape the acceleration site and thus become cosmic rays, possibly very high energy cosmic rays, and/or they interact with ambient backgrounds to produce high energy photons and possibly neutrinos over a wide range of energies. In this regard, the physics of relativistic shocks is a key element in the arena of astroparticle physics and high energy astrophysics.

How such shocks form, how particle acceleration takes place, how magnetic turbulence is excited and to what level, remain fundamental open questions. Reviews on the physics, the numerical simulation and the phenomenology of relativistic collisionless shocks exist in the literature, which describe the progress achieved so far, see in particular Bykov et al. (2012b); Sironi et al. (2015); Marcowith et al. (2016). The present review takes a different stance and tries to describe analytically the physics of these phenomena to some level of detail. Due to the limited space, not all details could be presented nor could all regimes be investigated. The emphasis has been deliberately placed on unmagnetized relativistic shocks, which are representative, for example, of external shocks in gamma-ray bursts; other types of shocks are nonetheless briefly addressed to provide a more general picture.

The paper is laid out as follows; in Sec. 2, we provide essential notions on relativistic collisionless shock waves, discussing in particular the shock jump conditions and the particle kinematics; Sec. 3 desribes how a relativistic collisionless shock forms, depending on the ambient magnetization and shock velocity; Sec. 4 discusses the physics of the precursor, in which the superthermal particle population mixes with the background plasma, and where the latter is pre-heated and slowed down; Sec. 5 presents some phenomenological consequences and observational tests of the theory. Finally, we complete this analytical model with a description of Monte Carlo numerical simulation techniques in Sec. 6, which provides a complementary point of view on the physics of relativistic shocks. While Monte Carlo simulations must use as input some microphysics describing particle kinematics, which can be derived from the analytical theory presented earlier, they also provide a mesoscopic non-linear picture which is useful to extrapolate microphysical results to astrophysical scales of interest. Our main results are summarized in Sec. 8. Some details of a more technical nature are presented in Apps. A.1, A.2 and A.3.

2 The basics of Relativistic Shocks

Just as any other shock front, a relativistic shock in a plasma corresponds to a transition between two plasma states at a “front” which moves at velocity close to the velocity of light relative to some unshocked medium (the background plasma). Such a shock transition requires that the shock front moves at a velocity $v_{\rm sh}$ which is faster than the velocity at which causal disturbances can be transmitted in the background unshocked plasma. In the hydrodynamical limit, this velocity is $c_{\rm s}$, the sound velocity; if the background plasma is magnetized, however, various modes can propagate in the plasma, at different velocities, hence several types of shock waves can be envisaged; see Kirk and Duffy (1999) for a detailed discussion. This review will be interested in fast shock waves, which move faster than the fast magnetosonic velocity $\left(c_{\rm s}^{2}+v_{\rm A}^{2}\right)^{1/2}$, $v_{\rm A}$ corresponding to the Alfvén velocity in the background plasma. For simplicity, it will also be assumed that the shock front is planar and that it propagates into the background plasma along its normal.

Upon crossing the shock front, the ordered kinetic energy is converted into thermal energy mostly. A description of the physics of shock fronts is conveniently made in the so-called front frame that separates two homogeneous plasmas, namely the upstream plasma, the unshocked background plasma that inflows into the shock front, and the downstream plasma, which flows away from the shock at a subsonic velocity as a consequence of intense heating in the shock transition. In a collisional plasma, this transition is usually quite sharp, its size being measured in units of the mean free path of particles. In collisionless, non-relativistic plasmas, with an oblique magnetic field (i.e. meaning that the field lines have a finite angle with respect to the shock normal), the size of the transition is measured in term of the gyration radius of ions. In parallel shocks (i.e. with a magnetic field parallel to the shock normal), the transition is generally much more extended, because strong disturbances are generated and the absence of collisional dissipation allows them to travel over larger distances. In this case, the upstream and downstream plasmas close to the shock front need not be uniform. Nevertheless, on some scale, provided the downstream flow has been isotropized, the so-called Rankine-Hugoniot jump relations, which express the balance of the fluxes of matter, momentum and energy, are fulfilled.

In the next section, we recall these jump conditions and discuss the kinematics of particles around this shock front, which are peculiar to relativistic shock waves.

2.1 Shock jump conditions

In the case of a relativistic shock, it is convenient to derive these jump relations in the shock front frame by using the energy-momentum tensor $T_{\mu\nu}$ (e.g. Blandford and McKee, 1976; Double et al., 2004). A word of caution regarding notations: except otherwise noted, and except for proper quantities such as density $n$, energy density $e$ and pressure $p$ or enthalpy $h$, quantities will be expressed in the shock front rest frame; if expressed in the upstream or downstream rest frames, these quantities will be indexed with u or d respectively. Figure 1 illustrates the various frames used in this text and the corresponding flow velocities.

Assuming a plane front and a normal coordinate $x$, the relations express that the momentum flux $T_{x}^{x}$ and the energy flux $T_{x}^{t}$ keep the same values upstream and downstream. The momentum flux conservation reads:

[TABLE]

and the energy flux conservation reads:

[TABLE]

$u_{\rm u}$ (resp. $u_{\rm d}$) denote the $x-$component of the upstream (resp. downstream) plasma $4-$velocity in the shock frame, with corresponding Lorentz factor $\Gamma_{\rm u}$ (resp. $\Gamma_{\rm d}$). To make contact with our notations above, the shock velocity relative to the unshocked background plasma is of course $\beta_{\rm sh}\,=\,-\beta_{\rm u}$, with associated Lorentz factor $\gamma_{\rm sh}\,=\,\left(1-\beta_{\rm sh}^{2}\right)^{-1/2}$.

The above equations assume a cold upstream plasma. We introduce an equation of state $p_{\rm d}=w_{\rm d}e_{\rm d}$ for the downstream plasma; since the downstream plasma is relativistically hot, $w_{\rm d}=1/3$. Note that this does not require that the distribution function reflect thermal equilibrium, since this ratio holds for any isotropic distribution for which most of the particles are ultra-relativistic; in particular the distribution function can be composed of a thermal part with a high energy isotropized tail. One can generalize the above equation of state to that of a plasma containing a magnetic component set in the transverse directions to the shock normal (i.e. a magnetic field lying in the shock front plane), carrying a fraction $\sigma_{\rm d}$ of the energy density of the downstream plasma. Then,

[TABLE]

Now, the downstream velocity $\beta_{\rm d}$ is obtained by solving the following algebraic equation derived from the two conservation equations:

[TABLE]

For $\beta_{\rm u}^{2}\,\simeq\,1$, corresponding to the ultra-relativistic limit, the two roots are $1$ and $w_{\rm d}$. The downstream pressure of a strong relativistic shock is then

[TABLE]

From the mass flux conservation law, one derives easily that the apparent density in the shock front frame is amplified by a factor $w_{\rm d}^{-1}$, as measured in the front rest-frame; however the Lorentz transform from the background plasma rest-frame to the shock front frame amplifies the proper density by a factor $\gamma_{\rm sh}$, while the Lorentz transform from the shock frame to the downstream plasma rest-frame reduces the proper density by a factor $\gamma_{\rm d}\,\simeq\,3/(2\sqrt{2})$ in the hydrodynamic limit $\sigma_{\rm d}\,\ll\,1$. Thus the proper density eventually suffers an amplification by a compression ratio corresponding to the jump of the 4-velocity vector $r\equiv u_{\rm u}/u_{\rm d}=2\sqrt{2}\gamma_{\rm sh}$ (the latter equality applying to the hydrodynamical limit).

Combining this result with Eq. (5), one derives a mean energy per particle downstream of the shock of the order of $e_{\rm d}/n_{\rm d}\,\simeq\,\gamma_{\rm rel}$, with $\gamma_{\rm rel}\,=\,\gamma_{\rm u}\gamma_{\rm d}\left(1-\beta_{\rm u}\beta_{\rm d}\right)\,\simeq\,\gamma_{\rm sh}/\sqrt{2}$ the relative Lorentz factor between the up- and downstream frames. Superthermal particles share this mean energy per particle of order $\gamma_{\rm sh}$.

Let us consider now a mean field having some obliquity $\theta_{B|\rm u}$ (measured in the background frame). In the ideal magnetohydrodynamics (MHD) approximation, the normal ($x-$) component of the electric field that compensates the electromotive force $-\boldsymbol{\beta_{\rm u}}\times\boldsymbol{B}$ vanishes and thus the normal ($x-$) component $B_{x}$ of the magnetic field remains unchanged in the shock transition. In contrast, the transverse component of the electric field is continuous at the shock front, therefore $\beta_{\rm u}\boldsymbol{B_{{\rm u},\perp}}\,=\,\beta_{\rm d}\boldsymbol{B_{{\rm d},\perp}}$, the $\perp$ symbol denoting the component transverse to the shock normal; $B_{\perp}$ is thus amplified by a factor $w_{\rm d}^{-1}$ when crossing the shock front, in the shock front frame. However, taking into account the Lorenz factors involved when transforming from the flow co-moving frames to the front frame and vice-versa, the global amplification of the transverse magnetic field component is by the compression ratio $r\equiv u_{\rm u}/u_{\rm d}=2\sqrt{2}\gamma_{\rm sh}$ as for the density.

2.2 Particle kinematics

Because of this tremendous amplification of the transverse component of the magnetic field in ultra-relativistic shocks, relativistic shocks generically have a “superluminal” configuration, which means that one cannot find a Lorentz transform with a sub-luminal velocity along the shock surface that would cancel the electric field and have the flow moving along the magnetic field (the so-called de Hoffmann Teller frame, see de Hoffmann and Teller (1950) as well as Begelman and Kirk (1990)). Consider indeed a frame moving at the velocity $\boldsymbol{\beta_{\rm dHT}}$ along the direction of the projection $\boldsymbol{B_{\rm u,\perp}}$ on the shock surface; then in order to cancel the electric field, one needs to boost to a frame moving with a velocity $\boldsymbol{\beta_{\rm dHT}}$ such that $\boldsymbol{E}+\boldsymbol{\beta_{\rm dHT}}\times\boldsymbol{B_{\rm u,\parallel}}=0$; inserting $\boldsymbol{E}=-\boldsymbol{\beta_{\rm u}}\times\boldsymbol{B_{\rm u,\perp}}$, one obtains

[TABLE]

Since $\beta_{\rm dHT}$ must be smaller than one, the condition for having a Lorenz transform which wipes out the electric field in the front is very restrictive, because it requires

[TABLE]

When $\gamma_{\rm sh}\sin\theta_{B,|\rm u}>1$, the intersection point of a field line on the shock plane moves at a velocity larger than the velocity of light, hence the word “super-luminal” for qualifying such a shock. Quite generically, ultra-relativistic shocks are super-luminal and the matter flow always crosses the field lines and experiences an electric field.

As for superthermal particles, if they cross the shock front from downstream to upstream, it means that the $x-$component of their $3-$velocity $\beta_{x}$ is positive in the shock front frame, or equivalently larger than $\beta_{\rm sh}$ in the background plasma rest frame. Thus, in this frame, the pitch angle $\alpha$ of the superthermal particle motion with respect to the shock normal is such that $\sin\alpha<1/\gamma_{\rm sh}$; the particle distribution of superthermal particles is extremely anisotropic in the rest frame of the background plasma, although it can be considered as nearly isotropic in the shock rest frame.

The trajectory of superthermal particles upstream of the shock is a crucial point in the theory of relativistic shocks, since its extension determines the size of the precursor, which corresponds to that region in which the superthermal and background populations intermix. To a good approximation, one can consider that a superthermal particle, while traveling upstream with momentum $p_{|\rm u}$, is caught up by the shock front once its momentum is deflected by an angle $\gtrsim\,1/\gamma_{\rm sh}$, because the $x-$component of the particle $3-$velocity then falls below $\beta_{\rm sh}$ and the shock is trailing right behind the particle. In the presence of an oblique background magnetic field, this takes place on a timescale $t_{\rm u|u}\,\simeq\,\omega_{\rm L,0}^{-1}/\gamma_{\rm sh}$, with $\omega_{\rm L,0}\,\equiv\,eB_{\rm|u}/p_{\rm|u}$ the gyrofrequency in the background field. Hence the extension of the precursor (Achterberg et al., 2001; Milosavljević and Nakar, 2006; Pelletier et al., 2009):

[TABLE]

As $B_{\rm|u}\,\rightarrow\,0$, the size of the precursor becomes controlled by the scattering of the superthermal particles in the self-generated turbulence; this case will be discussed further on.

3 What is a Collisionless Relativistic Shock?

The concept of collisionless shock relies on the development of collective effects, including the Landau effect, that control the dissipation process, i.e. the transformation of the flux of kinetic energy through the generation of several kinds of entropy. The main entropy generation is the formation of a thermal distribution of ultra-relativistic temperature, with (proper) temperature downstream $T_{\rm d}\,\sim\,\gamma_{\rm sh}mc^{2}$ (in units in which $k_{\rm B}\,=\,1$) or merely the isotropization of the distribution function, as mentioned previously in the case of ultra-relativistic shocks.111Of course shocks do not have to be ultra-relativistic and can have any supersonic speed $\lesssim c$. This is a critical aspect of gamma-ray burst afterglows where the jet shock slows from ultra-relativistic to trans-relativistic to non-relativistic speeds. Trans-relativistic shocks are difficult to describe analytically but can be treated directly with PIC or Monte Carlo techniques. This isotropic distribution function can contain a superthermal component in the form of a quasi-power-law tail; this is what is expected for the generation of very high energy particles in astrophysics to account for gamma-ray radiation, for the production of very high energy cosmic rays, especially ultra-high energy cosmic rays, or even high energy neutrinos. The generation of a superthermal tail does not require more energy, of course; for a given global energy budget, a fraction of the incoming energy is transferred into the pressure of a cosmic ray component corresponding to a significant fraction of the thermal pressure:

[TABLE]

Note that the above equation neglects factors of order unity (see in particular Eq. 5). Typically a conversion factor $\xi_{\rm cr}$ of order 10 percent is found in numerical simulations and is expected to account for the observations.

Such particle acceleration cannot occur without the excitation of some electromagnetic turbulence that also takes a significant fraction of the incoming energy:

[TABLE]

How this turbulent component is generated will be discussed in Sec. 4. Of course a preexisting mean magnetic field or background turbulence may be present in addition to self-generated turbulence.

In the next section we describe two limiting cases: that of a magnetized shock with a strong mean ambient field, and a weakly magnetized shock where the self-generated turbulence dominates the mean field. We then present a general phase diagram for the physics of these relativistic collisionless shock waves and discuss the consequences for Fermi acceleration.

3.1 With a mean field

Early work on the theory and simulation of a relativistic shock with a transverse mean field in pair plasmas, with an obvious application to the termination shock of Pulsar Wind Nebulae, was worked out by J. Arons and collaborators in several steps (Alsop and Arons, 1988; Hoshino et al., 1992; Gallant and Arons, 1994) (see also Kennel and Coroniti, 1984). In this 1-D simulation, the incoming pair plasma is partially reflected back at a magnetic barrier, which is formed as the transverse magnetic field is compressed in the slow-down of the flow. A fraction of the incoming plasma is transmitted whereas a significant fraction is reflected making two opposite loops in phase space that end up in the downstream flow, producing a synchrotron resonant electromagnetic field that propagates both forwards and backwards in the form of an extraordinary mode. The backward mode, because of its compressive nature, heats the downstream flow.

When the flow is composed of protons and electrons, the intense coherent forward wave produces, through a ponderomotive force, an electrostatic wake of large amplitude, that brings the incoming electrons to a relativistic temperature up to equipartition with the protons. This eventually quenches the formation of the electrostatic wake by the ponderomotive force, which is proportional to $q^{2}/m$, once the relativistic mass of the electrons reaches that of the proton. Once electrons have been heated to near equipartition with the protons, the formation of the shock becomes similar to the case of a pair plasma.

This scenario depends on a crucial parameter called “magnetization,” which is defined as the ratio of the flux of magnetic energy across the shock to the flux of incoming matter energy:

[TABLE]

Note that $B_{\perp}$ is defined in the shock frame, following our notation convention. The magnetization parameter is sometimes defined relative to the incoming kinetic energy $\gamma_{\rm sh}\left(\gamma_{\rm sh}-1\right)n_{\rm u}mc^{2}$, but in the ultra-relativistic limit $\gamma_{\rm sh}\,\gg\,1$, these two definitions are equivalent.

Apart from the angular dependence, $\sigma$ depends essentially on the magnetization of the ambient plasma, since $B_{\perp}=\gamma_{\rm sh}B_{\perp|\rm u}$. In the typical interstellar medium, $\sigma\,\sim\,10^{-9}$, thus a ultra-relativistic shock in the interstellar medium corresponds to a very low magnetization.

3.2 Weakly magnetized shocks

In this subsection, we indicate the requirement for producing micro-turbulence in a weakly magnetized plasma. At very low magnetization, a self-consistent process develops such that a fraction of incoming particles is reflected back into the upstream region. These reflected particles trigger some form of streaming instability which generates precursor electromagnetic turbulence which eventually produces enough scattering to reflect particles back to the subshock.

Recall indeed from the previous section that particles that outrun the shock form a beam with opening angle $\,\simeq\,1/\gamma_{\rm sh}$ in the rest frame of the background plasma. Such a configuration is prone to beam-plasma instabilities.

The shock thus behaves as a self-sustaining self-generating structure. Moreover this process can allow for the development of a superthermal tail, which expands the size of the precursor (which increases as some power of gyroradius $r_{\rm g}$, hence of particle energy). In this description, superthermal particles originate from shock-reflected particles; since both populations share a common mean energy per particle of $\gamma_{\rm sh}$, we will not distinguish one from the other in the following and we will use “returning particles” or “superthermal particles” interchangeably.

Consider a streaming instability, triggered by the beam of returning particles, with a growth rate $g_{\rm inst|u}$ defined in the rest frame of the background plasma. In the present situation, an instability is relevant only if its growth time is shorter than the time spent by the background plasma in the shock precursor. Assume also that, although weak, the background magnetic field nevertheless sets the size of the precursor: $\ell_{\rm prec|u}\,\sim\,\omega_{\rm L,0}^{-1}/\gamma_{\rm sh}^{3}$. Therefore a streaming instability emerges as a good candidate for generating turbulence in the precursor if its growth rate is such that

[TABLE]

where $\omega_{\rm c}=eB_{\rm|u}/mc$ is the cyclotron frequency of the background plasma; the above uses $\omega_{\rm L,0}\,=\,\omega_{\rm c}/\gamma_{\rm sh}^{2}$ since $p_{\rm|u}\,\simeq\,\gamma_{\rm sh}^{2}mc$ for a superthermal particle in the background plasma rest frame. This is a severe constraint at a relativistic shock; expressed for a proton plasma as a function of magnetization, this constraint becomes:

[TABLE]

with $\omega_{\rm pi}\,\equiv\,\left(4\pi n_{\rm u}e^{2}/m_{p}\right)^{1/2}$. The process not only requires a very low magnetization, it also requires instabilities that work on micro-physical scales since the criterium (12) excludes slower MHD dynamics. This latter criterium is the essential reason for which the phenomenology of ultra-relativistic shocks displaying intense magnetic fields and superthermal tails, like the termination shocks of gamma-ray bursts and pulsar wind nebulae, unavoidably requires a micro-physical description.

3.3 Classification in term of the magnetization. Phase diagram

In this subsection, we present the micro-instabilities that are possible candidates for the generation of the magnetic turbulence. The maximum growth rates of standard beam-plasma instabilities are (within a factor of order unity) $(\omega_{\rm pb}^{2}\omega_{\rm ep}/\gamma_{\rm b|u}^{2})^{1/3}$ for the electrostatic two stream instability ($\gamma_{\rm b|u}$ represents the typical Lorentz factor of the beam of superthermal particles in the upstream plasma frame, see below), $(\omega_{\rm pb}^{2}\omega_{\rm ep})^{1/3}$ for the oblique two-stream instability, and $\omega_{\rm pb}$ for the Weibel instability, where $\omega_{\rm pb}$ is the (relativistic) plasma frequency of the beam of returning particles:

[TABLE]

and $\omega_{\rm ep}$ represents the electron plasma frequency of the background plasma. Using Eq. (9) with $n_{\rm cr|u}\,=\,\gamma_{\rm sh}P_{\rm cr}/\left(\gamma_{\rm sh}mc^{2}\right)$, and $\gamma_{\rm b|u}\,\simeq\,\gamma_{\rm sh}^{2}$, one finds $\omega_{\rm pb}\,\simeq\,\sqrt{\xi_{\rm cr}}\omega_{\rm p}$, where the conversion parameter $\xi_{\rm cr}$ is defined by relation (9) and $\omega_{\rm p}$ is the electron plasma frequency $\omega_{\rm pe}$ in a pair plasma, or $\omega_{\rm pi}$ in a proton electron plasma.

Note that a growth time scale of the order of $\omega_{\rm pi}^{-1}$ can be compatible with condition (12) if the magnetization is low enough, but not the proton cyclotron time $\omega_{\rm ci}^{-1}$. A detailed analysis reveals that the only magnetized waves which could be excited are the electron whistler waves (Lemoine and Pelletier, 2010).

The two stream instability and the oblique two stream instability grow through a Čerenkov resonance of the electrostatic (Langmuir) wave with the background electrons. When the background electrons are heated to relativistic temperature, the electrostatic waves become superluminal and thus these instabilities are quenched by resonance suppression (Bret, 2009; Lemoine and Pelletier, 2011; Shaisultanov et al., 2012). Nevertheless, the oblique two stream instability certainly plays an important role at the beginning of the precursor by heating the background electrons.

A phase diagram involving the two parameters $\sigma$ and $\gamma_{\rm sh}$ can be drawn to delimit the region in which micro-turbulence can be self-generated. For the most often discussed instability, namely the filamentation/Weibel instability, the criterion is simple: $\gamma_{\rm sh}^{2}\sigma<\xi_{\rm cr}$. This limit is a first order approximation to the development of Weibel-like turbulence in the precursor of relativistic shocks; in the following, we show that a refined description of the physics of the precursor introduces a “buffer effect” which modifies somewhat this limit at moderate magnetization.

3.4 Fermi acceleration

Similarly to what has been proposed for non-relativistic shocks, relativistic shocks can develop a superthermal tail thanks to the Fermi process. Indeed the excitation of magnetic turbulence that extends in the precursor and is also transmitted to the downstream flow allows for fast scattering of superthermal particles that can cross the shock front back and forth and thus gain energy. However at relativistic shocks, the Fermi process has some specific features.

We assume that the magnetic disturbances that scatter particles are almost frozen in the plasma flow both upstream and downstream (this will turn out to be true with Weibel turbulence). The particle energy gain through the Fermi process depends on the relative velocity $\beta_{\rm rel}$ between the upstream and downstream flows. Let a particle, having a pitch angle cosine $\mu_{1}$ before scattering, cross the shock front, be scattered and come back with a pitch angle cosine $\mu_{2}$. Its energy gain is

[TABLE]

As defined above, $\mu_{1}$ and $\mu_{2}$ are expressed in the upstream rest frame and $\beta_{\rm rel}$ represents the velocity of the downstream plasma relative to the upstream frame; changing the sign convention for $\beta_{\rm rel}$, one would obtain the same formula with $\mu_{1}$ and $\mu_{2}$ now both expressed in the downstream rest frame. In equation (15), and for the following discussion, it is implicitly assumed that the shock is unmodified by Fermi acceleration with a discontinuous jump between the upstream and downstream regions.

Choosing the downstream reference frame, $\beta_{\rm rel}\simeq 1-1/\gamma_{\rm sh}^{2}$ and the pitch angle cosine of a particle crossing the front towards downstream has $-1<\mu_{1}<\beta_{\rm sh}$; a particle coming back from downstream to upstream has $\beta_{\rm sh}<\mu_{2}<1$.

The average energy gain can be calculated easily in the case of non-relativistic shocks only, because the distribution functions are almost isotropic and $\mu_{1}$ and $\mu_{2}$ independent; one then finds:

[TABLE]

Moreover, because of the almost isotropic distributions, the probability of escape through advection in the downstream plasma can also be easily estimated: $P_{\rm esc}=4\left|\beta_{\rm d}\right|$. The index of the power law distribution ${\rm d}N/{\rm d}p\,\propto\,p^{-s}$ is then obtained as Bell (1978):

[TABLE]

i.e. $s=2$ for the strong adiabatic shock where $\left|\beta_{\rm rel}\right|=3\left|\beta_{\rm d}\right|$.

For relativistic shocks, the strong anisotropy of the distribution functions makes these calculations more complicated, (e.g. Ellison et al., 1990a; Gallant and Achterberg, 1999; Kirk et al., 2000; Achterberg et al., 2001) see however Vietri (2003), Lemoine and Pelletier (2003), or Keshet and Waxman (2005) for an alternative point of view. If one were to consider a random incident pitch angle cosine $\mu_{1}$, then $\mu_{2}\simeq 1-{\cal O}(1/\gamma_{\rm sh}^{2})$ would imply $\langle G\rangle\sim\gamma_{\rm sh}^{2}$. However, this maximum limit can only be achieved for particles crossing from downstream to upstream for the first time. Indeed, particles participating in a further Fermi cycle are caught up by the front upstream with $\mu_{1}\simeq 1-{\cal O}(1/\gamma_{\rm sh}^{2})$ so that the gain falls to a few only; precisely it is always close to 2 after the first Fermi cycle (e.g. Achterberg et al., 2001; Lemoine and Pelletier, 2003).

The index can be defined through the equation

[TABLE]

and because $G$ remains close to 2 at each cycle, the index turns out to be close to

[TABLE]

The numerical estimate is $1-P_{\rm esc}\,\simeq\,0.4$ leading to $s\simeq 2.3$. An illustration of the spectrum of accelerated particles for an unmagnetized, unmodified shock with $\gamma_{\rm sh}\,=\,100$ is shown in Fig. 3.

As we describe in Sec. 6, basic requirements of energy and momentum conservation mean the test-particle power law shown in Fig. 3 will be modified to some extent if shocks accelerate cosmic rays efficiently. If the scattering mean free path is an increasing function of particle momentum as expected, the accelerated spectrum will obtain a concave shape being harder at high energies before a turnover produced by finite shock effects.

In addition, if the shock is sweeping up non-relativistic matter, such as the ISM, the injection of highly anisotropic thermal particles must be described consistently with CR production. Currently, the only techniques that can handle thermal particle injection and efficient non-linear CR production in relativistic shocks are PIC simulations and Monte Carlo techniques, described in greater detail in Sec. 6.

4 Physics of the precursor

The main physical processes that develop in the precursor of a relativistic collisionless shock undergoing Fermi acceleration are the following:

(i) The incoming precursor plasma slows before the sharp transition at the viscous subshock due to transverse currents. This occurs with ordered or disordered magnetic fields for both pair and ionic plasmas and produces an important buffer effect in the precursor where a highly relativistic shock far upstream turns into a mildly relativistic one deep in the precursor, close to the shock front. If Fermi acceleration is efficient the backpressure of accelerated particles produces an additional shock modification forced by conservation of energy and momentum.

(ii) We argue that the Weibel/filamentation instability is the dominant generator of precursor magnetic turbulence; this instability does not quench when the electrons are heated to relativistic temperatures, but it saturates when the beam dispersion becomes too large.

(iii) Then we develop the important concept of the Weibel frame for both pair and ionic plasmas. This is the frame where the turbulence is purely magnetic and where particles undergo a simple pitch-angle scattering process.

(iv) We discuss the crucial heating process of background electrons and relate it to the drag force experience by the background particles.

(v) We discuss some aspects of the non-linear theory for the formation of the filaments.

(vi) We then return to the physics of particle scattering in Weibel turbulence and derive it from the the theory of Hamiltonian chaos (the details are given in App. A.1).

4.1 Dynamics in the precursor with a mean field; buffer effect

The slow-down of the incoming plasma with respect to the shock front will be analyzed in a dedicated Section 4.5 for the case of an unmagnetized plasma. In this section we examine the simpler case of the slow-down due to a transverse mean field.

4.1.1 The case of a pair plasma

Consider first the beam of returning particles. Because of the transverse mean magnetic field (aligned say along $\boldsymbol{z}$), electrons and positrons are deflected in opposite directions (along $\boldsymbol{y}$) during their orbit around $B$ and thus a current transverse to the shock normal and to the mean field is generated, with density $j_{\rm b}=n_{\rm cr}ec$ (measured here in the front frame). One can show that this current is compensated to an excellent accuracy and on a short timescale by a current carried by the incoming background plasma, once it enters the precursor. As the density of this background current density can be written $j_{\rm bg}=\gamma_{\rm sh}n_{\rm u}e\beta_{y}c$ in the front frame, with $n_{\rm cr}\,=\,\xi_{\rm cr}\gamma_{\rm sh}n_{\rm u}$, one infers a transverse velocity $\beta_{y}\,\simeq\,\xi_{\rm cr}$ (Lemoine et al., 2014b, a).

The main slow-down effect on the incoming plasma is due to the deviation of background electrons and positrons due to the development of the compensating current. Indeed, one can show that the total Lorentz factor of the background plasma is conserved, while the development of a transverse velocity $\pm\xi_{\rm cr}$ for electrons/positrons implies that the center-of-mass of the background plasma moves solely along $\boldsymbol{x}$ with a Lorentz factor:

[TABLE]

Thus any shock front of Lorentz factor larger than $1/\xi_{\rm cr}$ produces a deceleration of the incoming plasma such that the center of mass moves with a Lorentz factor $1/\xi_{\rm cr}$. This is a crucial buffer effect that transforms any ultra-relativistic shock into a mildly one in its precursor from the point of view of the growth of instabilities.

As discussed in Lemoine et al. (2014b, a), this modification of the motion of the center-of-mass indeed modifies the criterion for the growth of instabilities in the precursor of a shock, due to time dilation effects. One can show in particular that the Weibel instability can now grow provided $\sigma<\xi_{\rm cr}^{3}$, independently of the shock Lorentz factor, as illustrated in Fig. 2. Moreover, the above configuration of a compensated transverse current turns out to be unstable on a short timescale, leading to the growth of a current-driven filamentation instability with growth rate $\sim\,\omega_{\rm pe}$, i.e. faster than the Weibel/filamentation rate. This latter instability can thus grow provided $\sigma<\xi_{\rm cr}^{2}$ at any value of the shock Lorentz factor, thereby bridging the gap between weakly and moderately/strongly magnetized shock waves.

4.1.2 The case of an electron-ion shock

We consider first an electron-ion plasma where electrons have been rapidly heated to a relativistic temperature of order $m_{p}c^{2}$ (this heating will be explored in a next section). Their inertia is thus comparable to that of the protons, so that a significant fraction of the electrons is reflected back on the shock front along with protons, leading to a charge neutral beam of returning particles. In such a case, the situation is similar to the case of a pair plasma discussed above. If the electrons are relativistically hot, but below equipartition with the protons, i.e. $T_{e}=\xi_{\rm th}m_{p}c^{2}$ with $\xi_{\rm th}>m_{e}/m_{p}$, one can recover similar results, as will be presented elsewhere.

If the population of incoming electrons remains cold because the precursor is too short, then electrons are not reflected back and the beam of returning current carries a net electric charge $\rho_{\rm cr}=\xi_{\rm cr}\gamma_{\rm sh}n_{\rm u}e$ (front frame). As the background plasma compensates this electric charge, a force density develops along the y-direction, transverse to the shock normal and to the mean field $B$; in the relativistic MHD approximation, this force density reads:

[TABLE]

The above equation uses a zero current density along $x$, which is guaranteed by the assumption of one-dimensional stationary dynamics in the shock front frame. The work of this force over the length of a Larmor gyration is $\xi_{\rm cr}\beta_{\rm sh}\gamma_{\rm sh}^{2}n_{\rm u}m_{p}c^{2}$. Thus protons are deviated in the transverse direction with a velocity $\xi_{\rm cr}$ and the bulk motion is diminished as previously.

4.2 Weibel turbulence

Because the mean magnetic field must be very weak and the growth of the instability faster than the cyclotron frequency of electrons in a pair plasma, or of protons in a protonic plasma, the expected turbulence must be excited by a fast micro-instability on short length scales. The typical scale must be at most the inertial scale of protons: $\delta_{i}\equiv c/\omega_{\rm pi}$; MHD develops at larger scales than this inertial scale.

The Weibel instability, which oversteps the other streaming instabilities and survives to electron heating, is mostly an electromagnetic instability with a very low phase velocity, and thus produces slowly propagating magnetic turbulence, quite suitable for particle scattering and thus for a Fermi process. The mechanism of the instability is simple: the beam of returning particles is supposed to be charge neutralized, although when a transverse magnetic wave has been produced, its quiver force separates electric charges; this generates an intense electric current, which in turn amplifies the primary magnetic wave; then it blows up (see e.g. Medvedev and Loeb, 1999; Wiersma and Achterberg, 2004; Lyubarsky and Eichler, 2006; Milosavljević and Nakar, 2006; Achterberg and Wiersma, 2007; Achterberg et al., 2007). Its growth rate is fast, of the order of $\omega_{\rm pb}$, with respect to the background frame, and the instability is purely growing in this frame with a weak electrostatic component. A sketch of the mechanism of this instability in the precursor of a relativistic shock is presented in Fig. 4.

We already mentioned that relativistically hot electrons do not quench this Weibel instability (to be discussed shortly); however a large angular dispersion of the beam can stabilize the Weibel instability, as discussed e.g. in Rabinak et al. (2011); Lemoine and Pelletier (2011). A simple description of that saturation effect by the beam dispersion consists in treating the transverse dispersion of the beam as a temperature effect; this introduces an effective sound velocity $\beta_{*}c\sim c/\gamma_{\rm sh}$, see Eq. (69), which modifies the hydrodynamic response of the beam fluid to the quiver force that separates opposite charges. The instability then disappears (in this fluid description) when $\Im\omega\,\lesssim\,\beta_{*}\omega_{\rm p}$, i.e. when $\gamma_{\rm sh}$ becomes smaller than $1/\xi_{\rm cr}^{1/2}$. A more detailed analysis with some consequences of this saturation effect will be addressed in a forthcoming paper.

4.3 The Weibel frame in a pair plasma

For the estimation of the efficiency of scattering and thus of the Fermi process, and also for estimating the radiation efficiency of an ultra-relativistic shock, it is very important to know what is the proper frame of Weibel turbulence.

The Weibel frame is defined as the frame where the electrostatic potential vanishes. Its determination is detailed in App. A.1, where it is shown that it propagates at the speed $\xi_{\rm cr}$ with respect to the background plasma. This has two important consequences, namely that the high energy returning particles are scattered with respect to that frame and that the background plasma suffers a drag force which slow it down and heats the electrons. The Weibel frame is the suitable frame to calculate the growth rate of the instability, which is purely magnetic and purely growing, since the charge separation in the beam is exactly compensated by an opposite charge separation in the background plasma.

In the front frame, a strong electrostatic potential $\delta\Phi$ is set-up; its magnitude can be determined by the Lorentz transform from the Weibel frame to the front frame: $\delta\Phi\,=\,\gamma_{\rm w}\beta_{\rm w}\delta A_{x|\rm w}$, defining $\beta_{\rm w}$ (resp. $\gamma_{\rm w}$) as the velocity (resp. Lorentz factor) of the frame in which $\delta\Phi$ vanishes, relatively to the shock front frame. Similarly, the $x-$component of the magnetic potential reads $\delta A_{x}\,=\,\gamma_{\rm w}\delta A_{x|\rm w}$, so that $\delta\Phi\,=\,\beta_{\rm w}\delta A_{x}\,\simeq\,-\beta_{\rm sh}\delta A_{x}$. In the front frame, there must be a net total electric charge density fluctuation $\delta\rho$, also determined through a Lorentz transform: $\delta\rho=\gamma_{\rm w}\beta_{\rm w}\delta j_{x|\rm w}/c$; the net current density fluctuation is $\delta j_{x}=\gamma_{\rm w}\delta j_{x\mid\rm w}$ so that $\delta\rho\,\simeq\,-\beta_{\rm sh}\delta j_{x}/c$. It is assumed here that in the Weibel frame, the total charge density perturbation $\delta\rho_{|\rm w}$ vanishes. The Maxwell equation that determines the electromagnetic field as a function of the current density, namely $\Delta\delta A_{x}=-4\pi j_{x}/c$, is then the same as the Poisson relation between the electrostatic potential and the electric charge.

To conclude this section on the Weibel instability in a pair shock, let us show that it is quite convenient to calculate the growth rate of the Weibel instability in the Weibel frame. Indeed in this frame $\delta\Phi_{\rm|w}=0$ and $\delta\rho_{\rm tot|w}=0$, which allows one to write two simple coupled equations to derive the instability growth rate in the cold approximation. Switching over to Fourier space, assuming $k_{x}\,=\,0$ as before, we insert $\delta\rho_{\rm b|w}=-\delta\rho_{\rm bg|w}$ and $\beta_{\rm bg|w}=-\xi_{\rm cr}$ in the dynamical equation $\partial_{\alpha}\partial^{\alpha}\delta A^{\mu}\,=\,-4\pi\delta j^{\mu}$ and in the background dynamical equation, to obtain:

[TABLE]

The plasma frequency of the background plasma emerges in the first equation because of the response of the conduction current to the electromagnetic perturbation. In details:

[TABLE]

and the conduction current density perturbation $\rho_{\alpha}\delta\beta_{x\alpha}$ is obtained through the dynamical equation for $\delta\beta_{x\alpha}$:

[TABLE]

The strong dependence on $\gamma_{\alpha}$ justifies the neglect of the conduction response of the beam of superthermal particles; furthermore, in the Weibel frame, $\gamma_{\rm bg|w}\,\simeq\,1$.

One thus derives the dispersion relation to leading order in $\omega/kc$:

[TABLE]

This is a purely growing magnetic mode. If we take into account the case $k_{x}\neq 0$, two modifications are introduced: a complete Laplacian operator in the Maxwell equation and a Doppler effect in the dynamical equation of the background plasma, with $\partial_{t}^{2}\mapsto(\partial_{t}-\xi_{\rm cr}c\partial_{x})^{2}$ which leads to a simple modification of the dispersion equation:

[TABLE]

The growth rate is invariant by parity; thus in the Weibel frame, forward and backward waves propagate with a frequency $\xi_{\rm cr}k_{x}c$ with the same growth rate. A strong Doppler effect is produced in the front frame, so that both forward and backward waves are caught up by the shock front.

4.4 The Weibel instability in an electron-ion shock

We now assume a plasma made of protons and hot electrons and extend the calculations to include the Landau effects in the electromagnetic response. The detailed calculations are presented in App. A.2.

As for the pair shock, one can derive the velocity of the Weibel frame by analyzing the generation of electric charge by the quiver force, taking into account the susceptibilities of the various plasma components. To circumvent the technical difficulty associated with the Lorentz transform of these susceptibilities from one frame to another, it is more convenient to rely on a Lorentz covariant formulation of linear response theory, as formulated by Melrose (1986) and Achterberg and Wiersma (2007), which involves the polarization tensor of species $s$ $\alpha^{\mu}_{s,\nu}$ defined through $\delta j_{s}^{\mu}=\alpha^{\mu}_{s,\nu}\delta A^{\nu}$. The combination of Maxwell equations in the Lorentz gauge, namely $\partial_{\mu}\partial^{\mu}\delta A^{\nu}\,=\,-\sum_{s}\delta j_{s}^{\nu}$ and of the above response then leads to the definition of the Weibel frame as that in which

[TABLE]

The calculations of $\alpha^{01}_{s}$ in a generic frame, for streaming hot electrons, are too lengthy to be reported here; these calculations will be presented elsewhere. Nevertheless, one obtains the same result as before, $\beta_{\rm w|u}\,\sim\,\xi_{\rm cr}$.

The concept of Weibel frame is essential in the prospect of particle acceleration and radiation, since this is the frame in which particles are purely scattered (at least almost, for there is an induction electric field significantly smaller than the magnetic field). The conclusion of this section is that particles are scattered by the magnetic field of Weibel turbulence with respect to a frame moving forward at sub-relativistic speed ($\sim\,\xi_{\rm cr}$) with respect to the background plasma.

4.5 Electron heating and the drag force

Electrons of the background plasma can be heated through two possible electric fields generated by the Weibel instability: one is the strong motional electrostatic field measured in the front frame due to the Lorenz transform, oriented both perpendicular to the magnetic field fluctuation and to the shock normal; the other is the induction electric field, which derives from the growth of the magnetic field energy density, and which is significant during the stage of instability growth.

The electrostatic field $\delta E_{\perp}$ is associated with a potential fluctuation which has been obtained in the front frame through a Lorentz transform: $\delta\Phi\,=\,\beta_{\rm w}\delta A_{x}$. The heating due to this field is caused by the scattering of electrons in the Weibel frame, which leads to a kind of second order Fermi process in the front frame, because a variation of $x-$momentum in the Weibel frame produces a variation of energy in the front frame. This is the strongest electric field when measured in the front frame with $\delta E_{\perp}\sim\delta B$.

In the Weibel frame, the inductive electric field derives from:

[TABLE]

Therefore $\delta E_{x}^{2}\,\sim\,\xi_{\rm cr}\delta E_{\perp}^{2}$ in the front frame and $\xi_{\rm cr}\,<\,1$; this inductive electric field nonetheless provides an efficient source of heating (Gedalin et al., 2012).

4.5.1 Scattering and heating in the electrostatic field

In a first approximation, background electrons are scattered by a static magnetic field in the Weibel frame that moves at speed $\beta_{\rm w|u}\,\simeq\,\xi_{\rm cr}$ with respect to the background frame. Therefore, in the background frame, their temperature increases and they suffer a drag force. Because their energy does not change during scattering in the Weibel frame (as long as one neglects the induction field effect), they suffer a modification of their $x-$momentum $\Delta p_{x|\rm u}$ and energy $\Delta\epsilon_{|\rm u}$ in the background frame.

Together with $\Delta\epsilon_{|\rm w}\,=\,0$, the Lorentz transform $\Delta\epsilon_{|\rm w}=\gamma_{\rm w|u}(\Delta\epsilon_{|\rm u}-\beta_{w}c\Delta p_{x|\rm u})$ provides an essential relation between energy gain and momentum variation:

[TABLE]

The stochastic force $\Delta p_{x|\rm u}/\Delta t_{|\rm u}$, after statistical averaging, leads to a macroscopic drag force $F_{\rm d|u}$ and the distribution of background electrons, isotropic in its rest-frame, gets heated at a rate (background frame):

[TABLE]

The drag force is best evaluated in the Weibel frame through a Fokker-Planck term:

[TABLE]

with a momentum diffusion coefficient of the form

[TABLE]

with $h_{ij}^{\perp}\,=\,\delta_{ij}-p_{i|\rm w}p_{j|\rm w}/p_{\rm|w}^{2}$ the projector orthogonal to momentum. The timescale $\tau_{\rm c|w}$ represents the coherence timescale of the random force associated to the Weibel turbulence. Using the Lorentz invariance of $f$ and the fact that $f$ depends only on $\epsilon_{\rm|u}/T$ in the background frame, one can derive

[TABLE]

up to a prefactor of order unity.

Using the above, we can estimate the temperature that can be achieved and the slow-down of the background flow that ensues. As discussed further below – see Eq. 45 in particular– the precursor length is approximately $\delta_{i}/\xi_{B}$, as measured in the background rest-frame, and therefore the background plasma experiences the drag force for a duration $\xi_{B}^{-1}\omega_{pi}^{-1}$. The temperature achieved is then roughly of order

[TABLE]

In the background frame, the dragging effect is seen as a force that pushes the plasma forward, with an intensity $F_{\rm d|u}=F_{\rm d|w}$. Over the precursor crossing timescale, the upstream plasma is displaced in energy by

[TABLE]

indicating that in the front frame, the variation in momentum is of order unity if $\tau_{\rm c|w}\omega_{\rm pi}\,\sim\,1$.

Therefore this scattering process acting on the background plasma is able to slow-down the whole incoming plasma, with possible consequences regarding the saturation of the Weibel instability in regards of the angular dispersion of the beam.

A more complete calculation should take into account the spatial profiles of the various quantities in a dynamical system; however the former estimates provide useful orders of magnitude regarding heating and slow down.

4.5.2 Heating in the inductive electric field

The induction field acts during the phase of instability growth and its heating effect can be estimated with quasi-linear theory which leads to (in the Weibel frame)

[TABLE]

with $\Gamma_{xx}\,\simeq\,(g_{\rm inst|w}\omega_{\rm pi})^{2}e^{2}<\delta B_{\perp}^{2}>\tau_{\rm c|w}$. The heating rate is thus a factor $\omega_{\rm pi}^{-1}g_{\rm inst|w}/\xi_{\rm cr}$ times that in the stochastic electrostatic field; in the cold plasma limit, both are comparable because $g_{\rm inst|w}\,\sim\,\xi_{\rm cr}\omega_{\rm pi}$, but at larger temperatures, heating in the inductive field is slightly larger. However, this heating process in the inductive field only applies to the linear growth phase; it is much reduced in the non-linear phase. Overall, one thus expects to reach a scaling given by Eq. (34).

One can also estimate the electron temperature from the phenomenology of the relativistic regime of oscillations in very intense fields: $T_{e}\sim eB\ell_{\rm c}$. The typical transverse scale of turbulence is the electron Debye length: $\ell_{\rm c}=\lambda_{\rm D}=\xi_{\rm th}^{1/2}\delta_{i}$ in the hot electron regime. Then we get another estimate compatible with the previous one:

[TABLE]

This is a remarkable result that suggests that the thermal energy density follows the magnetic energy density and that confinement structures can occur in the background plasma.

4.5.3 Weibel filaments

The above relation suggests that heating of the background plasma ensures that the thermal pressure inside the filaments counterbalance the magnetic stress. It is important here to note that the filaments are built out of the background plasma, not by filamentation of the beam of returning particles, as is often thought. Indeed, the transverse momentum of these particles in the background plasma frame is of the order of $\gamma_{\rm sh}mc$, which is too large to allow them to be confined in structures of transverse size $\delta_{i}$, since $e\delta B_{|\rm u}\delta_{i}\,\sim\,\xi_{B}^{1/2}mc^{2}$.

The filaments observed in PIC simulations are thus real, but they form in the background plasma. This is tightly related to the fact that the background plasma moves at a velocity $-\beta_{\rm w|u}$ with respect to the magnetic field. These filaments are some sort of coherent structures of a quasi bi-dimensional turbulence, where inverse cascade develops.

One can derive the general structure of these filaments in the transverse plane by considering the inhomogeneous equilibrium of the background plasma in the presence of a finite $A_{x}$ component (background plasma frame). Using a non-relativistic description, appropriate to the transverse dynamics, one finds from the force balance in the transverse direction:

[TABLE]

The possibly residual kinetic energy $E_{K}=mv^{2}/2$ is in the form of irrotational motions, the Lorenz driving being irrotational. One obtains Boltzmann type equilibria, assuming the same temperature for both populations of particles:

[TABLE]

In a pair plasma, one gets a charge separation such that

[TABLE]

Because in motion of velocity $-\xi_{\rm cr}$ with respect to the Weibel frame, its associated current generates a magnetic field such that:

[TABLE]

Defining $a\,\equiv\,\xi_{\rm cr}eA_{x}/T$, we obtain a sinh-Poisson equation, whose analytical solutions are known:

[TABLE]

The solutions of that equation are known to be in the form of a double lattice of filaments of alternating polarity. The size of the filaments appears clearly as $\lambda=\lambda_{\rm D}/\xi_{\rm cr}=\sqrt{\xi_{\rm th}}/\xi_{\rm cr}\delta_{i}$. This estimate does not account for the kinetic energy contribution that increases the Debye length. However its contribution is expected to remain not larger than $\xi_{\rm cr}eB\lambda\,\sim\,|a|T$, which is also the estimate of the achieved temperature. In a pair plasma that energy is sub-relativistic.

In the case of an electron-ion plasma, the solution are similar if the kinetic energy is negligible. Otherwise, it breaks the symmetry by giving more importance to the electron contribution. In that case instead of a sinh-Poisson equation, the filaments are described by a Poincaré equation, whose solutions are of Bennett equilibrium type with all the filaments having the same polarity and thus undergoing more merging than in the sinh-Poisson case:

[TABLE]

4.6 Scattering off Weibel turbulence

The size of the precursor is determined by the residence time of returning particles. As mentioned earlier, and as seen from the background frame, they are outrun by the shock front when their pitch angle has opened up to $1/\gamma_{s}$ due to a scattering at a frequency $\nu_{s}$ such that:

[TABLE]

The Lorentz factor of beam particles in the background plasma rest frame $\gamma_{\rm b|u}\,\simeq\,\gamma_{\rm sh}^{2}$. The above estimate derives from a random phase approximation (Milosavljević and Nakar, 2006; Pelletier et al., 2009) and turns out to be valid even with a more elaborated method based on Hamiltonian chaos, described in App. A.3.

Measured from the background frame one obtains:

[TABLE]

assuming that $\tau_{\rm c|u}\,\sim\,\omega_{\rm pi}^{-1}$. Thus the travel time of the background plasma across the precursor is about $(\xi_{B}\omega_{\rm pi})^{-1}$, which is longer than the growth time of the instability $(\sqrt{\xi_{\rm cr}}\omega_{\rm pi})^{-1}$, while the residence time of superthermal particles of Lorentz factor $\gamma_{\rm|u}$ in the precursor is of the order of

[TABLE]

This residence directly controls the acceleration timescale, which thus scales as the energy squared of the superthermal particle. This result has been confirmed by PIC simulations (Sironi et al., 2013).

The value of the coherence time for the superthermal particles is not a trivial issue, as mentioned in Achterberg et al. (2007). Indeed, in an idealized Weibel turbulence with $k_{x}\,=\,0$, the conjugate momentum $p_{x}-eA_{x}/c$ is a conserved quantity of the particle trajectories. Since $p_{x}\,\sim\,\gamma_{\rm sh}mc$ while $eA_{x}/c\,\sim\,\xi_{B}^{1/2}\gamma_{\rm sh}mc$ in the front frame, the former largely dominates, hence superthermal particles would never be able to go back to the shock if the turbulence were exactly $x-$independent. Consequently, the finite $k_{x}$ must play a key role in the transport of superthermal particles.

The general scaling of the residence time with energy $\epsilon$, hence of the acceleration timescale $t_{\rm acc}\,\propto\,\epsilon^{2}$, has been obtained in the framework of a random phase approximation model in Pelletier et al. (2009). Nevertheless, one can also derive this scaling using a model of Hamiltonian chaos applied to the trajectories of superthermal particles. This leads to a validation of the scattering law in a more general situation than the random phase approximation. In App. A.3, we present a simple model of such chaotic dynamics; a more complete model will be presented in a forthcoming publication.

5 Phenomenological consequences

5.1 Scattering efficiency and maximal energy

5.1.1 Onset of acceleration

As discussed in Sec. 2.2, an oblique background magnetic field in the upstream becomes essentially transverse in the downstream flow, leading to a super-luminal shock configuration. As discussed in Lemoine et al. (2006) and Pelletier et al. (2009), such a configuration inhibits Fermi acceleration unless the scattering frequency in the turbulence exceeds the gyrofrequency in the background field (which controls the rate at which particles are advected in the downstream plasma). In Lemoine et al. (2006), it was actually speculated that the development of Weibel micro-turbulence should lead to an efficient Fermi process and this point of view has been confirmed by PIC simulations.

For scattering in micro-turbulence in the downstream flow, the above condition for the development of Fermi cycles can be written:

[TABLE]

with $\tau_{\rm c\mid d}\sim\omega_{\rm pi}^{-1}$. Here $\gamma$ denotes the Lorentz factor of the accelerated particle in the downstream frame; for the first Fermi cycle, $\gamma\,\sim\,\gamma_{\rm sh}$. The above can be written in terms of the parameters $\sigma$ and $\xi_{B}$ as

[TABLE]

It thus indicates that only weakly magnetized shock waves should be able to accelerate particles efficiently. For typical values $\xi_{B}\,\sim\,10^{-2}$ in the vicinity of the shock front, the above bound suggests that particle acceleration should take place in shocks with $\sigma\,\lesssim\,10^{-4}$. This point of view, developed in Lemoine et al. (2006) and Pelletier et al. (2009) (see also Niemiec et al. (2006) for numerical simulations) has been beautifully confirmed by PIC simulations, e.g. Sironi et al. (2013).

5.2 Maximal energies

The maximal energy that can be reached by the Fermi process in a shock of a given size and age depends on several factors, in particular the scattering rate [Eq. (44)], the energy loss or escape timescale, and the background magnetic field. The above discussion indicates that in the front frame, the upstream residence time typically dominates over the downstream residence time; this effect arises because the upstream fluctuations move at high velocity relative to the front, while the downstream perturbations are essentially static and isotropic in the shock front frame. However, as discussed in detail by Plotnikov et al. (2013), the inclusion of a small but finite background magnetic field changes the scale of the precursor.

To simplify the discussion, one can use the downstream residence time as a proxy for the acceleration timescale; one would obtain similar results for a more careful calculation taking into account the effect of a background magnetic field comparable to that of the interstellar medium and $\gamma_{\rm sh}\,\sim\,100$, representative of gamma-ray bursts external shocks.

One can then show that the maximal proton energy is of order $10^{16}\,$eV at the external shock of a gamma-ray bursts. In contrast, the maximum electron energy is limited by synchrotron energy losses to values of the order of $10\,$TeV, independent of the magnetic field and shock Lorentz factor (Kirk and Reville, 2010; Plotnikov et al., 2013; Wang et al., 2013). This can give rise to GeV synchrotron photons at early times $\sim\,10^{2}-10^{3}\,$s, as observed by the Fermi satellite in some gamma-ray bursts.

5.3 Radiative diagnostics

Weakly magnetized shocks, which are prone to Fermi acceleration, are those in which intense micro-turbulence can be excited in the shock precursor by streaming instabilities. The fate of this turbulence downstream of the shock, on timescales $\gg\,\omega_{\rm pi}^{-1}$ has been the subject of intense debate. The linear dissipation rate of microturbulence is $\gamma_{k}\,\simeq\,k^{3}c^{3}/\omega_{\rm pi}^{2}$ in a relativistic plasma (e.g. Chang et al., 2008), a result which extends into the mildly non-linear regime (Lemoine, 2015a), suggesting that turbulence on scales $\sim\,\omega_{\rm pi}^{-1}$ is dissipated on $\omega_{\rm pi}^{-1}$ timescales. In contrast, the typical (comoving) size of the blast of a gamma-ray burst is of the order of $R/\gamma_{\rm sh}\,\sim\,10^{7}\,c/\omega_{\rm pi}$ (for an external density $n\,=\,1\,$cm*-3*), orders of magnitude larger, and this length scale gives the order of magnitude over which most electrons radiate their energy through synchrotron radiation. An important question therefore is whether Weibel turbulence, self-generated in the shock precursor, can survive over the above time/length scales in order to provide the turbulence in which synchrotron photons are produced, see e.g. Medvedev and Loeb (1999); Gruzinov and Waxman (1999); Chang et al. (2008); Keshet et al. (2009); Lemoine (2013, 2015a, 2015b) for discussions of this issue.

While additional instabilities or inverse cascade might potentially help sustain or amplify this turbulence, one may also conceive a simpler scenario, in which the power spectrum of magnetic fluctuations at the shock contains long wavelength modes, with a long lifetime $\propto k^{-3}$. As discussed in Lemoine (2013), these modes could be excited by the highest energy particles accelerated at the shock front, which propagate over large distances in the upstream. The gradual dissipation of the micro-turbulence by phase mixing would then lead to a $\xi_{B}$ which now depends on some power of the distance to the shock. The synchrotron spectrum produced by a population of particles interacting with such a dissipative turbulence bears some unique features, which could potentially be used as a diagnostic of Weibel turbulence (Lemoine, 2013, 2015a). A generic prediction, in particular, is a significantly larger ratio of flux at high energy to flux at low energy, because high energy photons are produced by high energy electrons, which radiate on a shorter timescale, hence in a stronger turbulence. Such a signature may have been detected in some gamma-ray bursts, whose multi-wavelength light curve can be nicely explained by a synchrotron spectrum resulting from a dissipative turbulence with $\xi_{B}\,\propto\,(\omega_{\rm pi}x/c)^{-0.5}$.

6 Macroscopic Simulations of Relativistic Shocks

The physics of relativistic collisionless shocks together with its application to the phenomenology of powerful astrophysical sources, are complex non-linear multi-scale problems, which have been addressed with a variety of methods (see e.g. Bykov and Treumann, 2011; Marcowith et al., 2016). Significant progress in our understanding of these phenomena has come through particle-in-cell (PIC) simulations, which probe the full non-linear relationship between the particles and the electromagnetic fields, see Sironi et al. (2015) and references therein. Such simulations are however computationally intensive and cannot yet explore the physics of shock waves on the long macroscopic time and space scales on which the source evolves, or on which radiation is produced. In this context, the analytical theory that we have proposed in the previous section is useful because it can be benchmarked against these PIC simulations on small physical scales and extrapolated to the macrophysical scales of interest. An alternative method is the numerical Monte Carlo technique which we describe in this Section.

Monte Carlo methods offer a modeling technique that lies between PIC simulations and semi-analytic models. In contrast to semi-analytic descriptions which make a global diffusion approximation, particles pitch-angle scatter locally with the Monte Carlo technique (e.g., Ellison and Eichler, 1984). The important difference is that for the diffusion approximation to be valid, particle distributions must be nearly isotropic in all frames and large flow gradients must be avoided. While these approximations are good for superthermal particles in non-relativistic shocks, they are never good for particles crossing relativistic shock front. In the Monte Carlo code, particles undergo a small pitch-angle scattering event after a time small compared to a gyroperiod. Thus the scattering is “local” and highly anisotropic distributions in sharp flow gradients can be treated directly. Of course, the properties of the scattering are not calculated from first principles, as done with PIC simulations. But scattering can be directly parameterized, as so far done for relativistic shocks (e.g., Ellison et al., 2013), or determined more fundamentally from CR-induced plasma instabilities, magnetic field amplification (MFA), and an analytic determination of scattering rates from the self-consistently determined magnetic turbulence power spectrum in non-relativistic shocks (e.g., Vladimirov et al., 2008; Bykov et al., 2014).

The analytical model described in the previous sections can be seen as an effective two-fluid theory, describing how self-generated electromagnetic turbulence affects and couples a fluid of superthermal particles with a background plasma fluid. As such, it can accurately model the back-reaction of the accelerated particles and of the background plasma on the turbulence. In contrast, the Monte Carlo model is a one-component theory. The pitch-angle scattering approximation treats all particles the same and Fermi accelerated particles arise smoothly from the background plasma which is a self-consistently described collection of particles rather than a separate fluid. The non-linear relationship between the Fermi accelerated population and the bulk plasma is determined by iteratively solving the energy-momentum tensor by forcing energy and momentum conservation across the shock.

Although the relativistic Monte Carlo simulations that follow do not make a distinction between the rest frame of the one-fluid plasma and the rest frame of the turbulence, in contrast to Sec. 4.2, they will be generalized in the near future to a more realistic setting in which the growth of turbulence and its rest frame are dictated by microphysical considerations. Monte Carlo simulations of non-relativistic shocks which incorporate these effects are described in the following as an illustration.

Most Monte Carlo models thus far have assumed plane, steady-state shocks with Fermi acceleration limited by free escape boundaries. Parallel (e.g., Ellison and Eichler, 1984; Ellison et al., 2016), oblique (e.g., Baring et al., 1993; Summerlin and Baring, 2012) and perpendicular (e.g., Takamoto and Kirk, 2015) shock geometries have been modeled.

The main features of the Monte Carlo technique are: (1) a large dynamic range can be modeled with reasonable computing resources following particles consistently from thermal energies to the highest CR energies with no constraint on the anisotropy of the particle distribution; (2) the non-linear (NL) feedback on the shock structure from efficient Fermi acceleration can be determined consistently with injection and the production of ultra-high-energy CRs; (3) since particle anisotropies are treated consistently, no distinction in the particle transport need be made between non-relativistic, trans-relativistic, or ultra-relativistic shocks; (4) electrons and ions can be simultaneously injected and accelerated allowing effects depending on particle mass to be determined self-consistently; and (5) absolutely normalized radiation, from electrons and ions, can be calculated as a straightforward function of model and environmental parameters.

The main drawbacks of this method are that it is intrinsically steady state, generally done for plain-parallel shocks, and thus far too computationally slow to include in hydrodynamic evolutionary models such as those for young SNRs using semi-analytic descriptions of NL diffusive shock acceleration (e.g., Ellison et al., 2007; Lee et al., 2012; Ellison et al., 2012; Zirakashvili and Ptuskin, 2012).

We describe here some physical results obtained with the Monte Carlo technique which may help to understand the backreaction effects of the accelerated particles on the shock flow at large scales which can affect the particle spectra and radiation and which are difficult to simulate yet with the PIC simulations.

6.1 PIC vs Monte Carlo simulations

At present, the self-consistent problem of relativistic collisionless shock structure, turbulence production, and particle acceleration can only be solved with full scale PIC simulations. A very important result of these simulations is the demonstration that the Weibel instability, as well as current-driven filamentation instabilities, produce short-wavelength turbulence which can scatter accelerated particles (see e.g. Spitkovsky, 2008; Plotnikov et al., 2011; Sironi et al., 2013; Lemoine et al., 2014c). In the simulations performed thus far, the self-generated turbulence is seen to fill a layer of a few hundred ion skin depths close to the subshock. Another important result is that PIC simulations have shown that, regardless of the background magnetic field geometry, unmagnetized relativistic shocks can be efficient particle accelerators converting $\gtrsim 10$% of their ram pressure into accelerated particles (e.g. Sironi et al., 2013).

Despite these successes, the limitations on box size, run time, and dimensionality (see Jones et al., 1998; Vladimirov et al., 2008) inherent in PIC simulations force the consideration of less fundamental descriptions. This is especially true for non-relativistic and trans-relativistic shocks where the precursors may extend over the scales of hundreds or thousands of gyroradii of maximal energy CRs. By parameterizing the particle transport, the Monte Carlo technique we discuss can simulate the non-linear shock structure on astrophysically significant scales for arbitrary shock speeds (see Ellison et al., 1995; Summerlin and Baring, 2012; Ellison et al., 2016). In the case of non-relativistic shocks, the Monte Carlo model can simultaneously include the non-linear shock structure, accelerated particle production, and magnetic turbulence generation accounting for both short scale and long-wavelength instabilities produced by the anisotropic CR particle distributions (e.g. Vladimirov et al., 2009; Bykov et al., 2014). It is even possible to include super-diffusive propagation (i.e. the so-called Lévy-walk models, e.g. Zimbardo and Perri, 2013) consistently in the shock precursor (Bykov et al., 2017).

In principle, microphysical parameters determined with PIC simulations can be implemented in the Monte Carlo code to allow the analysis of large-scale CR modified shocks. While this will allow for a large increase in dynamic range, it remains uncertain whether or not self-generated CR-driven instabilities occur on large scales, particularly for fully relativistic shocks (see e.g. Sagi and Nakar, 2012; Lemoine et al., 2014c). In the next section we present some Monte Carlo results showing the large-scale macroscopic structure of relativistic shocks. We also include a discussion of trans-relativistic and non-relativistic shocks to emphasize that the differences between shocks of varying speed is mainly one of analysis; mathematical techniques separate shocks more than the underlying physics which must be continuous through the trans-relativistic range.

7 Non-linear Shock Structure

It is reasonably certain from direct observations that non-relativistic shocks can inject and accelerate ions efficiently, while the electron acceleration efficiency, at least in interplanetary shocks, is still uncertain (e.g., Dresing et al., 2016). While no direct observations exist for relativistic shocks, many models, particularly those for $\gamma$-ray bursts (GRBs), assume Fermi acceleration is efficient with $\gtrsim 10\%$ of the shock ram kinetic energy placed in nonthermal particles. In many cases, this energy is assumed to be given to electrons since they produce readily observable synchrotron and inverse-Compton emission. To our knowledge, there are no self-consistent models of first-order Fermi shock acceleration that place more energy in electrons than ions. Unless shocks are accelerating a pair dominated plasma whatever energy fraction is assumed for electrons will be equaled or exceeded by ions even if they do not produce observable radiation.

If Fermi acceleration is efficient, the shock structure must be modified by the back pressure of CRs as they scatter in the shock precursor. This is illustrated in Fig. 5 which is adapted from figure 7 in Ellison et al. (2016). The shocks shown in Fig. 5 are relativistic with a Lorentz factor $\gamma_{\rm sh}=10$. For a given acceleration efficiency, non-relativistic shocks will have a more pronounced precursor structure because it is easier for CRs to scatter upstream in a non-relativistic flow. The Monte Carlo code determines the momentum and energy conserving shock structure by iteration, as described in Ellison et al. (1996, 2013).

7.1 Iteration of Shock Structure

If we assume that the shock is parallel and if we ignore the effects of escaping momentum and energy fluxes, the shock-frame, steady-state flux conservation relations are:

[TABLE]

[TABLE]

and

[TABLE]

Here, $F_{\mathrm{n}}^{0}$, $\Phi_{P0}$, and $F_{\mathrm{en}}^{0}$ are the far upstream number, momentum, and energy fluxes, respectively, and $P_{\mathrm{xx}}$ is the $xx$-component of the pressure tensor, $\mathcal{P}$. Given the input quatities, $F_{\mathrm{n}}^{0}$, $\Phi_{P0}$, and $F_{\mathrm{en}}^{0}$, these three equations have four unknowns: $n(x)$, $\gamma(x)$, $e(x)$, and $P_{\mathrm{xx}}(x)$, so an additional constraint is needed for solution. We use the adiabatic and gyrotropic equation of state:

[TABLE]

where $e(x)$ is the total energy density, $\rho(x)c(x)^{2}$ is the rest mass energy density, $\Gamma(x)$ is the adiabatic index, and $P(x)=P_{\mathrm{xx}}$ is the isotropic scalar pressure, i.e., $P=Tr(\mathcal{P})/3$ (see Double et al., 2004, for more details). All of these quantities, along with the bulk flow Lorentz factor (also called the shock profile) $\gamma(x)$, vary with position $x$ in the precursor.

The solution is accomplished in the Monte Carlo  code by iteration. Starting with a shock profile, i.e., $\gamma(x)$, the simulation is run and $F_{\mathrm{n}}(x)$, $F_{\mathrm{px}}(x)$, and $F_{\mathrm{en}}(x)$ are measured at every $x$. If the measured fluxes don’t match the known values $F_{\mathrm{n}}^{0}$, $\Phi_{P0}$, and $F_{\mathrm{en}}^{0}$ to within some error limits, $P(x)$, which is a local frame quantity, is determined using Eq. (7.1) or (51). Using this new $P(x)$, Eqs. (49)–(51) are used to predict the next shock profile, i.e., the next $\gamma(x)$. This process is iterated until the fluxes are conserved. An example from data presented in Ellison et al. (2016) is shown in Fig. 5. In contrast with semi-analytic techniques where multiple solutions can be obtained (see Drury, 1983, for a review), unique solutions always result with the Monte Carlo approach because of the non-linear (NL) feedback between the thermal injection and the shock structure.

The above procedure becomes significantly more complicated with an oblique magnetic field. In this case, there are $x$-dependent parallel and perpendicular components of the magnetic field and equations (49)—(51), are replaced by six equations which include magnetic pressure. This system has been derived by Double et al. (2004). Self-consistent solutions have been found for NL oblique non-relativistic shocks (e.g., Ellison et al., 1996), oblique relativistic unmodified shocks (e.g., Ellison and Double, 2004), but we know of no solutions, other than from PIC simulations, for NL oblique relativistic shocks.

7.2 Particle Transport and Magnetic Turbulence

Pitch-angle particle scattering and convection are determined in one of two ways in the Monte Carlo code. In the simplest scheme the mean free path $\lambda(x,p)$ is a fully parameterized function of particle momentum $p$ and position $x$ relative to the subshock at $x=0$. This method has been used extensively in previous non-relativistic applications (e.g., Berezhko and Ellison, 1999) and is currently used for relativistic shocks.

In Fig. 6 we show examples of the scattering mfp for relativistic shocks ($\gamma_{\rm sh}=10$) where the momentum dependence has been varied as follows (see Ellison et al., 2016, for full details):

[TABLE]

Here $p$ is the particle momentum in the local frame, $\alpha_{L}\leq 1$, $\alpha_{H}\geq 1$, $p_{d}$ is a dividing momentum between the low- and high-momentum ranges, $r_{g}(p_{d})=p_{d}c/(eB_{0})$ is the gyroradius for a proton with local-frame momentum $p_{d}$ in the background field $B_{0}$, and $\eta_{\mathrm{mfp}}\geq 1$ is a parameter that determines the strength of scattering. More complicated forms for $\lambda_{\mathrm{mfp}}(p)$ can be used in the Monte Carlo simulation to model semi-analytic and/or PIC results where available. Note that $\lambda_{\mathrm{mfp}}(p)$ was taken to be position independent in the results shown in Figs. 6 and 7.

The Bohm limit is described by $\eta_{\mathrm{mfp}}=\alpha_{L}=\alpha_{H}=1$, and the conditions $\alpha_{L}\leq 1$ and $\alpha_{H}\geq 1$ ensure that $\lambda_{\mathrm{mfp}}\geq\eta_{\mathrm{mfp}}r_{g}$ for all $p$.222The case with $\eta_{\mathrm{mfp}}=\alpha_{L}=\alpha_{H}=1$ is often referred to as “Bohm diffusion” however particle trajectories are calculated in the Monte Carlo code without making a diffusion approximation. In the simple geometry of a plane-parallel, steady-state shock, results scale with $\eta_{\mathrm{mfp}}$. Large values of $\eta_{\mathrm{mfp}}$ imply weak scattering with long length and time scales.

In Fig.6, $\alpha_{H}=2$ is used as a typical value predicted by the analytical theory of Sec. 4.2 of particle transport in small-scale micro-turbulence.

Regardless of the form for $\lambda_{\mathrm{mfp}}$, pitch-angle scattering is modeled as follows: after a time $\delta t\ll\tau_{g}$ ($\tau_{g}$ is the gyro-period) a particle scatters isotropically and elastically in the local plasma frame through a small angle $\delta\theta\leq\delta\theta_{\mathrm{max}}$. The maximum scattering angle in any scattering event is given by (see Ellison et al., 1990a)

[TABLE]

where $\delta t=\tau_{g}/N_{g}$, $t_{c}=\lambda_{\mathrm{mfp}}/v$ is the collision time, $v$ is the particle speed in the local frame, and $N_{g}$ is a free parameter. Typically, $N_{g}$ is chosen large enough to produce fine-scattering results that do not change substantially as $N_{g}$ is increased further (see Summerlin and Baring, 2012, for examples with values of $N_{g}$ producing large-angle scattering).

All particles are injected far upstream with a thermal distribution and scatter and convect into and across the subshock into the downstream region. A critical approximation used to date in the Monte Carlo simulation is that the subshock is transparent, i.e., no attempt is made to describe the effects of a cross-shock potential, amplified magnetic turbulence, or other effects that may occur in the viscous subshock layer. This approximation has been relaxed in recent work with the parameterization of energy transfer between ions and electrons at the subshock and an injection threshold (Warren et al., 2015a).

Upon interacting with the downstream plasma, some particles will gain energy sufficient to allow them to scatter back across the subshock and be further accelerated. The fraction of particles that are injected, i.e., do manage to re-cross the subshock is determined stochastically and constitutes a “thermal leakage injection” model once the further assumption is made that the subshock is transparent. This injection model requires no additional parameters or assumptions once the scattering mean free path is defined in equation (53).

The additional parameter implemented by Warren et al. (2015a) to model radiation spectra from unmagnetized relativistic shocks is $f_{\mathrm{ion}}$ which is the fraction of far upstream ion ram kinetic energy transferred to electrons. The far upstream kinetic energy of an ion is $(\gamma_{\rm sh}-1)m_{i}c^{2}$, where $m_{i}$ is the ion mass. When an ion crosses the subshock from upstream to downstream for the first time, a part of this energy equal to $f_{\mathrm{ion}}(\gamma_{\rm sh}-1)m_{i}c^{2}$ is removed from it. The total energy removed from all of the ions is divided equally among electrons and added to their energy as they cross the subshock into the downstream region for the first time. This parameterisation was implemented to account for the effect of the energy exchange between the ions and electrons in the subshock which is clearly seen in Fig. 11 from Sironi et al. (2013) where PIC simulations of a low magnetization ($\sigma<10^{-5}$) ion-electron plasma shock were presented.

7.3 Particle Escape

While the majority of shocked particles convect downstream and don’t re-cross the subshock to experience Fermi acceleration, a small fraction will continue to be accelerated until they run out of acceleration time or their diffusion length becomes comparable to the shock size. Both the acceleration rate and the effective shock size depend critically on self-generated turbulence which will always have a finite extent. The finite extent of self-generated turbulence guarantees that particle escape will play a role in Fermi acceleration.

In quasi-steady-state shocks, such as the Earth bow shock, or long-lived shocks, such as SNR blast waves, the effective shock size limits acceleration (except at earliest times for SNRs) and the escape of particles from the shock becomes an important part of the Fermi process. In fact, in high Mach number non-relativistic shocks, it isn’t possible to conserve energy unless escape occurs and this has been modeled with a free escape boundary (FEB). The steeper spectrum produced by Fermi acceleration in relativistic shocks obviates the energy conservation argument, but the finite extent of self-generated turbulence means particle escape will occur here as well.

Because self-generated turbulence will decay in the downstream region, particle escape can occur from downstream as well as upstream from the precursor.333In the often used approximation of an infinite, plane shock, particles do not “escape” downstream unless a boundary is imposed. Particles convect away from the shock but always have a finite, but increasingly small, probability of scattering back to the subshock. The Monte Carlo simulation limits acceleration with an upstream and/or downstream FEB, although lateral escape can also be modeled. Lateral escape may be important in objects as diverse as the Earth bow shock and GRB jets. How escape actually occurs will depend on the detailed geometry of the shock but all these scenarios produce a corresponding $p_{\mathrm{max}}$ with no differences large enough to be discerned from current observations. Even though the Monte Carlo model is steady-state, the acceleration time of individual particles is calculated and Fermi acceleration can also be limited by setting a maximum time (see Achterberg et al., 2001; Drury, 2011, for a discussions of particle escape).

An upstream FEB has been used extensively in models of CR production in young SNRs (e.g., Ellison et al., 2007; Morlino et al., 2009; Caprioli et al., 2010) where there is typically a SNR age where a transition occurs between $p_{\mathrm{max}}$ being determined by the remnant age to being determined by the finite shock radius. We note that there is direct observational evidence for upstream escape at the quasi-parallel Earth bow shock (e.g., Trattner et al., 2013).

Using Eq. (53) to vary the momentum dependence of $\lambda_{\mathrm{mfp}}(p)$ we show some examples in Fig. 7. These are shock-frame, phase-space distributions for relativistic shocks with $\gamma_{\rm sh}=10$. For these examples, the size of the precursor (i.e., the position of the upstream FEB) has been varied to produce approximately the same $p_{\mathrm{max}}$. The model parameters are described in detail in Ellison et al. (2016) but it is sufficient to say that the precursor size depends strongly on the momentum dependence of $\lambda_{\mathrm{mfp}}(p)$. For Model (J) with $\alpha_{H}=1$, the upstream FEB is at $-25\,r_{g0}$, while Model (L) with $\alpha_{H}=2$ and essentially the same $p_{\mathrm{max}}$, has a FEB at $-2\!\times\!10^{6}\,r_{g0}$. Here $r_{g0}=m_{p}u_{0}c/(eB_{0})$ where $u_{0}$ is the shock speed, $B_{0}$ is the far upstream magnetic field, and $m_{p}$, $c$, and $e$ have their using meanings. A value of $\alpha_{H}>1$ implies weak turbulence but this also means that energetic particles will be able to scatter well into the precursor opening the possibility of self generation of precursor turbulence.

7.4 Thermal Particle Injection

A collisionless shock moving through ambient material will heat all of the material and inject some fraction of that material into the Fermi acceleration mechanism. While it is clear that a background superthermal population will always be accelerated (e.g., Blandford and Ostriker, 1978), there is compelling evidence that some fraction of the cold unshocked gas, upon being heated, will be Fermi accelerated as well (e.g., Blandford and Eichler, 1987; Ellison et al., 1990b; Sironi et al., 2013). We characterize any process that turns cold, unshocked material into superthermal particles, thermal particle injection.

Thermal particle injection (TPI) has been considered a particularly difficult part of Fermi acceleration from both observational and theoretical aspects. By definition, thermal particles are always highly anisotropic at shocks of any speed, making injection difficult to describe with semi-analytic models based on the diffusion-advection equation. Injection can be treated more directly with the Monte Carlo code since there are no constraints on anisotropy. This makes it possible to describe injection in an internally self-consistent way within the assumptions made to describe particle transport with no additional assumptions.

However, as seen in PIC simulations, the subshock layer will be highly turbulent (e.g., Sironi et al., 2013; Caprioli et al., 2015) making the injection process fundamentally complex. While the simple approximations for particle scattering made by the Monte Carlo model cannot describe the full complexity of TPI, Monte Carlo results have been shown to be reasonably consistent with PIC simulations of relativistic shocks. Using similar parameters, the Monte Carlo simulation can reasonably reproduce the NL shock structure, particle spectral shape, and overall acceleration efficiency obtained by the PIC simulation (see Ellison et al., 2013, for details). We know of no semi-analytic method that can do this for relativistic shocks. In effect, the Monte Carlo injection approximations average over the complex microscopic details while ensuring that global momentum and energy conservation is maintained. These approximations allow parameterizations of PIC results that can then be applied on larger CR momentum, length, and time scales than are currently possible with PIC simulations.

7.5 Precursor Structure

In order for Fermi acceleration to work, some downstream particles must scatter back across the shock into the upstream region creating a precursor. While the length of this precursor is a strong inverse function of shock Lorentz factor, it must be present on some scale if CRs are accelerated beyond the initial boost they obtain in their first interaction with the shock.

We characterize the precursor scale with a penetration depth, $D_{\mathrm{pen}}(p_{\mathrm{peak}})$, where $p_{\mathrm{peak}}$ is the momentum where $p^{4.23}f(p)$ peaks. Examples are shown in Fig. 9 for an unmodified (UM) shock (top panel) and the NL shock (bottom panel) where momentum and energy are conserved and the shock is smoothed. While it is clear that the momentum of the peak, $p_{\mathrm{peak}}$, is not precise—the precursor distributions are not delta functions—there is a clear difference in $D_{\mathrm{pen}}(p_{\mathrm{peak}})$ between the two cases.

In Fig. 10 we plot $D_{\mathrm{pen}}$ vs. $p_{\mathrm{peak}}$ for shocks with varying Lorentz factors. In the UM shocks (top panel), $D_{\mathrm{pen}}\propto p_{\mathrm{peak}}^{2}$ in all cases. Note that the curves in the top panel are scaled to the $\gamma_{\rm sh}=10$ result with a normalization factor $(\gamma_{\rm sh}/10)^{2}$ as shown. For a given $p$, the penetration depth scales as $1/\gamma_{\rm sh}^{2}$, as expected for an UM flow. The NL shocks show something very different. For $\gamma_{\rm sh}=10$, 30, and 60, $D_{\mathrm{pen}}\propto p_{\mathrm{peak}}^{1/2}$. It is much easier for low momentum particles to diffuse upstream once the precursor develops. At the highest $p$ shown in Fig. 10, the penetration depth becomes similar for the UM and NL cases. Note that the trans-relativistic UM $\gamma_{\rm sh}=2$ shock obeys the $D_{\mathrm{pen}}\propto p_{\mathrm{peak}}^{2}$ scaling but doesn’t conform to the $1/\gamma_{\rm sh}^{2}$ scaling.

Fermi acceleration in relativistic shocks has been suggested as a possible source of ultra-high-energy CRs (e.g., Keshet and Waxman, 2005), and relativistic shocks may be responsible for non-thermal radiation observed from a number of astrophysical sources including $\gamma$-ray bursts (GRBs), type Ibc supernovae, and extra-galactic radio jets. In order to explain this CR production and radiation, electron and ion acceleration to ultra-relativistic energies must be modeled. Particle-in-cell simulations concentrate on microphysical aspects of relativistic shocks but lack the dynamic range needed to model nonlinear shock structure and Fermi acceleration from thermal injection energies to CR energies above $\sim 10^{15}$ eV. On the other hand, semi-analytic descriptions of magnetic field amplification and particle scattering can cover a large dynamic range but have not yet included non-linear effects from efficient Fermi acceleration. We believe the Monte Carlo technique offers a way to combine insights gained from semi-analytic studies with PIC results, such as the transfer of ion energy to electrons, in a NL, broadband model capable of producing radiation that can be compared with observations (e.g., Warren et al., 2015b, 2017).

7.6 Lessons from non-relativistic shock models:

Scattering Center Speed

Thus far, the calculation of self-generated turbulence in NL relativistic shocks has only been done with PIC simulations. Semi-analytic calculations have been performed in unmodified relativistic shocks, as discussed earlier. As a prelude to future relativistic shock Monte Carlo calculations, we show here some examples for NL non-relativistic shocks performed recently by Bykov et al. (2014). Despite the physical differences between relativistic and non-relativistic shocks discussed above, the non-relativistic results can serve as a useful guide.

The NL, non-relativistic Monte Carlo model includes magnetic field amplification from resonant instabilities, as well as short-wavelength (e.g., Bell, 2004; Pelletier et al., 2006; Amato and Blasi, 2009) and long-wavelength (e.g., Bykov et al., 2012a; Schure et al., 2012) non-resonant instabilities. In Fig. 11, taken from Bykov et al. (2014), we show downstream particle spectra (top panel) for examples with and without the long-wavelength, non-resonant instability. The examples shown were simulated with and without the turbulence cascade, as indicated. It is important to note that in these non-relativistic results, the scattering mean free path is self-consistently determined and has an arbitrary $p$ and $x$ dependence. The self-consistently determined diffusion coefficient is shown in the bottom panel in Fig. 11, again taken from Bykov et al. (2014). It is clear from Fig. 11 that the mfp does not have a simple Bohm-like momentum dependence. A very important distinctive feature of these non-relativistic results is the spectral hardness at high $p$. Such hard spectra are not expected in relativistic shocks, as shown in Fig. 8. The hard CR spectrum in NL, non-relativistic shocks provides the free energy needed for strong amplification of magnetic turbulence in the shock precursor. It remains to be seen how turbulence will generated in relativistic shocks.

Self-generated magnetic turbulence is essential for Fermi acceleration to work since particles must scatter nearly elastically off this turbulence to cross the subshock several times and be accelerated above thermal energies. An important aspect of this is that the turbulence may move through the bulk plasma so the effective compression ratio for Fermi acceleration, $R_{\mathrm{eff}}(x)=v_{\mathrm{scat}}(x)/u_{2}$, may differ from the compression ratio of the bulk plasma, $R=u(x)/u_{2}$, where $u_{2}$ is the downstream flow speed in the shock rest frame.

Fig. 12 shows an example from a non-relativistic shock calculation where $v_{\mathrm{scat}}$ has been determined self-consistently with the Monte Carlo simulation. While most semi-analytic calculations assume $v_{\mathrm{scat}}$ is the Alfvén speed, $v_{a}(x)$, calculated either with the ambient magnetic field or the amplified field, the self-consistent results obtained in Bykov et al. (2014) show that it can differ substantially from $v_{a}(x)$ regardless of the magnetic field assumed.

The velocity of the scattering centers, i.e. the Weibel frame, has been determined analytically in the relativistic regime in Sec. 4.2. This motion is of course implicitly included self-consistently in PIC simulations, and it will be incorporated in future Monte Carlo simulations of relativistic shocks.

7.7 Radiation from Monte Carlo shock modeling

Within the assumptions of the TPI model, the Monte Carlo technique calculates the thermal and nonthermal plasma consistently. Therefore, the distribution functions, such as those shown in Fig. 7, contain the background plasma as well as the accelerated population. This means the ion and electron spectra are absolutely normalized for a given set of environmental and model parameters and the synchrotron, bremsstrahlung, and pion-decay emission is also absolutely normalized without requiring additional parameters. The current Monte Carlo code does not include synchrotron-self-Comptonization so an external photon field is required for inverse Compton. Additional parameters will also be required for the other radiation mechanisms if radiation is calculated from shock accelerated particles interacting with an external medium.

In Fig. 8, distributions from Warren et al. (2015a) are shown where the fraction of energy transferred from ions to electrons, $f_{\mathrm{ion}}$, has been set to two different values as indicated. In both cases, the NL shock structure was calculated including electrons, protons, and He2+. Fig. 13, also from Warren et al. (2015a), shows the broadband continuum radiation produced in these shocks. The transfer of ion energy to electrons radically influences the synchrotron emission. The emission at $\sim 1$ Mev with $f_{\mathrm{ion}}=0.4$ is $\sim 100$ times that for $f_{\mathrm{ion}}=0.1$. On the other hand, the pion-decay emission in the GeV-TeV region drops for $f_{\mathrm{ion}}=0.4$ since the ions have less energy. See Warren et al. (2015a) for a full discussion of these effects.

8 Summary

In this article, we have provided an analytical description of the main components of a relativistic collisionless shock, building on earlier work and borrowing from the results of particle-in-cell simulations. We have completed this discussion with a review of Monte Carlo techniques, which offer a complementary method to approach the mesoscale physics of these shock waves. As an executive summary, we review here some of the most salient points of the physics of relativistic collisionless shock waves:

(i) In unmagnetized shocks where the background magnetic field is dominated by self-generated turbulence, the shock forms thanks to a partial reflection of the incoming particles on a front produced by the enhancement of the electromagnetic energy density.

(ii) The electromagnetic precursor is generated by a streaming micro-instability. The Weibel instability emerges as the leading micro-instability in the unmagnetized regime. At weak but finite magnetization, a current driven instability is expected to play an important role as well.

(iii) At magnetization levels not far below unity (and higher), the shock is formed by the compression of the background magnetic field.

(iv) As we have argued, the Weibel turbulence which is excited in the precursor propagates forwards at a sub-relativistic speed (typically $\xi_{\rm cr}$ with respect to the background plasma, $\xi_{\rm cr}$ being the ratio of the cosmic ray pressure over the incoming energy density in the front frame).

(v) Electron heating in the precursor occurs through two processes: through the motion of the scattering centers with respect to the background, and through the induction electric field.

(vi) The motion of the Weibel frame with respect to the background, in addition to its heating effect, produces a slowing of the incoming plasma that might reduce its speed to a mildly relativistic one, thereby allowing a saturation of the Weibel instability by the dispersion effect of the beam.

(vii) Weibel turbulence gives rise to filaments in the background plasma, whereas the beam of superthermal particles does not suffer filamentation due to its high inertia.

(viii) The particle scattering rate decreases with the energy of the considered particle in proportion to $\gamma^{-2}$. This result first established with a random phase approximation has been confirmed by a Hamiltonian chaos investigation, and has also been confirmed by PIC simulations.

(ix) Because of this fast decrease of the scattering frequency, even a tiny mean field can stop the Fermi acceleration process; acceleration can take place in shock waves for which the magnetization $\sigma\,\lesssim\,\xi_{B}^{2}$, i.e., $\sigma\,\lesssim\,10^{-4}$. In particular, if the Weibel instability dominates turbulence production on all scales into the non-linear regime, protons cannot be accelerated beyond $10^{16}\,$eV in the termination shocks of gamma-ray bursts. Electrons can be efficiently accelerated up to $10\,$TeV or so in similar conditions and give rise to GeV synchrotron emission at early times.

(x) If relativistic shocks are as efficient in producing superthermal particles as suggested by PIC simulations, and often assumed in astrophysical models, non-linear effects from efficient CR production must be considered.

(xi) With efficient CR production, the shock structure, turbulence production, Fermi acceleration, and particle injection from the thermal background are coupled. By using phenomenological particle scattering, Monte Carlo methods can address this multi-scale coupling over large dynamic ranges in an internally consistent fashion.

(xii) The analysis of relativistic shocks is substantially more difficult than non-relativistic ones because particle distributions are not nearly isotropic in all frames and the diffusion approximation cannot be used directly. Nevertheless, the underlying physics of collisionless shock formation, particle injection, wave generation, and Fermi acceleration must be continuous from non-relativistic speeds through the trans-relativistic regime. Monte Carlo techniques can address important non-linear aspects of this transition.

(xiii) A major direction for future work is to include magnetic field amplification for relativistic and trans-relativistic shocks. Additional work can be directed toward non-linear oblique shocks and using the Monte Carlo code in astrophysical applications.

Appendix A Appendices of the first part of the paper

A.1 Determination of the Weibel frame in a pair plasma

We consider the transverse perturbed motion in the linear approximation in the background plasma and in the relativistic beam of returning particles, assuming here a pair shock. For a species indexed by $\alpha$, moving at velocity $\beta_{\alpha}$ in a frame to be specified later, with proper enthalpy density $h_{\alpha}$ and number density $n_{\alpha}$, the perturbed dynamical and continuity equations can be written as

[TABLE]

It is understood that $\delta p_{\alpha}$, $\boldsymbol{\delta u_{\alpha}}$, $\delta n_{\alpha}$ are first order perturbations while $h_{\alpha}$, $n_{\alpha}$ etc. are zeroth-order uniform and constant in time variables. The above implicitly assumes that the electromagnetic perturbations are characterized by $\delta A^{0}$ and $\delta A^{x}$ only, and that $\partial_{x}\cdot\,=\,0$ for all quantities.

For the sake of simplicity, we assume cold species and neglect the pressure perturbation. The more general case of warm plasmas is discussed further below. Then for each fluid (one fluid being composed of $+$ and $-$ species), meaning the superthermal fluid and the background plasma fluid, we derive from the above two equations an uncoupled equation for the evolution of its apparent charge perturbation:

[TABLE]

where apparent charge density means $\rho_{\alpha}\,\equiv\,\gamma_{\alpha}e(n_{\alpha+}-n_{\alpha-})$; $\omega_{\rm p,\alpha}$ represents the plasma frequency defined in the rest frame of the considered fluid, such that

[TABLE]

with the index $\alpha$ running over the superthermal (returning) particles and over the background plasma. For a plasma temperature $T\,\ll\,mc^{2}$ (cold non-relativistic plasma), $h_{\alpha}\,\simeq\,n_{a}mc^{2}$ and $\omega_{\rm p\alpha}^{2}=4\pi n_{\alpha}e^{2}/m$. This is the case of the background plasma, for which $n_{\alpha}$ designs the density $n_{\rm u}$ of the ambient plasma. For a plasma of relativistic temperature $T\,>\,mc^{2}$, $h_{\alpha}\,\simeq\,4p_{a}\,\simeq\,4/3\,\Delta\gamma_{\alpha}n_{\alpha}mc^{2}$ and $\omega_{\rm p\alpha}^{2}=3\pi n_{\alpha}e^{2}/\left(\Delta\gamma_{\alpha}m\right)$. This is the case of the plasma of returning particles for which: $n_{\alpha}\,\simeq\,\xi_{\rm cr}\gamma_{\rm sh}n_{\rm u}$, while the thermal Lorentz factor is $\Delta\gamma_{\alpha}\,\sim\,\gamma_{\rm sh}$ (front frame); thus its plasma frequency is such that $\omega_{\rm pb}\,\simeq\,\xi_{\rm cr}\omega_{\rm p}$, $\omega_{\rm p}$ the plasma frequency of the background plasma.

We define the Weibel frame as that in which the instability becomes purely magnetic, i.e. where the electrostatic potential and the total electric charge density vanish. The evolution of the total charge density is governed by

[TABLE]

In the Lorentz gauge, $\delta\Phi$ can be expressed directly in terms of the total charge perturbation; switching over to Fourier space, one infers

[TABLE]

Thus there is no source of electrostatic excitation when $\sum_{\alpha}\beta_{\alpha}\omega_{p\alpha}^{2}=0$; this defines the Weibel frame, that always exists. For the case we are treating, the relation defining the Weibel frame is simply:

[TABLE]

with $\beta_{\rm b|w}$ (resp. $\beta_{\rm bg|w}$) the velocity of the beam (resp. the background plasma) relative to this Weibel frame. The beam moves relativistically with respect to the background plasma, therefore either $\beta_{\rm b|w}\,\sim\,1$ or $\beta_{\rm bg|w}\,\sim\,-1$; we will argue further on that the plasma filamentation is carried by the background, which motivates the former choice. Thus the Weibel frame moves forward (in the same direction as the shock front) relative to the background plasma, but at a low velocity

[TABLE]

One can also directly calculate the speed of the Weibel frame with respect to the front, written $\beta_{\rm w}$ (with other velocities are also written in the front frame):

[TABLE]

which leads to

[TABLE]

It can be easily checked that the latter result is also obtained by a Lorentz transform of the result in the background frame.

A.2 Determination of the Weibel frame in a ionic plasma with hot electrons

We assume now a plasma made of protons and hot electrons. When electrons are relativistically hot, they are in a quasi-Boltzmann equilibrium and also they are responsible for strong Landau effects in the electromagnetic response functions. The transverse and longitudinal response functions of background electrons are defined through

[TABLE]

with

[TABLE]

and let us recall the notation: $\xi_{\rm th}\,\equiv\,kT_{e}/(m_{p}c^{2})$. Using a cold response for the background ions in their rest frame, i.e. $\delta j_{x,i}\,=\,-\omega_{\rm pi}^{2}/(4\pi)\delta A_{x|\rm u}$, inserting these responses in the equation for $\delta A_{x|\rm u}$, one obtains the dispersion relation:

[TABLE]

This dispersion relation neglects terms of order $\omega^{2}/k^{2}$ in the second term on the l.h.s., and it neglects a contribution of the electrostatic potential on the r.h.s, which is of order $\beta_{\rm w|u}\,\sim\,\xi_{\rm cr}$ relative to the above term. The effective sound velocity $\beta_{*}c$ which enters the equation of motion for the charge perturbation of the beam is of order $1/\gamma_{\rm b}$, i.e. $\beta_{*}^{2}c^{2}=c_{\rm s,b}^{2}/\gamma_{\rm sh}^{2}$ with $c_{\rm s,b}=1/\sqrt{3}$.

Let us write the complex frequency of the Weibel modes as follows: $\omega=iz\xi_{\rm th}kc$, where $z\in\mathbb{C}$. The modification of the transverse response by the electronic Landau effect is the main modification of the dispersion relation for Weibel modes, at least off the saturation regime, that reads (omitting terms of order $1/\gamma_{\rm sh}^{2}\xi_{\rm cr}^{2}$):

[TABLE]

where $a\equiv 1+k^{2}c^{2}/\omega_{\rm pi}^{2}$. The physical situation of interest corresponds to $\xi_{\rm cr}\gg\xi_{\rm th}$ where

[TABLE]

We thus find an instability with a growth rate:

[TABLE]

This growth rate is not significantly different from the cold one for $kc=\omega_{\rm pi}$ and $\xi_{\rm th}\,\sim\,\xi_{\rm cr}$. The threshold due to beam dispersion is found by keeping all the terms except the third degree term and assuming $\xi_{\rm cr}-a\beta_{*}^{2}$ positive but arbitrary small. Then we find the unstable root:

[TABLE]

which leads to the saturation condition: $\gamma_{\rm bg}=\xi_{\rm cr}^{-1/2}$ in the front frame.

A.3 The scattering law derived from the Hamiltonian chaos

We thus consider returning cosmic rays in the Weibel frame with a very high Lorentz factor $\gamma$, such that their motion in the $x-$direction is almost constant ($v_{x}=c$ with a very good accuracy); then, only their non-relativistic transverse motion matters with $v_{\perp}\,\sim\,c/\gamma_{\rm sh}$. We assume that the magnetic modes carry an $x-$wave number $k_{0}$ which is small compared to the transverse ones $k_{0}\,\ll\,k_{\perp}\,\sim\,\delta_{i}^{-1}$, as it appears to be in numerical simulations. Thus the cosmic rays can suffer resonances with Weibel modes: $k_{\perp}v_{\perp}\pm k_{0}c=0$.

The most simple description consists in considering the transverse dynamics only, the contribution of the x-coordinate amounting to a time dependent driving with $k_{0}x(t)=k_{0}ct$, with a frequency $\omega_{0}\,\equiv\,k_{0}c\,\ll\,\omega_{pi}$. In a 2D description, with $y$ denoting the transverse coordinate, this simple description involves a Lagrangian of the form:

[TABLE]

where the small amplitude parameter is defined as

[TABLE]

The corresponding Hamiltonian describes a Chirikov standard map, with the small $\epsilon$ parameter characterizing the perturbation which will provide the transition to chaotic dynamics.

Let us define a local coordinate of phase space $\theta_{n}\,\equiv\,k_{n}y(t)\mp\omega_{0}t$. A resonance is defined by $\dot{\theta}_{n}=0$. Around such a resonance, the Lagrangian can be approximated as the time average:

[TABLE]

The conjugate momentum of $\theta_{n}$ is

[TABLE]

and the corresponding Hamiltonian is:

[TABLE]

This time independent Hamiltonian describes the energy invariant around each resonance. The $n-$resonance occurs for $J_{n}\,=\,J_{n}^{(r)}\,=\,\pm\omega_{0}/k_{n}^{2}$ (which is nothing but $k_{n}v_{n}=\pm\omega_{0}$). Assuming for the time being that the resonances are sufficiently isolated, then setting $J_{n}\,=\,J_{n}^{(r)}+\delta J_{n}$ around such a resonance, one finds a pendulum like Hamiltonian, with oscillations described by

[TABLE]

The oscillations in $\theta_{n}$ have a pulsation $\Omega_{n}$ in the bottom of the well such that:

[TABLE]

The separatrix being defined as the region which separates oscillations from rotations, one derives the maximum amplitude of the separatrix around the $n-$resonance in the $J_{n}$ direction of phase space as $\Delta J_{n}=2\sqrt{2}\Omega_{n}/k_{n}^{2}$. In contrast, the distance between two neighboring resonances is $\omega_{0}\left(k_{n}^{-2}-k_{n-1}^{-2}\right)\,\simeq\,2\omega_{0}\delta k_{n}/k_{n}^{3}$. The simple heuristic threshold for chaos, proposed by Chirikov (1979), is that the separatrix amplitude be larger than the spacing between resonances, giving: $\Omega_{n}>\omega_{0}\delta k_{n}/k_{n}^{2}$, which is always easily satisfied in our case, the transverse spectrum of Weibel instability being a continuum.

Therefore, above a very low threshold, the resonances overlap and the transverse motions of cosmic rays randomly jump from one resonance to another one. Using Eq. (79), one infers the evolution of $J_{n}$ after $i$ jumps:

[TABLE]

with $\Delta t\,\simeq\,\omega_{0}^{-1}$ the time taken to change resonance, which plays the role of coherence time. The above assumes a random $\theta_{n}$ in resonance $n$ at step $i$. Now, the variation in $J_{n}$ is related to the variation of the pitch angle $\alpha$ of the particle, since $\Delta J_{n}\,=\,\Delta v_{\perp}/k_{n}\,\simeq\,\Delta\alpha/k_{n}$. As $t\,=\,i\,\Delta t\,\rightarrow+\infty$, one thus infers a scattering frequency

[TABLE]

Acknowledgements.

G. Pelletier and M. Lemoine were financially supported by the Programme National Hautes Énergies (PNHE) of the C.N.R.S. and by the ANR-14-CE33-0019 MACH project. A. M. Bykov was supported by RSF grant 16-12-10225.