Physics of relativistic collisionless shocks: The scattering center frame

TL;DR

This paper investigates the microphysics of unmagnetized relativistic collisionless shocks, focusing on the Weibel instability and introducing the concept of the scattering center (Weibel) frame, with theoretical and simulation-based analysis.

Contribution

It introduces and characterizes the Weibel frame in relativistic shocks, providing kinetic and nonlinear models that describe its velocity relative to plasma and shock front.

Findings

The Weibel frame moves at subrelativistic speeds relative to the plasma.

Its velocity scales with the beam pressure parameter .

Simulation results agree with theoretical predictions.

Abstract

In this first paper of a series dedicated to the microphysics of unmagnetized, relativistic collisionless pair shocks, we discuss the physics of the Weibel-type transverse current filamentation instability (CFI) that develops in the shock precursor, through the interaction of an ultrarelativistic suprathermal particle beam with the background plasma. We introduce in particular the notion of "Weibel frame", or scattering center frame, in which the microturbulence is of mostly magnetic nature. We calculate the properties of this frame, using first a kinetic formulation of the linear phase of the instability, relying on Maxwell-J\"uttner distribution functions, then using a quasistatic model of the nonlinear stage of the instability. Both methods show that: (i) the Weibel frame moves at subrelativistic velocities relative to the background plasma, therefore at relativistic velocities relative to the shock front; (ii) the velocity of the Weibel frame relative to the background plasma scales with $\xi_{\rm b}$, i.e., the pressure of the suprathermal particle beam in units of the momentum flux density incoming into the shock; and (iii), the Weibel frame moves slightly less fast than the background plasma relative to the shock front. Our theoretical results are found to be in satisfactory agreement with the measurements carried out in dedicated large-scale 2D3V PIC simulations.

Equations

$$p_{\mathrm{b}}={\xi }_{\mathrm{b}}{F}_{\infty }.$$

$$w_{\alpha }{u_{\alpha }}^{\mu }{\partial }_{\mu }{u_{\alpha }}^{\nu }+{h_{\alpha }}^{\mu \nu }{\partial }_{\mu }{p}_{\alpha }=q_{\alpha }{n}_{\alpha }{u_{\alpha }}_{\mu }{F}^{\nu \mu },$$

$$w_{\alpha }{u_{\alpha }}^{\mu }{\partial }_{\mu }{\delta u_{\alpha }}^{\nu }+{h_{\alpha }}^{\mu \nu }{\partial }_{\mu }\delta p_{\alpha }=q_{\alpha }{n}_{\alpha }{u_{\alpha }}_{\mu }\delta F^{\nu \mu }.$$

$${\gamma }_{\alpha }^2{w}_{\alpha }{\partial }_t{\delta {\beta }_{\alpha }}^x+{\beta }_{\alpha }{\partial }_t\delta p_{\alpha }$$

$${\gamma }_{\alpha }^2{w}_{\alpha }{\partial }_t\delta {\mathbf{\beta }}_{\perp \mathbf{\alpha }}+{\nabla }_{\perp }\delta p_{\alpha }$$

$${\gamma }_{\alpha }{\partial }_t\delta n_{\alpha }+n_{\alpha }{\partial }_t\delta {\gamma }_{\alpha }+{\gamma }_{\alpha }{n}_{\alpha }{\nabla }_{\perp }\cdot \delta {\mathbf{\beta }}_{\perp \mathbf{\alpha }}=0.$$

$$\delta p_{\alpha }={\hat{\Gamma }}_{\alpha }\frac{p_{\alpha }}{n_{\alpha }}\delta n_{\alpha }.$$

$$c_{\alpha }^2\equiv {\hat{\Gamma }}_{\alpha }\frac{p_{\alpha }}{w_{\alpha }},$$

$${\partial }_t\delta p_{\alpha }=\frac{c_{\alpha }^2}{{\gamma }_{\alpha }^2\left(1-c_{\alpha }^2{\beta }_{\alpha }^2\right)}\left[q_{\alpha }{\gamma }_{\alpha }{n}_{\alpha }{\beta }_{\alpha }{\partial }_t\delta A^x+\frac{{\gamma }_{\alpha }{w}_{\alpha }}{n_{\alpha }}{\partial }_t\delta N_{\alpha }\right].$$

$$c_{\mathrm{eff}\alpha }^2=\frac{c_{\alpha }^2}{{\gamma }_{\alpha }^2\left(1-c_{\alpha }^2{\beta }_{\alpha }^2\right)}.$$

$${\gamma }_{\alpha }^2{w}_{\alpha }{\partial }_t\delta {\mathbf{\beta }}_{\perp \mathbf{\alpha }}+$$

$$-q_{\alpha }{\gamma }_{\alpha }{n}_{\alpha }\left({\nabla }_{\perp }\delta \Phi -{\beta }_{\alpha }{\nabla }_{\perp }\delta A^x\right).$$

$${\partial }_t^2\delta N_{\alpha }-$$

$$\frac{q_{\alpha }{n}_{\alpha }^2}{w_{\alpha }}\left[△\delta \Phi -{\beta }_{\alpha }\left(1-c_{\mathrm{eff}\alpha }^2\right)△\delta A^x\right],$$

$$\delta {\rho }_{\alpha }=\frac{{\Omega }_{\mathrm{p}\alpha }^2}{4\pi }\frac{k^2}{{\omega }^2-c_{\mathrm{eff}\alpha }^2{k}^2}\left[\delta \Phi -{\beta }_{\alpha }\left(1-c_{\mathrm{eff}\alpha }^2\right)\delta A^x\right],$$

$${\Omega }_{\mathrm{p}\alpha }^2\equiv \frac{4\pi n_{\alpha }{q}_{\alpha }^2}{w_{\alpha }/n_{\alpha }}.$$

$${\delta j_{\alpha }}^x={\rho }_{\alpha }{\delta {\beta }_{\alpha }}^x+{\beta }_{\alpha }\delta {\rho }_{\alpha }.$$

$$\displaystyle{\delta j_{\alpha}}^{x}=\frac{\Omega_{\rm p\alpha}^{2}}{4\pi}\biggl{\{}$$

$$\displaystyle-\left[1-\beta_{\alpha}^{2}\left(1-c_{\rm eff\,\alpha}^{2}\right)\frac{\omega^{2}-k^{2}}{\omega^{2}-c_{\rm eff\,\alpha}^{2}k^{2}}\right]\delta A^{x}\biggr{\}}\,.$$

$$\sum _{\alpha }{\Omega }_{\mathrm{p}\alpha }^2{\beta }_{\alpha |\mathrm{w}}\frac{1-c_{\mathrm{eff}\alpha }^2}{{\zeta }^2-c_{\mathrm{eff}\alpha }^2}=0,$$

$$n_{\mathrm{p}}{\beta }_{p|w}\left(1+\frac{{\hat{\Gamma }}_{\mathrm{p}}{T}_{\mathrm{p}}}{{\gamma }_{p|w}^2{\zeta }^2}\right)+\frac{n_{\mathrm{b}}{\beta }_{b|w}}{T_{\mathrm{b}}}\left(1-\frac{1}{{\hat{\Gamma }}_{\mathrm{b}}}\right)\simeq 0,$$

$${\beta }_{w|p}\equiv -{\beta }_{p|w}\simeq \frac{n_{\mathrm{b}}{\beta }_{b|w}}{n_{\mathrm{p}}{T}_{\mathrm{b}}}\left(1-\frac{1}{{\hat{\Gamma }}_{\mathrm{b}}}\right)$$

$${\beta }_{w|p}\simeq {\beta }_{b|w}{\xi }_{\mathrm{b}}\frac{{\gamma }_{\mathrm{p}}}{{\gamma }_{\infty }}\frac{{\hat{\Gamma }}_{\mathrm{b}}-1}{{\kappa }_{T_{\mathrm{b}}}^2{\hat{\Gamma }}_{\mathrm{b}}}.$$

$${\beta }_{w|p}\simeq {\xi }_{\mathrm{b}}\frac{{\gamma }_{\mathrm{p}}}{{\gamma }_{\infty }}\frac{{\hat{\Gamma }}_{\mathrm{b}}-1}{{\kappa }_{T_{\mathrm{b}}}^2{\hat{\Gamma }}_{\mathrm{b}}}.$$

$${\beta }_{w|p}\simeq \frac{1}{4{\kappa }_{T_{\mathrm{b}}}^2}{\xi }_{\mathrm{b}}.$$

$${\beta }_{\mathrm{w}}=\frac{{\beta }_{w|p}+{\beta }_{\mathrm{p}}}{1+{\beta }_{w|p}{\beta }_{\mathrm{p}}}\simeq {\beta }_{\mathrm{p}}\left(1-\frac{1}{4{\kappa }_{T_{\mathrm{b}}}^2}\frac{{\xi }_{\mathrm{b}}}{{\gamma }_{\mathrm{p}}^2}\right).$$

$${\delta j_{\mathrm{p}}}_{|w}^x\simeq \frac{{\omega }_{\mathrm{p}}^2}{4\pi }\delta A_{|w}^x$$

$$\left({\omega }^2-k^2-{\omega }_{\mathrm{p}}^2\right)$$

$$\left({\omega }^2-k^2{c}_{\mathrm{eff}\mathrm{b}}^2\right)\left({\omega }^2-k^2-{\omega }_{\mathrm{p}}^2\right)\simeq {\Omega }_{\mathrm{pb}}^2{\beta }_{b|w}^2{k}^2,$$

$${\omega }^2\simeq k^2{c}_{\mathrm{eff}\mathrm{b}}^2-{\Omega }_{\mathrm{pb}}^2{\beta }_{b|w}^2\frac{k^2}{k^2+{\omega }_{\mathrm{p}}^2}.$$

$${\epsilon }_{yy}({\epsilon }_{xx}-1/{\zeta }^2)={\epsilon }_{xy}^2,$$

$${\epsilon }_{ij}(\omega ,k)=$$

Full text

I Introduction

I.1 Motivations and objectives

The Weibel-type current filamentation instability (CFI) that develops in anisotropic plasma flows Weibel (1959); Fried (1959); Davidson et al. (1972); Califano et al. (1997); Achterberg and Wiersma (2007), has gained a lot of attention in past decades because of its relevance to various fields of physics. Through the generation of skin-depth-scale current filaments surrounded by toroidal magnetic fields, which gradually isotropize the counterstreaming plasmas, this instability is a key process in the formation of astrophysical collisionless shock waves Moiseev and Sagdeev (1963). Because it also arises in the precursor of weakly magnetized shock waves, i.e., the upstream region where the unshocked plasma flows against a beam of suprathermal particles, it further sustains the shock transition and likely generates the magnetized turbulence that is required for particle acceleration and radiation. Such microphysics may well underpin the outstanding phenomenon of gamma-ray burst afterglows, in which a substantial fraction of the $\sim 10^{52}\,$ergs liberated by the cataclysmic event is radiated away by shock-accelerated electrons, on hours to months timescales Medvedev and Loeb (1999); Gruzinov and Waxman (1999).

In laser-driven high-energy-density physics and laboratory astrophysics, the CFI also plays a central role Sentoku et al. (2003); Adam et al. (2006); Allen et al. (2012); Debayle et al. (2010); Masson-Laborde et al. (2010); Mondal et al. (2012); Quinn et al. (2012); Fiuza et al. (2012); Ruyer et al. (2015a). In particular, it has been recently observed to grow in the interpenetration of laser-ablated plasmas Fox et al. (2013); Huntington et al. (2015), marking a first step towards the generation of collisionless shocks in the laboratory Drake and Gregori (2012); Chen et al. (2015); Lobet et al. (2015); Ruyer et al. (2016); *Ruyer_2017.

The development of the CFI in symmetric colliding plasmas has been studied theoretically in both the linear  Medvedev and Loeb (1999); Wiersma and Achterberg (2004); Lyubarsky and Eichler (2006); Achterberg and Wiersma (2007) and nonlinear Medvedev et al. (2004); Milosavljevic and Nakar (2006); Achterberg et al. (2007); Bret et al. (2013); Ruyer et al. (2015b); Vanthieghem et al. (2018) regimes, and numerically using particle-in-cell (PIC) simulations, Davidson et al. (1972); Lee and Lampe (1973); Silva et al. (2003); Frederiksen et al. (2004); Jaroschek et al. (2005); Nishikawa et al. (2009); Shvets et al. (2009); Bret et al. (2010a). In such a symmetric configuration, the lab frame appears as a privileged frame in which to discuss the physics of the CFI and it is found, indeed, that the CFI generates an essentially magnetic structure in that frame.

However, in the precursor of relativistic collisionless shocks, the CFI develops in a highly asymmetric configuration: in the reference frame in which the shock front lies at rest (the shock rest frame $\mathcal{R}_{\rm s}$), the background plasma appears cold, and streams at relativistic velocities through a quasi-isotropic, ultrarelativistically hot gas of suprathermal particles. Conversely, in the background plasma rest frame $\mathcal{R}_{\rm p}$, the suprathermal particles form a dense beam of angular dispersion $\ll 1$, carrying high inertia, and moving through a tenuous gas of sub- or mildly- relativistic temperature. The growth of the CFI in the precursor of relativistic collisionless shocks has been discussed in a number of numerical studies Kato (2007); Spitkovsky (2008); Martins et al. (2009); Keshet et al. (2009); Sironi and Spitkovsky (2009); Sironi et al. (2013); Haugbølle (2011); Kumar et al. (2015), but only a few theoretical studies have discussed its properties in such conditions Lemoine and Pelletier (2010); Rabinak et al. (2011); Lemoine and Pelletier (2011); Shaisultanov et al. (2012). In particular, the notion of a frame in which this CFI is essentially of magnetic nature has not, to the best of our knowledge, received attention so far; the only exception being Ref. Ruyer et al. (2016); *Ruyer_2017 where this frame was introduced to simplify the calculation of the CFI in the precursor of a nonrelativistic electron-ion shock.

The present paper is the first of a series in which we discuss the microphysics of unmagnetized, relativistic collisionless pair shocks. Here, we address in detail this notion of a “Weibel frame”, in which the electromagnetic configuration is essentially magnetic in nature. As shown in accompanying papers of this series, and in particular Lemoine et al. (2019a), this reference frame plays a fundamental role in the physics of collisionless shock waves. In Paper II Lemoine et al. (2019b), it is shown how the noninertial character of the “Weibel frame” controls the heating and slowdown of the background plasma. In Paper III Lemoine et al. (2019c), this reference frame is used to calculate the scattering rate of suprathermal particles; the latter quantity controls the residence time of suprathermal particles in the upstream, hence the acceleration timescale, thus the maximum energy etc. Therefore, this notion of a “Weibel frame” has direct phenomenological consequences and potential observable radiative signatures. Finally, in Paper IV Vanthieghem et al. (2019), we discuss the growth of the filamentary turbulence in the precursor of a collisionless shock. In the following, as in all accompanying papers, we support our present theoretical findings through detailed comparisons to dedicated large-scale 2D3V PIC simulations of relativistic collisionless shocks.

The present paper is laid out as follows. Section II introduces the notion of the “Weibel frame” in a fluid model of the linear phase of the current filamentation instability. Section III then extends those results to a fully kinetic model of the CFI; the details of these calculations are provided in Apps. A, B and C. Section IV introduces the notion of the “Weibel frame” in a nonlinear model of the filamentation phase, in which the filaments are modeled in quasistatic equilibrium at all points in the precursor. Finally, Sec. III.4 and IV.2 provide detailed comparisons of those various models to dedicated PIC simulations. A summary and conclusions are provided in Sec. V. We use units such that $k_{\rm B}=c=1$. The metric signature is $(-,+,+,+)$.

I.2 Setup

In this work, we consider an unmagnetized, relativistic collisionless shock wave propagating through an electron-positron plasma. In the shock rest frame $\mathcal{R}_{\rm s}$, the background plasma is initially incoming at relativistic velocity $\beta_{\infty}\simeq-1$ (Lorentz factor $\gamma_{\infty}\gg 1$) along the $x-$axis from $+\infty$. Inside the precursor, its mean velocity is written $\beta_{\rm p}$ (Lorentz factor $\gamma_{\rm p}$). The precursor is defined as the region permeated by a beam of suprathermal particles, which were reflected off the shock surface or accelerated through a Fermi-like process, possibly up to high energies.

Our conventions are as follows: the index p refers to the background plasma, while the index b refers to the suprathermal particle population. For species $\alpha$, $n_{\alpha}$ represents the density, $w_{\alpha}$ the enthalpy density, $p_{\alpha}$ the pressure, and $T_{\alpha}$ the temperature, all measured in the species rest frame, while ${u_{\alpha}}^{\mu}=\gamma_{\alpha}\left(1,\bm{\beta_{\alpha}}\right)$ represents its four-velocity. All throughout, $\alpha$ will denote a species (‘p’ or ‘b’) and not a space-time index. Both the background plasma and the suprathermal beam are composed of electrons and positrons, so that there are actually four species. However, unless otherwise noted, we do not make any particular distinction between electrons and positrons in a given population, and therefore drop any reference to the charge. The above proper hydrodynamic quantities nevertheless correspond to one charged species (electrons or positrons) of one population (either the background plasma or the suprathermal beam). As in other papers of this series, we rely on the following frames: the shock front rest frame $\mathcal{R}_{\rm s}$, the background plasma rest frame $\mathcal{R}_{\rm p}$, the downstream rest frame $\mathcal{R}_{\rm d}$ and the “Weibel frame” $\mathcal{R}_{\rm w}$, to be determined hereafter. Quantities evaluated in one or the other reference frame are indicated with the respective subscripts |s, |p, |d and |w. By default, however, frame-dependent quantities that lack subscripts are defined in the lab frame, which coincides with the shock rest frame $\mathcal{R}_{\rm s}$. Finally, thermodynamic moments $n_{\alpha}$, $p_{\alpha}$, $T_{\alpha}$ and $w_{\alpha}$ are always defined as proper, unless otherwise stated; they are thus defined in the rest frame of species $\alpha$.

According to the fluid shock jump conditions Blandford and McKee (1976), the typical temperature of shocked particles is $T_{\rm b}=\kappa_{T_{\rm b}}\gamma_{\infty}m_{e}$, where $\kappa_{T_{\rm b}}=(\widehat{\Gamma}_{\rm b}-1)(2-\widehat{\Gamma}_{\rm b})^{1/2}\widehat{\Gamma}_{\rm b}^{-1/2}$, in terms of the adiabatic index $\widehat{\Gamma}_{\rm b}$ of the shocked gas; the latter is relativistically hot, so $\widehat{\Gamma}_{\rm b}=4/3$ in 3D, but $\widehat{\Gamma}_{\rm b}=3/2$ in 2D, which must be considered when drawing comparison with a 2D3V (2D in configuration space, 3D in momentum space) PIC simulation; thus $\kappa_{T_{\rm b}}=1/(3\sqrt{2})$ in 3D and $\kappa_{T_{\rm b}}=1/(2\sqrt{3})$ in 2D. Strictly speaking, the temperature of suprathermal particles is by definition larger than that of the shocked plasma, defined above. Yet, as this temperature of suprathermal particles always scales with $\gamma_{\infty}$, we retain the above definition, emphasizing that $\kappa_{T_{\rm b}}$ is generically larger than the above, by about an order of magnitude, so that $\kappa_{T_{\rm b}}\sim\,$a few. This is illustrated in Fig. 1, which shows the profiles of the main hydrodynamic properties of the background plasma and suprathermal beam for two reference PIC simulations, with respective Lorentz factors $\gamma_{\infty}=17$ and $\gamma_{\infty}=173$.

Our PIC simulation are described further on in Sec. III.4. Let us note here, however, that the reference frame of these 2D3V simulations (2D in configuration space, 3D in momentum space) coincide with the downstream rest frame $\mathcal{R}_{\rm d}$. Hence, a Lorentz factor $\gamma_{\infty}=17$ (resp. $\gamma_{\infty}=173$) in the $\mathcal{R}_{\rm s}$ frame corresponds to a simulation frame Lorentz factor $\gamma_{\infty|\rm d}=10$ (resp. $\gamma_{\infty|\rm d}=100$). Simulations are conducted in the $x-y$ plane, with $\widehat{\bm{x}}$ oriented along the shock normal, and $\widehat{\bm{y}}$ defining the transverse dimension.

The ultrarelativistic suprathermal particle gas is nearly isotropic in the shock rest frame, meaning a drift Lorentz factor $\gamma_{\rm b}\sim 1$ in the shock rest frame (see Fig. 1). This population is further characterized by its pressure $p_{\rm b}$, written $\xi_{\rm b}$ in units of the incoming momentum flux at infinity $F_{\infty}=\gamma_{\infty}^{2}\beta_{\infty}^{2}n_{\infty}m_{e}$:

[TABLE]

The leading instability that develops in the precursor of weakly magnetized ultrarelativistic shocks is the Weibel-type CFI described above. In principle, this instability is defined in all momentum space $(k_{\parallel},\,\mathbf{k_{\perp}})$ where $k_{\parallel}=\mathbf{k}\cdot\widehat{\bm{x}}$ and $\mathbf{k_{\perp}}=k_{\perp}\widehat{\bm{y}}$ in our PIC simulations. The purely transverse modes correspond to $k_{\parallel}\ll k_{\perp}$, while the so-called oblique modes correspond to the limit $k_{\parallel}\simeq k_{\perp}\gg\omega_{\rm p}$, with $\omega_{\rm p}=\left(4\pi n_{\infty}e^{2}/m_{e}\right)^{1/2}$ the plasma frequency of (one charged species of) the unperturbed far-upstream background plasma. Oblique modes are likely relevant far in the precursor, but most likely Landau damped once the background plasma heats up Bret et al. (2010a); Lemoine and Pelletier (2011); Shaisultanov et al. (2012). We thus assume that the transverse modes dominate in most of the precursor and restrict ourself to the limit $k_{\parallel}\rightarrow 0$ for simplicity.

Our 2D3V PIC simulations indicate that this is a rea- sonable approximation. Consider indeed Fig. 2, which plots the energy densities in electromagnetic components $\delta E_{x}$, $\delta E_{y}$ and $\delta B_{z}$ in our two reference PIC simulations, normalized to the incoming momentum flux at infinity $F_{\infty}$. The purely transverse CFI generates a transverse $\delta B_{z}$ and its associated electrostatic $\delta E_{y}$ component; by definition of the “Weibel frame” $\mathcal{R}_{\rm w}$ (see further on), $\delta E_{y|\rm w}=,0$, and hence $\delta E_{y}=\beta_{\rm w|\rm d}\delta B_{z}$ in the $\mathcal{R}_{\rm d}$ simulation frame. The growth of this instability also generates an inductive longitudinal electric field $\delta E_{x}$. Nonstrictly transverse CFI modes (with $k_{\parallel}\neq 0$) further generate an electrostatic longitudinal $\delta E_{x}$ component. Figure 2 reveals that the energy density in the transverse field components dominates that in the longitudinal electric field everywhere in the precursor in the $\mathcal{R}_{\rm d}$ frame, and that the electric field is smaller than the magnetic field in the near precursor where the ratio can be measured accurately (up to $x$ of the order of a few hundred to a thousand $c/\omega_{\rm p}$). Finally, the transverse magnetic field energy remains larger than the longitudinal electric field component when deboosted back to the $\mathcal{R}_{\rm w}$ frame. Therefore, this “Weibel frame” appears well defined in the 2D3V simulations.

Let us emphasize that the limit $k_{\parallel}\neq 0$ certainly remains of interest; as a matter of fact, a nontrivial structure along the streaming axis turns out to be a mandatory requirement to achieve pitch-angle diffusion of suprathermal particles Achterberg and Wiersma (2007); Lemoine et al. (2019c). Yet, we expect that the main features of the (already nontrivial) calculations that follow will remain valid in a limit $0<k_{\parallel}\ll k_{\perp}$. We thus consider purely transverse electromagnetic perturbations, i.e., $\delta A^{\mu}=\left(\delta\Phi,\delta A^{x},0,0\right)$, and $\partial_{x}=0$ for all quantities.

Figure 2 suggests that the ratio $\left\langle\delta E_{y}^{2}\right\rangle^{1/2}/\left\langle\delta B_{z}^{2}\right\rangle^{1/2}$, namely $\beta_{\rm w|d}$, depends on $x$, which implies that the “Weibel frame” is not globally inertial. This frame actually decelerates from the far to the near precursor, as discussed in Lemoine et al. (2019a). The noninertial nature of $\mathcal{R}_{\rm w}$ has important consequences for the physics of the shock, most notably the deceleration and heating of the background plasma, which form the focus of a subsequent paper in this series Lemoine et al. (2019b).

In the present paper, we characterize the velocity $\beta_{\rm w|d}$ (more specifically, $\beta_{\rm w|p}$) at each point $x$ of the precursor, given the physical conditions at this point. The discussion that follows thus relies on an implicit WKB-like approximation, which stipulates that the ‘Weibel frame” has time to adjust at each point to the local physical conditions on a deceleration length scale $\left|u_{\rm w}\,{\rm d}x/{\rm d}u_{\rm w}\right|$. Improving on this assumption would necessitate a proper inclusion of noninertial effects, which are characterized by the velocity profile ${\rm d}u_{\rm w}/{\rm d}x$, in the calculations that follow. However, this velocity profile is itself determined by the response of the background plasma and the suprathermal beam to the microturbulence. Such an endeavor thus represents a rather formidable task, well beyond the scope of any current study.

II The “Weibel frame” in a fluid model

Here we present a fluid derivation of the CFI and of its associated “Weibel frame”, in the context of the precursor of an unmagnetized, relativistic electron-positron shock. Although Sec. III.4 will demonstrate that, in actual relativistic shock precursors, kinetic corrections are mandatory to describe the physics of the CFI, the following fluid model retains the advantage of simplicity as well as a pedagogical virtue. A fully kinetic description of the CFI and its “Weibel frame” in the case of Maxwell-Jüttner plasma distribution functions will be provided in the forthcoming Sec. III.

To preserve covariance, we use a relativistic fluid formalism. The conservation of the total energy-momentum tensor can be written in the compact way Achterberg and Wiersma (2007):

[TABLE]

introducing ${h_{\alpha}}^{\mu\nu}=\eta^{\mu\nu}+{u_{\alpha}}^{\mu}{u_{\alpha}}^{\nu}$, which projects orthogonally to ${u_{\alpha}}^{\mu}$ (since ${u_{\alpha}}^{\mu}{u_{\alpha}}_{\mu}=-1$). The dynamical equation thus becomes, to first order in the perturbations,

[TABLE]

In the following, the four-velocity perturbation is decomposed as ${\delta u_{\alpha}}^{\mu}=\left(\delta\gamma_{\alpha},\gamma_{\alpha}^{3}{\delta\beta_{\alpha}}^{x},\gamma_{\alpha}{\delta\beta_{\alpha}}^{y},\gamma_{\alpha}{\delta\beta_{\alpha}}^{z}\right)$, with $\delta\gamma_{\alpha}=\gamma_{\alpha}^{3}\beta_{\alpha}{\delta\beta_{\alpha}}^{x}$. Using the short-hand notation $\delta\bm{\beta_{\perp\alpha}}=\left(0,{\delta\beta_{\alpha}}^{y},\,{\delta\beta_{\alpha}}^{z}\right)$ and $\bm{\nabla_{\perp}}=\left(0,\partial_{y},\,\partial_{z}\right)$, the system can be rewritten explicitly as

[TABLE]

Current conservation written to first order also yields

[TABLE]

The pressure perturbation can be related to the density perturbation through the adiabatic index $\widehat{\Gamma}_{\alpha}$:

[TABLE]

In the following, we use the isentropic sound speed squared

[TABLE]

so that $\delta p_{\alpha}=c_{\alpha}^{2}(w_{\alpha}/n_{\alpha})\delta n_{\alpha}$. Given the relation between $\delta\gamma_{\alpha}$ and ${\delta\beta_{\alpha}}^{x}$, Eq. (4) can then be used to express the pressure perturbation in terms of the perturbed apparent density $\delta N_{\alpha}\equiv\delta(\gamma_{\alpha}n_{\alpha})$:

[TABLE]

This equation involves the effective sound speed squared of the streaming plasma species,

[TABLE]

We have $c_{\rm eff\,\alpha}^{2}\simeq 5T_{\alpha}/(3\gamma_{\alpha}^{2})$ and $c_{\rm eff\,\alpha}^{2}\simeq 1/(2\gamma_{\alpha}^{2})$ in the nonrelativistic ($T_{\alpha}\ll m_{e}$, $\widehat{\Gamma}_{\alpha}\simeq 5/3$) and ultrarelativistic ($T_{\alpha}\gg m_{e}$, $\widehat{\Gamma}_{\alpha}\simeq 4/3$) thermal limits of a 3D gas, respectively. After integration, Eq. (9) provides a direct relationship between $\delta p_{\alpha}$ and $\delta A_{x}$ and $\delta N_{\alpha}$, which can be inserted in Eq. (5) to yield

[TABLE]

Finally, combining this equation with Eq (6), one obtains

[TABLE]

which determines the response of the apparent charge density to the electromagnetic perturbation: in Fourier space, with $\delta\rho_{\alpha}=q_{\alpha}\delta N_{\alpha}$, one finds

[TABLE]

where the relativistic (proper-frame) plasma frequency, $\Omega_{\rm p\alpha}$, is defined by

[TABLE]

Note that the far-upstream plasma frequency $\Omega_{\rm pp}$ coincides with $\omega_{\rm p}$ defined earlier.

The above also allows us to determine the response of the current density, ${\delta j_{\alpha}}^{x}$, which comprises both a perturbed conduction current density as well as a perturbed advection current density:

[TABLE]

Equations (4), (9) and (13) can then be combined to derive the response in Fourier space:

[TABLE]

We define the “Weibel frame” as that in which the linear instability becomes purely magnetic, i.e., the electrostatic potential and the total electric charge density vanish. From the response of the charge density, one finds that in this frame the following relation must be fulfilled:

[TABLE]

where $\zeta=\omega/k$ is the (complex) phase velocity . In the shock precursor, the background plasma, of proper density $n_{\rm p}$ and nonrelativistic temperature $T_{\rm p}\lesssim 1$, flows at a velocity $\beta_{\rm p|s}\simeq-1$ against the suprathermal beam, of proper density $n_{\rm b}$ and relativistic temperature $T_{\rm b}\gg 1$. To leading thermal corrections, the above equation can be recast in the form

[TABLE]

where we have assumed $|\zeta|^{2}\gg c_{\rm eff\,\alpha}^{2}$ for both populations, as is consistent within a hydrodynamic model (see Sec. III). This relation indicates that the velocity of $\mathcal{R}_{\rm w}$ is, in principle, mode-dependent. Yet the response of the beam proves to be more sensitive to thermal effects than that of the plasma since the inequality $\zeta^{2}\gg c_{\rm eff\,\alpha}^{2}$ implies $T_{\rm p}/\gamma_{\rm p|w}^{2}\zeta^{2}\ll 1$. As a consequence,

[TABLE]

is a good (mode-independent) approximation of the velocity of $\mathcal{R}_{\rm w}$ relative to the background plasma frame.

The beam is more conveniently characterized by its normalized pressure $\xi_{\rm b}$ and temperature $T_{\rm b}\,=\,\kappa_{T_{\rm b}}\gamma_{\infty}m_{e}$, so that $n_{\rm b}\,=\,\kappa_{T_{\rm b}}^{-1}\gamma_{\infty}\xi_{\rm b}n_{\infty}$, see Sec. I.2. Current conservation further implies $n_{\infty}/n_{\rm p}\,=\,\gamma_{\rm p}\beta_{\rm p}/(\gamma_{\infty}\beta_{\infty})$, with $\beta_{\rm p}\,\simeq\,\beta_{\infty}\,\simeq\,-1$, so that

[TABLE]

The beam moves relativistically with respect to the background plasma, therefore either $\beta_{\rm b|w}\sim 1$ or $\beta_{\rm p|w}\sim-1$. As $\xi_{\rm b}<1$, however, the former must hold, which provides the final result:

[TABLE]

In the case of negligible deceleration of the incoming plasma ($\gamma_{\rm p}\simeq\gamma_{\infty}$) and for a 3D adiabatic index $\widehat{\Gamma}_{\rm b}=4/3$,

[TABLE]

One can also directly calculate $\beta_{\rm w}$:

[TABLE]

The above indicates that: (i) the “Weibel frame” $\mathcal{R}_{\rm w}$ moves at a sub-relativistic velocity relative to the background plasma [i.e., $u_{\rm w|p}=\beta_{\rm w|p}(1-\beta_{\rm w|p}^{2})^{-1/2}<1$], and therefore at a relativistic velocity $\simeq\beta_{\rm p}$ relative to the shock front; (ii) $\beta_{\rm w|p}$ is of opposite sign to $\beta_{\rm p}$, which implies that, in magnitude, the “Weibel frame” moves slightly less fast than the background plasma relative to the shock front; (iii) the relative velocity $\beta_{\rm w|p}$ scales with $\xi_{\rm b}$.

Finally, because $\xi_{\rm b}$ is a function of distance $x$ to the shock, Eq. (23) indicates that $\beta_{\rm w}$ itself depends on $x$, i.e., the “Weibel frame” is not globally inertial, as already observed in Fig. 2. This observation has important implications with respect to the physics of the shock, which are addressed in detail in the follow-up paper Lemoine et al. (2019b).

In Sec. III, we extend the above calculations to a kinetic description. Although the expression for $\beta_{\rm w|p}$ will be found to take somewhat different values, the above three features will remain valid.

Using the above equations, it becomes straightforward to derive the dispersion relation of the CFI in this warm fluid description. Consider the CFI in $\mathcal{R}_{\rm w}$, where $\delta\Phi_{\rm|w}=0$. The response ${\delta j_{\rm p}}_{\rm|w}^{x}$ of the background plasma can be written

[TABLE]

to first order in $\xi_{\rm b}$, since the term that has been neglected here with respect to Eq. (II) is of the order of $\beta_{\rm p|w}^{2}\sim\xi_{\rm b}^{2}$. Note also that we have assumed that the background plasma temperature remains sub-relativistic over most of the precursor, as discussed in Paper II Lemoine et al. (2019b), so that $\Omega_{\rm pp}\simeq\omega_{\rm p}$. From the Maxwell equations written in the Lorentz gauge ($k^{\mu}\delta A_{\mu}=0$), i.e., $\mathop{}\!\mathbin{\Box}\delta A^{\mu}=-4\pi\sum_{\alpha}{\delta j_{\alpha}}^{\mu}$, one derives the dispersion relation in the “Weibel frame” (subscript w omitted for clarity):

[TABLE]

which, assuming $|\omega|^{2}\ll k^{2}$ and $c_{\rm eff\,b}^{2}\ll 1$, can be approximated by

[TABLE]

with solution

[TABLE]

The stabilizing effect of the beam dispersion is manifest in this equation; this effect has been noted before in Rabinak et al. (2011); Lemoine and Pelletier (2011).

III Relativistic kinetic model

In this Section, we evaluate the “Weibel frame” velocity and the local instability growth rate within a rigorous kinetic formalism. As discussed further below, kinetic effects must indeed be taken into account when considering the development of the CFI in the shock precursor. The derivation of the kinetic dielectric tensor involves rather heavy calculations, which are relegated to Apps. A, B and C. Below, we describe the general method, summarize the approximate expressions of the relevant dielectric tensor elements, and provide the estimates of the growth rate of the CFI and the associated estimate of $\beta_{\rm w}$ in various limits. The latter results, in particular, are given in Sec. III.3.

For the time being, the reference frame in which we work is left unspecified. We first recall the linear dispersion relation fulfilled by the CFI modes Silva et al. (2002):

[TABLE]

with $\zeta=\omega/k$ as the complex phase velocity. In a fully relativistic framework, the dielectric tensor elements read ($i,j\,=\,1,\,2,\,3$) Ichimaru (1973)

[TABLE]

where $v_{y}=u_{y}/(\gamma m_{e})$, $\omega_{\rm p\alpha}^{2}=4\pi n_{\alpha}e^{2}/m_{e}$ represents the nonrelativistic plasma frequency squared of species $\alpha$ ($n_{\alpha}$ represents as before the proper density) and $f_{\alpha}^{(0)}(\mathbf{u})$ the corresponding unperturbed momentum distribution function, normalized such that $\int{\rm d}^{3}u\,f_{\alpha}^{(0)}=1$. If the non-diagonal tensor element $\varepsilon_{xy}$ happens to vanish, Eq. (28) implies either $\varepsilon_{yy}=0$ or $\zeta^{2}\varepsilon_{xx}-1=0$. These two dispersion relations describe, respectively, purely electrostatic modes (with $\mathbf{\delta E}\parallel\widehat{\mathbf{y}}$) and purely electromagnetic (or inductive) modes (with $\mathbf{\delta E}\parallel\widehat{\mathbf{x}}$). Assuming that $f_{\alpha}^{(0)}(\mathbf{u})$ is even in $p_{y}$, $\varepsilon_{xy}$ reduces to

[TABLE]

This expression is generally nonzero, and hence the CFI excites mixed electromagnetic/electrostatic fluctuations. This feature has been often overlooked in the past, $\delta E_{y}=0$ being assumed from the outset in a number of calculations Molvig (1975); Cary et al. (1981); Okada and Niu (1980); Hill et al. (2005). The electromagnetic ($\delta E_{x}$) and electrostatic ($\delta E_{y}$) components of the solutions to Eq. (28) verify Bret et al. (2004, 2007):

[TABLE]

In Sect. II, $\beta_{\rm p|w}$ was determined by solving $\varepsilon_{xy}=0$ in the fluid limit, exploiting the fact that, to leading order, this equation is independent of the complex frequency $\zeta$. In the general kinetic case, however, $\varepsilon_{xy}$ depends on $\zeta$, whose knowledge involves solving Eq. (28). In practice, we are interested in the fastest-growing mode ($\zeta_{\rm max}$), which should be calculated in the (unknown) “Weibel frame”. Instead of solving simultaneously Eq. (28) and $\varepsilon_{xy|\rm w}=0$ for $\zeta_{\rm max|w}$ and $\beta_{\rm p|w}$, we follow a different approach, noting that the $(\delta E_{y},\delta B_{z})$ fluctuations pertaining to a given mode in the plasma and “Weibel frames” are related through

[TABLE]

since $\delta E_{y|\rm w}=0$. The velocity of the “Weibel frame” relative to the plasma rest frame is therefore given by

[TABLE]

Making use of Eq. $\eqref{eq:EysEx}$ and of $\delta B_{z|\rm p}=-\delta E_{x|\rm p}/\zeta_{\rm max|p}$ , one obtains

[TABLE]

This formula has the advantage of involving only quantities measured in the plasma frame.

We apply the above formalism to the case of Maxwell-Jüttner momentum distribution functions Jüttner (1911):

[TABLE]

where $\beta_{\alpha}\equiv\langle u^{x}/\gamma m_{e}\rangle$ is the normalized mean drift velocity of species $\alpha$ (corresponding drift Lorentz factor $\gamma_{\alpha}$), $\mu_{\alpha}\equiv m_{e}/T_{\alpha}$ and $K_{2}$ denotes a modified Bessel function of the second kind.

Compact expressions of the tensor elements $\varepsilon_{lm}$ can be obtained in terms of one-dimensional integrals over the velocity parallel to the wave vector ($v_{\parallel}$) Wright and Hadley (1975); Bret et al. (2007). These calculations are detailed in Appendix A. In Sec. III.4, such expressions will be used for the numerical resolution of Eq. (35) along with Eq. (28). The parameters of the background plasma and suprathermal beam will be then extracted from PIC simulations of relativistic collisionless shocks. In the remainder of this section, we will derive analytic approximations of $\zeta_{\rm max|p}$ and $\beta_{\rm w|p}$, valid in distinct instability regimes for the plasma and beam particles.

The starting points of these calculations are the alternative expressions (123), (124) and (125) of the dielectric tensor. For instance, $\varepsilon_{xx}$ can be rewritten as (123)

[TABLE]

where

[TABLE]

and $\chi_{\alpha}=\gamma_{\alpha}\zeta/\sqrt{1-\zeta^{2}}$. The integrals involved in $\varepsilon_{yy}$ and $\varepsilon_{xy}$ can be put in a similar form, see (127) and (128). Introducing $\Delta s_{\alpha}$ the characteristic width of the integrand in Eq. (38) [except for the denominator $(\chi_{\alpha}-s)^{-1}$], two limiting cases can be considered for each plasma species:

• •

the ‘hydrodynamic’ limit: $|\chi_{\alpha}|\gg\Delta s_{\alpha}$;

• •

the ‘kinetic’ limit: $|\chi_{\alpha}|\ll\Delta s_{\alpha}$.

The dimensionless variable $s$ is introduced immediately before Eq. (123); it corresponds to $\gamma_{\alpha}\beta_{\parallel}\gamma_{\parallel}$, with $\beta_{\parallel}$ the component of the particle velocity along the wavevector. Since the wavevector is transverse to the streaming direction, the typical extent of $s$ in the above integral is, up to the resonant factor, controlled by the proper temperature $1/\mu_{\alpha}$; the parameter $\chi_{\alpha}$ itself corresponds to the apparent phase four-velocity of the mode. Therefore, the meaning of the hydrodynamic limit is that the apparent typical transverse momentum (normalized to $m_{e}$) exceeds the apparent phase four-velocity, while the kinetic limit corresponds to the opposite case. In the following, approximate formulas of the dielectric tensor will be derived in these two limits.

III.1 Evaluation of the dielectric tensor for the background plasma

III.1.1 Hydrodynamic limit

In the outermost part of the precursor, the background plasma is characterized by a nonrelativistic proper temperature, $\mu_{\rm p}\gg 1$ Lemoine et al. (2019b). In this limit, $\Delta s_{\rm p}\simeq\sqrt{2/\mu_{\rm p}}$, and hence the hydrodynamic response of the background plasma implies $\sqrt{\mu_{\rm p}/2}|\chi_{\rm p}|\gg 1$. This condition coincides with the large-argument limit ($\tilde{\chi}_{\rm p}\simeq|\chi_{\rm p}|\sqrt{\mu_{\rm p}/2}\gg 1$) of the $\mathcal{Z}$ and $\mathcal{Z}^{\prime}$ functions involved in the low-temperature expressions derived in Appendix B. These formulas can be further expanded to first order in $1/\mu_{\rm p}$ by making use of the asymptotic series $\mathcal{Z}(\eta)\simeq-1/\eta-1/2\eta^{3}-3/4\eta^{5}\cdots$ Huba (2013):

[TABLE]

In the rest frame of the background plasma, $\beta_{\rm p|p}=0$; further assuming the weak-growth limit, $|\zeta^{2}|\ll 1$, the above relations simplify to

[TABLE]

III.1.2 Kinetic limit

We now consider the limit $\tilde{\chi}_{\rm p}\ll 1$ of Eqs. (137), (139) and (141). Using the power series $\mathcal{Z}(\eta)\simeq i\sqrt{\pi}\exp(-\eta^{2})-2\eta\cdots$ Huba (2013) and assuming $|\zeta^{2}|\ll 1$, this yields, in the background plasma rest frame

[TABLE]

III.2 Evaluation of the dielectric tensor for the suprathermal particles

III.2.1 Hydrodynamic limit

In contrast to the background plasma, the beam particles are characterized by an ultrarelativistic drift velocity in the background plasma rest frame ($\gamma_{\rm b}\equiv\gamma_{\rm b|p}\gg 1$) and a relativistic proper temperature ($\mu_{\rm b}\ll 1$). As a result, the integrand of Eq. (38) presents the approximate width $\Delta s_{\rm b}\simeq 1$, so that the hydrodynamic response of the suprathermal particles requires $|\chi_{\rm b}|\gg 1$. The corresponding dielectric tensor, $\varepsilon_{lm,\rm b}$, is obtained by expanding $(\chi_{\rm b}-s)^{-1}\simeq\chi_{\rm b}^{-1}\left[1+(s/\chi_{\rm b})^{2}\right]$ in Eqs. (126)-(128), and evaluating the various resulting integrals. For $\varepsilon_{xx,\rm b}$, this gives

[TABLE]

where the functions $I(t,\lambda,b)$ and $J(t,\lambda,b)$ are defined by Eqs. (152) and (145), respectively. Working out the derivatives, we find

[TABLE]

Inserting this expression into Eq. (37) yields, to leading order,

[TABLE]

assuming $|\zeta^{2}|\ll 1$ and $\gamma_{\rm b}^{-2}\ll\mu_{\rm b}\ll 1$, and expanding the Bessel functions accordingly.

Applying the same procedure to Eqs. (127) and (128) leads to the hydrodynamic expressions

[TABLE]

III.2.2 Kinetic limit

The kinetic response of the beam particles can be readily obtained, to leading order in $|\chi_{\rm b}|$, from the expansion $(\chi_{\rm b}-s)^{-1}\simeq-i\pi\delta(s)-P(1/s)$ in Eq. (123), where $\delta(s)$ is the Dirac delta function and $P$ denotes the Cauchy principal value, which here vanishes. In general, however, the beam particles appear to be only marginally kinetic in PIC shock simulations, so it could be useful to go to the next order. The series expansions derived in Appendix C are well suited to this purpose. In the limit $|\chi_{\rm b}|\ll 1$, Eqs. (149), (155) and (158) reduce to

[TABLE]

Combining these approximate expressions with Eqs. (123)-(125) and expanding the Bessel functions in the small-argument limit gives

[TABLE]

where we have further assumed $\gamma_{\rm b}\ll 1$.

III.3 CFI growth rates and frame velocities in various plasma response regimes

The previous formulas may now be applied to the case of the precursor of a relativistic shock to derive the growth rate of the purely transverse CFI and the velocity of the “Weibel frame” $\mathcal{R}_{\rm w}$. We consider the various limits in which the plasma and/or the beam can be described in a fluid-like or kinetic approximation, keeping in mind that the most relevant limit for the precursor is that of both kinetic beam and background plasma.

For reference, let us recall that in the limit $|\zeta|=|\omega/k|\ll 1$, which is applicable here, the plasma can be described in the hydrodynamic regime iff $\gamma_{\rm p}|\zeta|\ll\sqrt{2/\mu_{\rm p}}$ (with $\mu_{\rm p}=m/T_{\rm p}$ assumed greater than unity). As for the beam, it can be described in the hydrodynamic regime iff $\gamma_{\rm b}|\zeta|\ll 1$.

In the following, we solve for the dispersion relation in the background plasma rest frame, in order to derive $\beta_{\rm w|p}$ according to Eq. (35). We also assume that the plasma frame is close to the “Weibel frame”, so that the off-diagonal term $\epsilon_{xy}^{2}$ can be neglected in the dispersion relation as written in the background plasma frame. The dispersion relation may then be approximated as

[TABLE]

Finally, in order to make contact with our previous notations, we will repeatedly use the substitution: $\omega_{\rm pb}^{2}\mu_{\rm b}/\omega_{\rm p}^{2}=\xi_{\rm b}(n_{\infty}/n_{\rm p})/\kappa_{T_{\rm b}}^{2}$. This notably implies $\omega_{\rm pb}^{2}\mu_{\rm b}/\omega_{\rm p}^{2}\ll 1$. We also have $\gamma_{\rm b|p}\simeq\gamma_{\rm p|s}$ and $\beta_{\rm b|p}\simeq 1$.

III.3.1 Hydrodynamic plasma and beam

In the hydrodynamic regime (and in the background plasma rest frame), $|\zeta|\ll\mu_{\rm p}^{-1/2}$ and $|\zeta|\ll 1/\gamma_{\rm b|p}$. Hence, the dispersion relation gives to leading order:

[TABLE]

and so the growth rate saturates at $\Gamma_{\rm max}\simeq\omega_{\rm pb}^{2}\mu_{\rm b}$ for $k\gg\sqrt{2}\omega_{\rm p}$.

Adding up the hydrodynamic plasma and beam contributions into Eq. (35) and retaining only leading order terms yields the “Weibel frame” velocity

[TABLE]

As $|\zeta|\ll\omega_{\rm p}/k$ according to Eq. (60), the expression for $\beta_{\rm w|p}$ boils down to

[TABLE]

As $\beta_{\rm b|p}\simeq\beta_{\rm b|w}\simeq 1$, we recover the formula derived within a fluid approach, Eq. (20), provided one sets in the latter the adiabatic index $\widehat{\Gamma}_{\rm b}=2$. This is the value expected for a gas with one degree of freedom in the relativistic limit; the reduced effective dimensionality for the beam response results from the assumption of a purely 1D transverse fluctuation spectrum.

Finally, the fully hydrodynamic regime holds as long as $\min(\sqrt{\mu_{\rm p}/2},\gamma_{\rm b|p})|\zeta_{\rm max}|\gg 1$. Now, expressing Eq. (60) at $k_{\rm max}\equiv\sqrt{2}\omega_{\rm p}$ gives

[TABLE]

so that another way of expressing the validity of the hydrodynamic regime is:

[TABLE]

III.3.2 Kinetic plasma and hydrodynamic beam

In this limit, $1/\gamma_{\rm b|p}\gg|\zeta|\gg\mu_{\rm p}^{-1/2}$; to leading order, the dispersion relation takes the form

[TABLE]

The growth rate reaches its maximum value $\Gamma_{\rm max}=\sqrt{\omega_{\rm pb}^{2}\mu_{\rm b}}$ for $k\lesssim k_{\rm max}\simeq(2\pi\omega_{\rm pb}^{2}\mu_{\rm p}\mu_{\rm b})^{1/6}\omega_{\rm p}^{2/3}$. We define $|\zeta_{\rm max}|=\Gamma_{\rm max}/k_{\rm max}\simeq\left(\omega_{\rm pb}^{2}\mu_{\rm b}/\sqrt{2\pi\mu_{\rm p}}\omega_{\rm p}^{2}\right)^{1/3}$.

Making use of Eqs. (46), (47), (51), and (52), we express Eq. (35) as

[TABLE]

which gives

[TABLE]

III.3.3 Hydrodynamic plasma and kinetic beam

In this limit, we now have $1/\gamma_{\rm b|p}\ll|\zeta|\ll\mu_{\rm p}^{-1/2}$. Using Eqs. (42) and (56), the dispersion relation reduces to

[TABLE]

Let us first assume that $\omega_{\rm p}^{2}/(\mu_{\rm p}\zeta^{2})$ can be neglected in front of $2\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{3}\zeta$. An unstable solution then exists provided $\gamma_{\rm b|p}^{2}>\omega_{\rm p}^{2}/\left(\omega_{\rm pb}^{2}\mu_{\rm b}\right)\sim 1/\xi_{\rm b}$. This is the same condition as encountered in the $4-$fluid system of Sec. II, see Eq. (27) where $c_{\rm eff,b}^{2}\sim 1/\gamma_{\rm b|p}^{2}$ and $\Omega_{\rm pb}^{2}\sim\omega_{\rm pb}^{2}\mu_{\rm b}$.

If this condition is fulfilled, the fastest-growing solution corresponds to $\Gamma_{\rm max}\simeq\frac{8}{9\pi}\sqrt{\frac{2}{3}}\frac{[\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{2}-\omega_{\rm p}^{2}]^{3/2}}{\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{3}}$ and $\zeta_{\rm max}\simeq\frac{8i}{9\pi\gamma_{\rm b|p}}\left(1-\frac{\omega_{\rm p}^{2}}{\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{2}}\right)$. Note that the conditions for a hydrodynamic plasma ($\sqrt{\frac{\mu_{\rm p}}{2}}|\zeta_{\rm max}|>1$) and a kinetic beam ($\gamma_{\rm b|p}|\zeta_{\rm max}|<1$) are then verified, albeit marginally.

Moreover, combining Eqs. (43), (44), (57), and (58) gives the “Weibel frame” velocity

[TABLE]

Inserting the above expressions of $\Gamma_{\rm max}$ and $\zeta_{\rm max}$, there follows

[TABLE]

In the opposite limit, in which $\omega_{\rm p}^{2}/(\mu_{\rm p}\zeta^{2})$ dominates over $2\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{3}\zeta$ in Eq. (68), the dominant mode satisfies $\Gamma_{\rm max}\simeq\sqrt{2\omega_{\rm p}^{2}/\mu_{\rm p}}$ and $\zeta_{\rm max}\simeq i[2\mu_{\rm p}(1-\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{2}/\omega_{\rm p}^{2})]^{-1/2}$ for $k_{\rm max}\equiv\sqrt{2\omega_{\rm p}^{2}-2\omega_{\rm pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{2}}$. To leading order, we thus derive

[TABLE]

III.3.4 Kinetic plasma and beam

Finally, we consider the case of a fully kinetic beam-plasma system. This regime is of particular importance since it is found to hold in most of the precursor region in long-time shock simulations (see Sec. III.4). Using the expressions (45) and (56), the dispersion relation writes

[TABLE]

The dominant CFI mode is then defined by

[TABLE]

and $k_{\rm max}\simeq\sqrt{\frac{2}{3}\mu_{\rm b}}\omega_{\rm pb}\gamma_{\rm b|p}$.

The corresponding expression for the “Weibel frame” velocity is obtained by combining Eqs. (46), (47), (57), and (58). After some algebra, one finds

[TABLE]

or, to leading order and in terms of our usual parameters,

[TABLE]

III.4 Comparison to PIC simulations

In this Section, we compare the relative velocity $\beta_{\rm w|p}$ estimated using the kinetic model of the CFI developed in Sec. III.3 with that extracted from 2D3V PIC simulations.

These simulations have been performed using the massively parallel, relativistic PIC code calder Lefebvre et al. (2003). The shock is generated by means of the standard mirror technique Spitkovsky (2008). The background pair plasma is continuously injected into the domain from the right-hand boundary, and is made to reflect specularly off the left-hand boundary ($x_{|\rm d}=0$). Electrons and positrons are injected with a Maxwell-Jüttner momentum distribution [Eq. (36)] of proper temperature $T_{\infty}=0.01m_{e}$ and mean drift velocity $\beta_{\infty|\rm d}=-(1-1/\gamma_{\infty|\rm d}^{2})^{1/2}$ in the simulation frame. As mentioned earlier, this simulation reference frame coincides with the downstream plasma rest frame (as a consequence of the use of the mirror technique). Our two reference PIC simulations use $\gamma_{\infty|\rm d}=10$ and $\gamma_{\infty|\rm d}=100$, which respectively correspond to shock Lorentz factors (with respect to the upstream) of $\gamma_{\infty}=17$ and $\gamma_{\infty}=173$.

As the simulation proceeds, the right-hand boundary (injector) is progressively displaced towards the right so as to keep the reflected ballistic particles inside the (increasingly large) domain, while speeding up the calculation at early times Spitkovsky (2008). In order to quench the numerical Cherenkov instability, which notoriously hampers simulations of relativistically drifting plasmas, we make use of the Godfrey-Vay filtering scheme, combined with the Cole-Karkkainnen finite difference field solver Godfrey and Vay (2014). We use a mesh size $\Delta x=\Delta y=0.1\,c/\omega_{\rm p}$ and a time step $\Delta t=0.99\Delta x/c$. Periodic boundary conditions are employed in the tranverse direction for both particles and fields. The initial domain size is $L_{x}\times L_{y}=2700\times 340\,(c/\omega_{\rm p})^{2}$ for $\gamma_{\infty|\rm d}=10$ and $L_{x}\times L_{y}=2000\times 200\,(c/\omega_{\rm p})^{2}$ for $\gamma_{\infty|\rm d}=100$. Each cell is initially filled with 10 macro-particles per species (electrons or positrons). The simulation is run until $t_{\rm max}=3600\,\omega_{\rm p}^{-1}$ (resp. $t_{\rm max}=6900\,\omega_{\rm p}^{-1}$) for our simulation with $\gamma_{\infty|\rm d}=10$ (resp. $\gamma_{\infty|\rm d}=100$).

In order to test our model of the CFI developing in the precursor through the interpenetration of suprathermal and background plasma populations, we need to carefully distinguish these two in the simulations. In order to do so, we track the particles according to the sign of their $x-$momentum and how many times this sign has changed, due only to interactions with the electromagnetic turbulence. We then define the background plasma as those particles that move toward negative values of $x$ and that do not have undergone turnarounds, i.e., any change of sign of $u^{x}$. We define the beam particles as those that move towards positive values of $x$, independently of their number of turnarounds. This definition leaves a minority of particles: those that move with $u^{x}<0$ and have undergone at least one turnaround. In this population, however, it becomes difficult to distinguish particles that originate from the right boundary of the simulation box from those that originate from the left boundary; these populations have different temperatures, so that a single-fluid description of this compound population would introduce errors.

Downstream of the shock, the left-moving and right-moving particle populations are identical, because the plasma has isotropized in this simulation frame. There, our definition of suprathermal particles only counts half of the particles, therefore our $\xi_{\rm b}\simeq 1/4$ in this region: in 2D3V simulations, the downstream pressure represents $1/2$ of the energy density, which itself amounts to the incoming energy flux into the shock rest frame $F_{\infty}$.

We extract hydrodynamic moments $n_{\alpha}$, $T_{\alpha}$, $u_{\alpha}$ for each of the beam and background plasma population, assuming that they obey Maxwell-Jüttner momentum distributions. The spatial profiles of these various hydrodynamic quantities have already been presented in Fig. 1, as extracted from the simulations with $\gamma_{\infty|\rm d}=10$ and $\gamma_{\infty|\rm d}=100$ at respective times $t\simeq 3600\,\omega_{\rm p}^{-1}$ and $t\simeq 6900\,\omega_{\rm p}^{-1}$. One can see that $n_{\rm p}$ and $\gamma_{\rm p|\rm d}$ vary weakly across the precursor region, except near the shock front where the incoming plasma slows down significantly and experiences compression. By contrast, the plasma temperature steadily increases from its far-upstream value ($T_{\rm p}=0.01m_{e}$) to unity and beyond when approaching the shock front. This heating results from the interaction with the beam particles, whose density rises by $\sim 4-5$ orders of magnitude across the precursor Lemoine et al. (2019b). The beam Lorentz factor in the simulation frame is close to unity, confirming that the beam drifts at a weakly relativistic velocity in the shock frame.

The general dispersion relation (28) is numerically solved using the reduced forms (109), (115) and (122) of the dielectric tensor elements, as detailed in Ref. Bret et al. (2010b). At each sampled location, we look for the growth rate ($\Gamma_{\rm max}$), wave number ($k_{\rm max}$) and wave phase velocity ($\xi_{\rm max}$) of the fastest-growing mode, and then use Eq. (35) to evaluate the corresponding value of the “Weibel frame” velocity ($\beta_{\rm w|p}$).

Figures 3 display the spatial profiles of the $\chi_{\rm b}$ and $\tilde{\chi}_{\rm p}$ parameters defined by Eqs. (129) and (134), respectively. Interestingly, both the background plasma and beam populations appear to lie in the kinetic CFI regime ($\tilde{\chi}_{\rm p}<1$ and $\chi_{\rm b}<1$) throughout the precursor region. However, whereas $\tilde{\chi}_{\rm p}$ shows weak variations around relatively low values ($\tilde{\chi}_{\rm p}\sim 0.05$), so that taking the kinetic plasma limit is well justified, the $\chi_{\rm b}$ values are larger by about an order of magnitude and present stronger variations. Where $\xi_{\rm b}$ becomes of the order of unity, the beam response is then only marginally kinetic; consequently, analytical approximations present an error of about a factor 2 with respect to the full numerical calculation of $\beta_{\rm w|p}$.

As discussed in Sec. I.2, we estimate the 3-velocity of the “Weibel frame” in PIC simulations as the ratio $\beta_{\rm w|d}=\left\langle\delta E_{y}^{2}\right\rangle^{1/2}/\left\langle\delta B_{z}^{2}\right\rangle^{1/2}$, where averaging is performed over the transverse dimension. Given the locally measured value of the background plasma velocity $\beta_{\rm p|d}$, we convert it to the instantaneous local plasma rest frame through the standard transform $\beta_{\rm w|p}=\left(\beta_{\rm w|d}-\beta_{\rm p|d}\right)/\left(1-\beta_{\rm w|d}\beta_{\rm p|d}\right)$. We then compare in Fig. 4 this measurement with our theoretical estimates of the relative velocity, $\beta_{\rm w|p}$, between the “Weibel frame” and the background plasma, using both the numerical calculation (red circle/dashed line) and the analytical approximation (purple square/dashed). Both for $\gamma_{\infty}=17$ and $\gamma_{\infty}=173$, it is seen that, as predicted, $\beta_{\rm w|p}$ remains subrelativistic throughout the precursor, increasing from $\beta_{\rm w|p}\sim 10^{-3}$ at the far end of the precursor up to $\beta_{\rm w|p}\simeq 0.1$ near the shock front. The theoretical estimates appear to provide reasonable match to the PIC data in the region where $\beta_{\rm w|p}$ can be measured accurately.

For both reference simulations, Fig. 4 reveals significant fluctuations in the measured values of $\beta_{\rm w|p}$, with increasing amplitude at large $x$, which deserves some discussion. Far from the shock front, the magnitude of our observable $\left\langle\delta E_{y}^{2}\right\rangle^{1/2}/\left\langle\delta B_{z}^{2}\right\rangle^{1/2}$ is close to unity (see e.g., Fig. 2), which implies that any amount of numerical noise can artificially push $|\beta_{\rm w|d}|$ to values larger than unity, even though its true value might be $<1$. Furthermore, when transforming values to the background plasma rest frame, any error in $\beta_{\rm w|d}$ is amplified by $\sim\gamma_{\rm p|d}^{2}$, which is large far from the shock.

In a given portion of the precursor, the “Weibel frame” can be considered well defined where $|\beta_{\rm w|d}|<1$ for most data points. In Fig. 4, values $|\beta_{\rm w|p}|<1$ correspond to values $|\beta_{\rm w|d}|<1$, but the use of raw numerical data, binned linearly, and plotted on a log-scale, precludes a clear identification of the region where this frame is well-defined, at least by eye. By rebinning the data, however, we infer that the “Weibel frame” is well defined in the range $x\lesssim 10^{3}\,c/\omega_{\rm p}$, for both PIC simulations, as claimed in Sec. I.2. At larger distances from the shock, $|\beta_{\rm w|p}|$ fluctuates too often on both sides of unity, either because of numerical noise, or because of the contribution of electrostatic modes. Given the magnitude of the fluctuations, we have chosen to plot the data in full color only in the region $x\leq 300\,c/\omega_{\rm p}$, for the sake of being conservative.

IV The “Weibel frame” in the nonlinear regime

IV.1 Theoretical model

We adopt here another perspective on the “Weibel frame”. Specifically, we consider the nonlinear evolution of the CFI, once the current filaments have formed, borrowing on the work of Ref. Vanthieghem et al. (2018). We assume that, at each point in the shock precursor, the CFI has developed a quasistatic, transversally periodic system of current filaments. By quasistatic, it is meant that an equilibrium approximately holds between magnetic and thermal pressures in the filaments, according to the physical conditions at the point considered, and that these conditions evolve slowly enough that this equilibrium has time to adapt from one point to another. We also consider a 2D configuration, without loss of generality, so that the plasma is periodic along the $y$-axis and the species drift along the $x$-axis; the magnetic field is then transverse to the $(x,y)$-plane. The dominant components of the four-potential now explicitly depend on $y$ and are given by $A^{\mu}=\left[\phi(y),A_{x}(y),0,0\right]$. We assume a drifting Maxwell-Jüttner distribution for each of the four species (corresponding to two counterstreaming pair distributions), so that

[TABLE]

where $q_{\alpha}$ denotes the charge of the species and $n_{\alpha 0}$ represents a normalization prefactor. We can inject these density profiles in the potential formulation of the Ampère-Maxwell and Gauss-Maxwell equations, leading to

[TABLE]

where $\omega_{\alpha}=\sqrt{4\pi n_{\alpha 0}e^{2}/m_{e}}$ scales to the plasma frequency of species $\alpha$. It is worth noting that this system can also be obtained in a 4-fluid framework with isothermal closure condition Vanthieghem et al. (2018).

Following Vanthieghem et al. (2018), we introduce the plasma nonlinearity parameter:

[TABLE]

In the weakly nonlinear limit, $\Xi_{\rm p}\ll 1$, the sinh functions in the above equations can be approximated to unity, so that the vanishing of the electrostatic component entails:

[TABLE]

where the subscript w has been introduced, because the above quantities are now defined in the “Weibel frame”, $\mathcal{R}_{\rm w}$, in which $\phi_{\rm|w}=0$.

In the weakly nonlinear limit, the velocity of $\mathcal{R}_{\rm w}$ can be computed exactly using relation (83). Writing $\beta_{\rm b|w}\,=\,(\beta_{\rm b}-\beta_{\rm w})/(1-\beta_{\rm b}\beta_{\rm w})$, $\gamma_{\rm b|w}\,=\,\gamma_{\rm b}\gamma_{\rm w}(1-\beta_{\rm b}\beta_{\rm w})$ etc., in any given frame, one finds that $\beta_{\rm w}$ is solution to the equation:

[TABLE]

with

[TABLE]

Writing $n_{\rm b}$ in terms of $\xi_{\rm b}$ as before, we expand the above solution to first order in $\xi_{\rm b}$ to obtain the relative velocity $\beta_{\rm w|p}$:

[TABLE]

where $\gamma_{\rm b|p}\sim\gamma_{\rm p}$ represents the relative Lorentz factor between the beam and the background plasma. Interestingly, Eq. (86) corresponds to our earlier expression Eq. (78) obtained from the linear dispersion relation of the CFI in the kinetic plasma – kinetic beam limit.

In Paper II Lemoine et al. (2019b), it is argued that $\gamma_{\rm p}^{2}\xi_{\rm b}$ is much smaller than unity in the far precursor, where the incoming background plasma maintains its initial velocity, i.e., $\gamma_{\rm p}\simeq\gamma_{\infty}$, and of the order of unity in the near precursor, where $\gamma_{\rm p}<\gamma_{\infty}$ due to deceleration. The above thus implies that the “Weibel frame”, in this nonlinear description, moves at subrelativistic velocities relative to the background plasma, as in the linear limit studied earlier. That $\beta_{\rm w|p}$ is positive means that the “Weibel frame” moves at slightly smaller absolute velocity towards the shock front than the background plasma.

In the near precursor, as the background plasma is heated up to relativistic temperatures, $\beta_{\rm w|p}$ increases in magnitude; this implies that the background plasma decouples from the “Weibel frame”, hence increasing the heating rate and leading to the shock transition.

One can also compute the first nonlinear correction in $\Xi_{\rm p}$ to the above velocity. To this effect, we recast Eq. (81) in terms of the nonlinearity parameters of the beam ($\Xi_{\rm b}$) and the plasma ($\Xi_{\rm p}$) in relation (83):

[TABLE]

As discussed in Paper III Lemoine et al. (2019c), the beam particles carry such inertia that they hardly participate in the filamentation, meaning $\Xi_{\rm b|w}\ll 1$. In $\mathcal{R}_{\rm w}$, $\Xi_{\rm b|w}$ indeed represents the ratio of the electromagnetic component $eA_{x|\rm w}$ to the typical momentum $T_{\rm b|w}=T_{\rm b}/\gamma_{\rm b|w}$ of the particles, and for suprathermal particles, this ratio is much smaller than unity. Assuming $\Xi_{\rm p|w}\lesssim 1$, we then obtain to lowest order

[TABLE]

In this configuration, the nonlinearity appears as a second-order correction to the “Weibel frame” speed. In detail, one obtains

[TABLE]

Note that the “Weibel frame” is not always well-defined in the strongly nonlinear limit: if $\Xi_{\rm p,b}\gg 1$, Eq. (81) decouples into a set of two relations between the temperatures and densities

[TABLE]

which overdetermine the system for a given set of parameters. From Eq. (81), we see that the error on the electrostatic fields at leading order in $\Xi_{\rm p|w}$ if we evaluate the “Weibel frame” from relation (83) evolves as $\delta\phi\propto\Xi_{\rm p|w}^{3}$. In the case of interest, however, PIC simulations indicate that the weakly nonlinear limit represents a good approximation in the precursor of relativistic shocks, see Sec. IV.2.

IV.2 Comparison to PIC simulations

In this section, we confront the result of Sec. (IV) with the PIC simulations presented in the Sec. III.4. A strong hypothesis underlying the above formulae is that of a marginally nonlinear plasma response $\Xi_{\rm p|w}<1$. We wish to motivate this hypothesis by measuring this nonlinearity parameter. Assuming Eq. (79) holds, one expects

[TABLE]

We can then estimate the nonlinearity parameter in the simulation through the following relation

[TABLE]

where $\langle\cdot\rangle$ denotes the mean value along the direction transverse to the drift.

Using relation (93), we present in Figs. 5 our estimate of $\Xi_{\rm p|w}$ for the two reference simulations ($\gamma_{\infty}=17$ in the top panel, $\gamma_{\infty}=173$ in the bottom panel). In both cases, the nonlinearity parameter tends to increase from values well below unity in the far precursor, up to near unity within a few hundred $c/\omega_{\rm p}$ to the shock front. This behavior suggests that one should use the present nonlinear equilibrium model to estimate $\beta_{\rm w|p}$ in the near precursor, say $x\lesssim 300c/\omega_{\rm p}$, and the linear estimate obtained in Sec. III at larger distances, where $\Xi_{\rm p|w}$ falls to values small compared to unity. One must, however, keep in mind that the present approach is based on a fluid model, while the former is fully kinetic, and that the analysis of Sec. III indicates that, from the point of view of the instability, both the beam and the background plasma should be treated kinetically. Fortunately, both the present nonlinear estimate of $\beta_{\rm w|p}$, Eq. (86), and its linear kinetic counterpart, Eq. (78), match one another. Hence, our theoretical estimate of $\beta_{\rm w|p}$ can be considered as rather well established throughout the precursor.

This is supported by Figs. 6, which compare the values of $\beta_{\rm w|p}$ extracted from the simulations with the theoretical estimate (86), assuming (weakly) nonlinear equilibrium filaments. As before, the PIC simulation data are light colored where they cannot be measured accurately. In both simulation cases, the theoretical estimates appear in reasonable agreement with the PIC results, especially in the near precursor $x\lesssim 10^{3}c/\omega_{\rm p}$, where the velocity of the “Weibel frame” is well defined.

V Conclusions

This paper belongs to a series of articles in which we build a theoretical model of unmagnetized, relativistic collisionless pair shocks, and compare it with dedicated PIC simulations. More specifically, we have discussed in the present work the physics of the purely transverse CFI that results from the interpenetration of the background plasma and the beam of suprathermal particles in the shock precursor. We have argued that there exists a frame $\mathcal{R}_{\rm w}$, referred to as the “Weibel frame” in which the instability is mostly magnetic in nature. We have derived the velocity of this frame at each point of the shock precursor, through a kinetic description of the linear stage of the instability, as well as through a quasistatic model of the nonlinear phase of the filamentation instability. This “Weibel frame” is of particular importance to the physics of the shock and of the acceleration process, because it represents the frame of the scattering centers.

We emphasize the following properties of this “Weibel frame”: (i) it is found to move at subrelativistic velocities relative to the background plasma, (i.e., with 4-velocity $u_{\rm w|p}<1$), and, therefore, at relativistic velocities toward the shock front. This result can be related to the strong asymmetry of the interation between the beam of suprathermal particles and the background plasma. We also find that, (ii), the 3-velocity $\beta_{\rm w|p}$ is of opposite sign to $\beta_{\rm p}$, the velocity of the background plasma in the shock rest frame, implying that the “Weibel frame” moves slightly less fast than the background plasma relative to the shock front. Furthermore, (iii), the relative velocity $\beta_{\rm w|p}$ typically scales as $\xi_{\rm b}$, which represents the beam pressure normalized to the incoming momentum flux at infinity. Finally, it is to be emphasized that the “Weibel frame” is not globally inertial, because its velocity $\beta_{\rm w}$ depends on the distance to the shock. As discussed in  Lemoine et al. (2019a), and more particularly Lemoine et al. (2019b), which addresses in detail the consequences of this noninertial nature, the “Weibel frame” decelerates from the far to the near precursor. In the present article, we have determined the velocity of this frame at each point in the precursor, given the local physical conditions, in a WKB-like approximation.

Our PIC simulations confirm these various features. In particular, they reveal that $\delta E_{\perp}<\delta B_{\perp}$ close to the shock front and over most of the precursor, and hence that the “Weibel frame” is well defined; relatively to the shock front, the “Weibel frame” velocity, measured through the ratio of $\delta E_{\perp}/\delta B_{\perp}$, is also found to be slightly below the background plasma velocity; finally, the relative velocity $\beta_{\rm w|p}$ generally decreases away from the shock front, like $\xi_{\rm b}$.

Our kinetic model relies on the use of Maxwell-Jüttner distribution functions for the beam and the background plasma, and it provides (at the expense of rather complex calculations) simple approximations to the velocity of the “Weibel frame” and to the growth rate of the purely transverse CFI. Our quasistatic model of the nonlinear stage of the instability describes the current filaments as periodic magnetostatic structures in the “Weibel frame”, in pressure equilibrium with the plasma. Our PIC simulations indicate that, over most of the precursor, those filaments are in a mildly nonlinear stage, with a nonlinearity parameter below unity. We expect that both models should capture the salient features of the instability and indeed, the resulting formulae turn out to bracket rather well the value $\beta_{\rm w|p}$ seen in PIC simulations.

In subsequent papers, it will be shown that the “Weibel frame” plays an essential role in shaping the microphysics of the shock transition, in particular the physics of heating and deceleration of the background plasma.

Acknowledgements.

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

Appendix A Reduced expressions of the kinetic dielectric tensor

The triple integrals involved into Eq. (III) after substitution of the model distribution (36) can be recast in the form of one-dimensional quadratures Bret et al. (2007). To do so, we mostly follow the lines sketched in Ref. Wright and Hadley (1975), where it is shown that two of the momentum integrations can be carried out in closed form.

Let us first consider the dielectric tensor element $\varepsilon_{xx}$, which is involved in the purely electromagnetic dispersion relation of the CFI (i.e., in a frame close to the “Weibel frame”). After straightforward algebra, it can be rewritten as

[TABLE]

where we have defined

[TABLE]

The first integral in the right-hand side of (94) can be exactly solved by noting that

[TABLE]

We then proceed by changing to velocity variables in cylindrical coordinates along the wave vector: $\mathbf{v}=(v_{\parallel},v_{\perp}\cos\theta,v_{\perp}\sin\theta)$. Making use of ${\rm d}^{3}u=\gamma^{5}{\rm d}^{3}v$, Eq. (94) then becomes

[TABLE]

with $\gamma_{\parallel}=(1-v_{\parallel}^{2})^{-1/2}$. Given the integral representation of the modified Bessel functions of the first kind Abramowitz and Stegun (1972),

[TABLE]

we obtain

[TABLE]

Here we have introduced

[TABLE]

From the definition $t=\gamma/\gamma_{\parallel}$ follow the relations $\gamma v_{\perp}=(t^{2}-1)^{1/2}$ and $v_{\perp}{\rm d}v_{\perp}={\rm d}t/\gamma_{\parallel}^{2}t^{3}$, allowing us to express $\mathcal{S}_{30}^{\alpha}$ and $\mathcal{S}_{32}^{\alpha}$ as

[TABLE]

which can be further recast as

[TABLE]

We now take advantage of the following formula Gradshteyn and Rizhik (1980)

[TABLE]

valid for $\Re\nu>-1$ and $\delta>\zeta$. There follows

[TABLE]

where $\Pi_{\alpha}=(\gamma_{\parallel}^{2}-\beta_{\alpha}^{2})^{1/2}$. After evaluation of the first and second-order derivatives, one finds

[TABLE]

Substituting Eq. (108) into Eq. (99) finally gives the following simplified expression for $\varepsilon_{xx}$:

[TABLE]

Let us now consider $\varepsilon_{yy}$, which writes

[TABLE]

After changing to cylindrical velocity coordinates, one obtains

[TABLE]

which, making use of Eq. (98), reduces to

[TABLE]

Making use of Eq. (105), one obtains

[TABLE]

and therefore

[TABLE]

Combining Eqs. (112) and (114) then readily yields

[TABLE]

Finally, let us rewrite the non-diagonal dielectric tensor element $\varepsilon_{xy}$ as

[TABLE]

Changing to cylindrical velocity variables leads to

[TABLE]

Using Eq. (98) gives

[TABLE]

where $\mathcal{S}_{21}^{\alpha}$ can be recast in the form

[TABLE]

Exploiting Eq. (105) yields

[TABLE]

After evaluation of the derivatives, one finds

[TABLE]

Combining Eqs. (118) and (121), one obtains the simplified expression

[TABLE]

It should be stressed that the compact expressions (109), (115) and (122) are strictly valid for $\Im\zeta>0$ only. They thus lend themselves readily to the numerical resolution of the dispersion relation (28) when searching for (purely growing) unstable modes only, as has been done in Refs. Bret et al. (2007, 2008, 2010b). If damped modes are also examined, care must be taken for the analytic extension to the lower complex $\zeta$-plane of the integrals, owing to the presence of branch points at $\zeta=\pm 1$ Mikhailovskii (1981); Bret et al. (2010b).

To obtain analytic approximations, it is convenient to make the change of integration variable $v_{\parallel}\to s$, such that $\gamma_{\parallel}=\sqrt{1+s^{2}/\gamma_{\alpha}^{2}}$. This gives the alternative expressions

[TABLE]

where the $\mathcal{A}_{lm}^{\alpha}$ terms are defined as

[TABLE]

and we have introduced

[TABLE]

Appendix B Low-temperature expression of the kinetic dielectric tensor

In the following we expand the dielectric tensor elements in the nonrelativistic thermal limit $\mu_{\alpha}\gg 1$, of particular relevance to the background plasma particles. Our starting point is Eqs. (123)-(128).

Let us first address $\varepsilon_{xx}$, by rewriting (126) as

[TABLE]

where the integral function $\mathcal{B}_{xx}^{\alpha}(\mu_{\alpha},\chi_{\alpha})$ is of the form

[TABLE]

Applying Laplace’s method to $\mathcal{B}_{xx}^{\alpha}(\mu_{\alpha},\chi_{\alpha})$ Bender and Orszag (1978), one obtains to first order in $1/\mu_{\alpha}$:

[TABLE]

Using $\varphi(0)=1$, $\varphi^{\prime\prime}(0)=1$, $\varphi^{(3)}(0)=0$ and $\varphi^{(4)}(0)=-3$, as well as the expansions $f_{\alpha}(0)\simeq\gamma_{\alpha}^{2}\beta_{\alpha}^{2}+(1+3\gamma_{\alpha}^{2}\beta_{\alpha}^{2})/\mu_{\alpha}$ and $f^{\prime\prime}_{\alpha}(0)\simeq-3\gamma_{\alpha}^{2}\beta_{\alpha}^{2}$, one gets

[TABLE]

where we have defined

[TABLE]

Introducing the well-known plasma dispersion function Huba (2013)

[TABLE]

and noting that $\mathcal{Z}^{\prime}(\eta)=-2\left[1+\eta\mathcal{Z}(\eta)\right]$, Eq. (133) can be conveniently expressed as

[TABLE]

By combining this equation with Eqs. (123) and (130), and using the small-argument expansion $K_{2}(x)\simeq\sqrt{\frac{\pi}{2x}}e^{-x}(1+\frac{15}{8x})$ Abramowitz and Stegun (1972), one obtains the following low-temperature approximation of $\varepsilon_{xx}$:

[TABLE]

Note that $\tilde{\chi}_{\alpha}$ is the correct relativistic equivalent of the standard argument ($\eta=\sqrt{\frac{\mu_{\alpha}}{2}}\zeta$) of the $\mathcal{Z}(\eta)$ function involved in the nonrelativistic CFI dispersion relation Davidson et al. (1972). A similar result was obtained in Ref. Schlickeiser and Kneller (1997) in the case of electrostatic plasma waves.

Likewise, the integral involved in $\varepsilon_{yy}$ (124) can be expanded to first order in $1/\mu_{\alpha}$

[TABLE]

so that

[TABLE]

Finally, we expand the integral involved in $\varepsilon_{xy}$ (125) as

[TABLE]

Combining this expression with Eq. (125) and identifying the $\mathcal{Z}$ and $\mathcal{Z}^{\prime}$ functions yields

[TABLE]

Expressions (137), (139) and (141) provide the sought-for low-temperature expansions of the CFI dielectric tensor. The involved $\mathcal{Z}$ and $\mathcal{Z}^{\prime}$ functions can be readily evaluated in the entire complex plane using, e.g., the fast solver developed in Ref. Weideman (1994).

Appendix C Series expansion of the dielectric tensor

In similar fashion to Ref. Schlickeiser and Kneller (1997), the dielectric tensor can be expanded in the form of an infinite series, which proves convenient for deriving approximations in the kinetic regime.

Let us first address $\varepsilon_{xx}$, remarking that when $\Im\chi_{\alpha}>0$ Eq. (126) can be rewritten as

[TABLE]

where we have exploited the identity

[TABLE]

Equation (142) can be put in the form

[TABLE]

where we have introduced Erdelyi et al. (1954)

[TABLE]

Substituting the above expression into Eq. (142) yields

[TABLE]

Expanding $e^{it\chi_{\alpha}}=\sum_{n=0}^{\infty}(it\chi_{\alpha})^{n}/n!$, and making use of the identities Gradshteyn and Rizhik (1980)

[TABLE]

we obtain after some algebra,

[TABLE]

The above series can be further simplified using the multiplication theorem of Bessel functions Abramowitz and Stegun (1972), giving finally

[TABLE]

Likewise, the integral involved in $\varepsilon_{yy}$ (124) can be recast as

[TABLE]

Here, we have defined

[TABLE]

which can be solved in closed form as Erdelyi et al. (1954)

[TABLE]

As in the previous derivation, we substitute Eqs. (153) into Eq. (151), expand the exponential factor and exploit Gradshteyn and Rizhik (1980)

[TABLE]

There follows

[TABLE]

Applying again the multiplication theorem of Bessel functions, we find the following expression for $\varepsilon_{yy}$

[TABLE]

Reiterating the previous procedure, we first express the integral term involved in $\varepsilon_{xy}$ (125) as

[TABLE]

Evaluating the derivatives of the function $J$ defined by Eq. (145), we obtain

[TABLE]

Finally, Bessel-function identities allow us to obtain the following alternative expression for $\varepsilon_{xy}$

[TABLE]