# Current-driven filamentation upstream of magnetized relativistic collisionless shocks

Martin Lemoine, Guy Pelletier, Laurent Gremillet, Illya Plotnikov

► To cite this version:

Martin Lemoine, Guy Pelletier, Laurent Gremillet, Illya Plotnikov. Current-driven filamentation upstream of magnetized relativistic collisionless shocks. Monthly Notices of the Royal Astronomical Society, 2014, 440, pp.1365-1378. 10.1093/mnras/stu213. insu-03618578

HAL Id: insu-03618578

https://insu.hal.science/insu-03618578

Submitted on 24 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Current-driven filamentation upstream of magnetized relativistic collisionless shocks

Martin Lemoine,<fnref n="1" />★ Guy Pelletier,<fnref n="2" /> Laurent Gremillet<fnref n="3" /> and Illya Plotnikov<fnref n="2" />

Accepted 2014 January 28. Received 2014 January 24; in original form 2013 October 29

## ABSTRACT

The physics of instabilities in the precursor of relativistic collisionless shocks is of great importance in high-energy astrophysics, because these instabilities build up the shock, control the particle acceleration process and generate the magnetic fields in which the accelerated particles radiate. Two crucial parameters control the micro-physics of these shocks: the magnetization of the ambient medium and the Lorentz factor of the shock front; as of today, much of this parameter space remains to be explored. In the present paper, we report on a new instability upstream of electron–positron relativistic shocks and we argue that this instability shapes the micro-physics at moderate magnetization levels and/or large Lorentz factors. This instability is seeded by the electric current carried by the accelerated particles in the shock precursor as they gyrate around the background magnetic field. The compensation current induced in the background plasma leads to an unstable configuration, with the appearance of charge neutral filaments carrying a current of the same polarity, oriented along the perpendicular current. This “current-driven filamentation” instability grows faster than any other instability studied so far upstream of relativistic shocks, with a growth rate comparable to the plasma frequency. Furthermore, the compensation of the current is associated with a slow-down of the ambient plasma as it penetrates the shock precursor (as viewed in the shock rest frame). This slow-down of the plasma implies that the “current- driven filamentation” instability can grow for any value of the shock Lorentz factor, provided the magnetization $\sigma \lesssim 10^{-2}$. We argue that this instability explains the results of recent particle-in-cell simulations in the mildly magnetized regime.

**Key words:** acceleration of particles – instabilities – plasmas – shock waves.

## 1 INTRODUCTION

The physics of particle acceleration at relativistic collisionless shock waves plays a key role in the description of a number of powerful astrophysical objects, e.g. blazars, pulsar wind nebulae, gamma-ray bursts, etc. One of the lessons learned in the past decade in this field of research is the importance of the non-linear relationship that ties the acceleration process and the generation of micro-turbulence in the shock vicinity. It was anticipated early on that the self-generation of micro-turbulence on length scales much smaller than the gyroradius of the accelerated particles is a necessary condition for the proper development of the relativistic Fermi process (Lemoine, Pelletier & Revenu 2006), in agreement with test particle Monte Carlo simulations (Niemiec, Ostrowski & Pohl 2006). This

small-scale nature of the turbulence comes with a number of important consequences, most notably the limited maximal energy of particles accelerated at ultra-relativistic shock waves [see e.g. Kirk & Reville 2010; Bykov et al. 2012; Plotnikov, Pelletier & Lemoine 2013).

The particle-in-cell (PIC) numerical simulations of Spitkovsky (2008a,b) have confirmed the validity of these arguments and offered a more exhaustive picture of the acceleration process in the ultra-relativistic unmagnetized limit. These simulations have shown that the accelerated (supra-thermal) particle population excites filamentation instabilities upstream of unmagnetized shock waves (meaning, shock waves propagating in an unmagnetized medium; see also Nishikawa et al. 2009); these instabilities build up a magnetic barrier on plasma scales $c/\omega_p$ and at the same time serve as scattering centres for the acceleration process. As the magnetic field energy density grows to an equipartition fraction $\epsilon_B \sim 10^{-1}$ ($\epsilon_B$ denotes the fraction of incoming kinetic energy flux in the shock front rest frame stored in magnetic energy), incoming particles can be

<footnote n="1">¹Institut d'Astrophysique de Paris, CNRS, UPMC, 98 bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="2">²UJF-Grenoble I/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, F-38041 Grenoble, France</footnote>

<footnote n="3">³CEA, DAM, DIF, F-91297 Arpajon, France</footnote>

<footnote>* E-mail: lemoine@iap.fr</footnote>
isotropized on a coherence length scale of the order of ~10c/ωp, thereby initiating the shock transition. The gyroradius of accelerated particles remains larger than this length scale and the Fermi acceleration process develops as anticipated. These simulations have been confirmed, and followed by further PIC simulations with different conditions, in particular regarding the degree of magnetization of the upstream (background) plasma, the obliquity of the magnetic field and the nature (pairs versus electron-proton) of the incoming flow (e.g. Keshet et al. 2009; Martins et al. 2009; Sironi & Spitkovsky 2009, 2011; Haugbølle 2011; Sironi, Spitkovsky & Arons 2013).

The physics of the electromagnetic instabilities that lead to the formation of a ultra-relativistic collisionless shock and to the self-sustenance of the shock have naturally received a lot of attention: see e.g. Hoshino & Arons (1991), Hoshino et al. (1992) and Gallant et al. (1992) for magnetized shock waves; for weakly magnetized shock waves, see e.g. Medvedev & Loeb (1999), Wiersma & Achterberg (2004), Lyubarsky & Eichler (2006), Milosavljević & Nakar (2006), Achterberg & Wiersma (2007), Achterberg, Wiersma & Norman (2007), Pelletier, Lemoine & Marcowith (2009), Lemoine & Pelletier (2010, 2011), Bret, Gremillet & Bénisti (2010), Rabinak, Katz & Waxman (2011) and Shaisultanov, Lyubarsky & Eichler (2012). To summarize in a few lines the current understanding, the Weibel/filamentation instability appears to play a leading role in the generation of the small-scale magnetic field in the weakly magnetized shock limit, although electrostatic oblique modes and Buneman modes retain their importance in preheating the electrons away from the shock front (see the discussion in Lemoine & Pelletier 2011). At strongly magnetized shock waves, the synchrotron maser instability is recognized as the leading agent of dissipation (see e.g. Hoshino & Arons 1991; Hoshino et al. 1992; Gallant et al. 1992).

However, at intermediate magnetizations and/or very large Lorentz factors, the physics remains poorly known. Indeed, the filamentation instability and other two stream modes cannot be excited in these regions of parameter space, because the time-scale on which the incoming particles cross the precursor becomes shorter than the time-scale on which such instabilities can be excited (Lemoine & Pelletier 2010, 2011). Therefore, how the shock is structured in such conditions remains an open question.

Following up on Lemoine et al. (2014), we report here on a new current-driven instability which is likely to emerge as the dominant instability in this range of magnetization and at very large Lorentz factors. The electric current is carried by the suprathelial particles (or shock reflected particles) and results from their gyration in the background magnetic field: assuming that the magnetic field is oriented along the z axis, while the incoming plasma flows along -x in the shock rest frame, the current is generated along -y, since the Lorentz force deflects positive and negative suprathelial particles in opposite directions. As the ambient plasma penetrates the precursor, it develops a compensating current along +y. This configuration is found to be unstable, because a current fluctuation can couple to a density fluctuation and excite a combination of extraordinary modes and compressive modes of the ambient plasma. This will be made explicit further on.

As viewed from the rest frame of the ambient plasma, this perpendicular electric current is extraordinarily large. If one writes $\xi_{cr}$ the fraction of incoming kinetic energy flux carried by the suprathelial particles – see equation (1) – with $\xi_{cr} \sim 0.1$ indicated by PIC simulations, $\gamma_{sh} \gg 1$ the Lorentz factor of the shock wave in the ambient plasma frame and $n_u$ the proper density of the ambient plasma, the induced current reads $j_{y,cr} \sim \gamma_{sh}\xi_{cr}n_u ec$. For $\gamma_{sh}\xi_{cr} \gtrsim 1$, as expected in ultra-relativistic shocks, this current cannot be

compensated by the ambient plasma at rest. As we will demonstrate, the latter is actually accelerated to relativistic velocities relatively to its initial rest frame and it is squashed to an apparent density ~$\gamma_{sh}\xi_{cr}n_u$ in the frame in which there is no bulk motion along x (denoted as $\mathcal{R}$ in the following); then, particle motion at relativistic velocities along y leads to current compensation.

In this work, we focus on an electron-positron shock; in electron-ion shocks, a similar current develops but excites other modes, in particular Whistler waves. This case will be discussed in a forthcoming paper. In Section 2, we discuss the physics of the instability at the linear level, using a relativistic two-fluid model for the incoming background plasma exposed to a rigid external current set by the suprathelial particles. In Section 3, we discuss the relevance of this instability in relativistic collisionless shocks and compare it to results of recent PIC simulations. We discuss the structure of the precursor in Appendix A and provide conclusions in Section 4.

## 2 CURRENT-DRIVEN FILAMENTATION INSTABILITY

We describe the shock precursor as follows, in the shock front frame. The incoming plasma flows with four-velocity $u_x < 0$, carrying magnetic field $B = B_z z$ and convective electric field $E = \gamma_{sh}\beta_{sh}B_u y$, with $\beta_{sh} < 0$ the velocity of the incoming background plasma in the shock rest frame in units of c, i.e. $\gamma_{sh} \equiv (1 - \beta_{sh}^2)^{-1/2}$. In principle, $B_z$ depends on x, while $B_u$ corresponds to the upstream magnetic field measured in the upstream rest frame well beyond the precursor. The precursor also contains a population of relativistic suprathelial particles, which rotate around B and thereby induce a current along y, $j_{cr} \sim -\gamma_{sh}\xi_{cr}n_u ec y$. The quantity $\xi_{cr}$ characterizes the fraction of the incoming particle energy carried by the suprathelial particles:

$$ \xi_{cr} \equiv \frac{e_{cr}}{\gamma_{sh}^2 n_u m c^2}, \qquad (1) $$

with $e_{cr} = n_{cr}\gamma_{sh}mc^2$ in the shock frame, assuming that the suprathelial particles carry a density $n_{cr}$ and typical Lorentz factor $\gamma_{sh}$; from equation (1), one derives $n_{cr} = \gamma_{sh}\xi_{cr}n_u$, whence the expression for the current density $j_{cr}$.

The spatial profile of this current and the overall structure of the precursor are described in detail in Appendix A; Fig. 1 offers a sketch of the precursor. The typical size of the precursor is $c/\omega_c$, with $\omega_c = eB_u/(mc)$ the upstream cyclotron frequency; this size also corresponds to the typical gyration radius $r_L$ of the suprathelial particles in the shock front rest frame, whose typical Lorentz factor $\sim \gamma_{sh}$.

As the incoming particles cross the precursor, they are deflected along y in order to compensate the cosmic ray perpendicular current. Positrons drift towards +y while electrons drift towards -y. The absolute value of the four-velocity y-component for both fluids is equal, $|u_y| \sim \gamma_{sh}\xi_{cr}$ (in units of c), hence $|u_y| \gtrsim 1$ is expected for relativistic shocks, possibly $|u_y| \gg 1$.

The deflection of the incoming flow along y implies a substantial deceleration of the flow along x, which has drastic consequences regarding the development of the instability. The profile of the velocity of the flow is discussed in detail in Appendix A, but one can apprehend this slow-down as follows: the total Lorentz factor of the flow remains large, in particular the total three-velocity $|\beta| \sim 1$, up to corrections of order $\gamma_{sh}^{-2}$; however, a transverse velocity develops with magnitude $|\beta_y| \simeq \xi_{cr}$; the combination of these two facts implies that $\beta_x$ deviates from unity by quantities of order $\gamma_{sh}^{-2}$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch of the precursor of a relativistic magnetized pair shock, as viewed in the shock front rest frame. The reflected/shock-accelerated suprathermal particles (in red and blue) gyrate in the background magnetic field and accelerate parallel to the convective electric field, thereby generating a net perpendicular current $j_{y,cr}$. The incoming plasma must compensate this current as it penetrates the precursor.</figcaption></figure>

or $\xi_{cr}^2$, whichever is larger. In other words, assuming that $\gamma_{sh}\xi_{cr} \gg 1$, as expected in ultra-relativistic shocks, leads to $|\beta_x| \simeq 1 - \xi_{cr}^2/2$. If $\gamma_{sh}\xi_{cr} \ll 1$, $\beta_x$ remains unchanged compared to the asymptotic value outside the precursor.

This is a quite remarkable feature: the compensation of the current slows down the incoming plasma down to the (longitudinal) velocity $\beta_x$; thus, the $\mathcal{R}$ frame which corresponds to the instantaneous rest frame of the plasma, in which there is no bulk motion along $x$, moves with velocity $\beta_{\mathcal{R}|sh} = \beta_x$ relative to the shock front rest frame. At large values of the current, $\gamma_{sh}\xi_{cr} \gg 1$, the relative Lorentz factor between the $\mathcal{R}$ frame and the shock front rest frame becomes of the order of $1/\xi_{cr}$, independent of the far upstream Lorentz factor. In this sense, the shock precursor plays the role of a buffer, with important consequences for the physics of the shock, discussed in Section 3.

The Lorentz factor that corresponds to the relative velocity between this new rest frame $\mathcal{R}$ and the far upstream rest frame is easily calculated and well approximated by:

$$ \gamma_{\mathcal{R}|u} \approx \max(1, \gamma_{sh}\xi_{cr}/2). \qquad (2) $$

In the following, we analyse the evolution of the instability in the linear regime by adopting a relativistic two-fluid description of the incoming plasma, where two-fluid refers to the electron and positron components of the background plasma. This means, in particular, that we neglect the response of the cosmic rays and we treat as external the current that these suprathermal particles carry. The latter assumption is discussed in Section 3. In this section, we assume that current compensation is achieved to high accuracy in the shock precursor, as motivated by our discussion in Section A1 (see also the discussion in section 5.3.1 of Lemoine & Pelletier 2011). This two-fluid description allows us to probe the physics of the instability up to the inertial scale of the incoming plasma, where the growth rate is found to peak.

We write and solve the system in the instantaneous rest frame $\mathcal{R}$ of the plasma, in which there is no bulk motion along $x$. In such a rest frame, the instability is expected to be absolute (versus convective), provided the growth rate exceeds the inverse crossing time of

the precursor. In the $\mathcal{R}$ frame, $u_{x|\mathcal{R}} = 0$ (henceforth, all quantities concern the incoming plasma), but the (unperturbed) background electric and magnetic fields read

$$ B_{z|\mathcal{R}} = \gamma_{\mathcal{R}|u} B_u, \quad E_{y|\mathcal{R}} = -\gamma_{\mathcal{R}|u} \beta_{\mathcal{R}|u} B_u. \qquad (3) $$

## 2.1 Linear analysis

For simplicity, we assume the plasma and the velocity profile to be uniform throughout the precursor. It is possible to incorporate the terms associated with the variation of the profile by writing the system first in the shock front frame, then boosting it to the instantaneous rest frame of the incoming plasma. The new terms that appear contain spatial derivatives (along $x$) of the various unperturbed quantities. The typical magnitude of these inhomogeneous terms relative to the other terms is of order $\omega_c/\omega$ in Fourier variables; therefore, the above assumption will be justified provided $|k_x| \gg \omega_c/c$. As we show in the following, the growth rate peaks at values close to $\omega_p$ on short wavelengths, i.e. $k \simeq \omega_p/c$; this therefore justifies the above approximation of a uniform precursor.

Our linear analysis is based on a relativistic two-fluid model of the background plasma subject to the external current imposed by the gyrating supra-thermal particles. We thus perturb all variables of the incoming flow and the electromagnetic structure. The unperturbed equations are:

$$ \begin{aligned} \partial_{\mu} (n_{\pm} u_{\pm}^{\mu}) &= 0 \\ \partial_{\mu} T_{\pm}^{\mu\nu} &= \pm en_{\pm} u_{\pm}^{\mu} F_{\nu}^{\mu}. \end{aligned} \qquad (4) $$

The indices $\pm$ refer to the positron/electron species of the background plasma, $u_{\pm}^{\mu}$ to the four-velocity and $T_{\pm}^{\mu\nu}$ to the corresponding energy-momentum tensors. The perturbed system then reads:

$$ \begin{aligned} & u_{\pm}^{\mu} \partial_{\mu} \left( \frac{\delta n_{\pm}}{n} \right) + \partial_{\mu} \delta u_{\pm}^{\mu} = 0 \\ & u_{\pm}^{\mu} \partial_{\mu} \delta u_{\pm}^{\nu} + \beta_{s}^{2} \partial^{\nu} \left( \frac{\delta n_{\pm}}{n} \right) = \pm \frac{e}{m} \delta u_{\pm}^{\sigma} F_{\sigma}^{\nu} \pm \frac{e}{m} u_{\pm}^{\sigma} \delta F_{\sigma}^{\nu}, \end{aligned} \qquad (5) $$

together with the Maxwell equations. We have implicitly assumed a cold background plasma limit, although we incorporate temperature effects through the sound velocity $\beta_s$.

We recombine the two fluid variables $\delta n_{\pm}$ and $\delta u_{\pm}^{\mu}$ into:

$$ \delta n \equiv \frac{\delta n_{+} + \delta n_{-}}{2}, \quad \delta \rho \equiv \frac{\delta n_{+} - \delta n_{-}}{2} \qquad (6) $$

$$ \delta u^{\mu} \equiv \frac{\delta u_{+}^{\mu} + \delta u_{-}^{\mu}}{2}, \quad \Delta u^{\mu} \equiv \frac{\delta u_{+}^{\mu} - \delta u_{-}^{\mu}}{2}. \qquad (7) $$

Of course, to zeroth order, $n_{-} = n_{+} \equiv n$, $u_{-}^{0} = u_{+}^{0} \equiv u^{0}$, $u_{+,-} = -u_{-,+} \equiv u_{y}$. Furthermore, $(u_{\pm}^{\mu} + \delta u_{\pm}^{\mu})(u_{\pm\mu} + \delta u_{\pm\mu}) = -1$ implies

$$ \delta u^{0} = \beta_{y} u^{0} \Delta u_{y}, \quad \Delta u^{0} = \beta_{y} u^{0} \delta u_{y}, \qquad (8) $$

with $\beta_y \equiv u_y/u^0$. In the $\mathcal{R}$ frame, in which we are working here, $u^0 = (1 + u_y^2)^{1/2}$; therefore $u_y \sim \gamma_{sh}\xi_{cr} \gg 1$ at large shock Lorentz factors implies $|\beta_y| \sim 1$. In the limit $\gamma_{sh}\xi_{cr} \gg 1$ (but $\xi_{cr} \ll 1$), the parameters $\gamma_{\mathcal{R}|u}/u^0 \simeq 1/2$ and $\beta_{\mathcal{R}|u} \simeq 1$.

The perturbed current $\delta j^{\mu} = \delta j_{+}^{\mu} + \delta j_{-}^{\mu}$ reads

<footnote n="1">¹ We use a metric with signature $(-, +, +, +)$.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Growth rate $\Im\omega/\omega_p$ versus $k_x c/\omega_p$ for $k_y = k_z = 0$. In solid lines, $\sigma = 10^{-4}$, $\beta_s = 0$ and from top to bottom, $\beta_y = 0.99$ (red, corresponding to $u_y = \gamma_{sh}\xi_{cr} \simeq 7.0$), $\beta_y = 0.1$ (blue, corresponding to $\gamma_{sh}\xi_{cr} = 0.1$), and $\beta_y = 0.03$ (orange, $\gamma_{sh}\xi_{cr} = 0.03$). In dashed lines, same as above for $\sigma = 10^{-2}$ (the growth is strongly suppressed for $\beta_y = 0.03$ in this case).</figcaption></figure>

$$ \delta j^0 = 2\mathrm{sec} (\Delta u^0 + u^0 \delta \rho / n), \qquad (9) $$

$$ \delta j_x = 2\mathrm{sec} \Delta u_x, \qquad (10) $$

$$ \delta j_y = 2\mathrm{sec} (\Delta u_y + \beta_y u^0 \delta n / n), \qquad (11) $$

$$ \delta j_z = 2\mathrm{sec} \Delta u_z. \qquad (12) $$

We define the plasma frequency following: $\omega_p^2 = \omega_{p+}^2 + \omega_{p-}^2 = 8\pi n e^2 / m_e$, and the magnetization parameter:

$$ \sigma = \frac{B_u^2}{8\pi nm_e c^2} = \frac{\omega_c^2}{\omega_p^2}. \qquad (13) $$

The full dispersion relation is calculated from the linear system discussed in Appendix B, by going through Fourier variables, then taking the determinant of the matrix using the Mathematica package. This dispersion relation is too lengthy to be reported here.

However, it can be given in the following form in the 1D approximation $k_y = k_z = 0$, cold plasma limit $\beta_s = 0$:

$$ \omega^5 + i \beta_{\mathcal{R}|u} \beta_y \kappa \omega_c \omega^4 - [k_x^2 c^2 + \kappa^2 \omega_c^2 + (1 - \beta_y^2) \omega_p^2] \omega^3 - i \beta_{\mathcal{R}|u} \beta_y \kappa k_x^2 c^2 \omega_c \omega^2 + (\kappa^2 k_x^2 c^2 \omega_c^2 - \beta_y^2 k_x^2 c^2 \omega_p^2 - i \kappa \beta_y^3 k_x c \omega_c \omega_p^2) \omega - i \beta_{\mathcal{R}|u} \beta_y^3 \kappa k_x^2 c^2 \omega_c \omega_p^2 = 0. $$

We recall here the definition $\kappa \equiv \gamma_{\mathcal{R}|u}/u^0$ (see Appendix B). The growth rate is represented as a function of $k_x$ for various values of the parameters $\beta_y$ and $\sigma$ in Fig. 2. The global trend that emerges is a maximal growth rate

$$ \Im\omega \sim \beta_y \omega_p \quad (k_x c \sim \omega_p, \quad \beta_y \gg \sqrt{\sigma}, \quad \beta_y \gg \beta_s). \qquad (14) $$

The growth rate collapses as soon as one of the conditions indicated in the brackets is no longer satisfied. The last condition $\beta_y \gg \beta_s$ is typical of current-driven instabilities: as the temperature rises and the thermal velocity exceeds the drift velocity, the instability disappears. However, we do not expect this situation in ultra-relativistic

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Contour plot of $\log_{10}(\Im\omega/\omega_p)$ assuming $k_y = 0$, for $\gamma_{sh}u_y = 7$ (i.e. $\beta_y = 0.99$), $\sigma = 10^{-3}$, $\beta_s = 0$.</figcaption></figure>

pair shocks with $\gamma_{sh}\xi_{cr} \gg 1$, since $\beta_y \sim 1$ in that limit, while the heating of the incoming flow inside the precursor remains limited to sub-relativistic velocities [see e.g. Lemoine & Pelletier (2011) for a discussion and Spitkovsky (2008a) for PIC simulations].

In the 2D $k_y = 0$, cold plasma ($\beta_s = 0$), and small current limit ($\xi_{cr}\gamma_{sh} \ll 1$, in which case $\beta_{\mathcal{R}|u} \sim 0$ and $u^0 \sim 1$), the dispersion relation also reduces to the compact form:

$$ \omega^6 - \omega^4 (\omega_p^2 + \omega_c^2 + k_x^2 c^2 + k_z^2 c^2 - \beta_y^2 \omega_p^2) + \omega^2 [(k_x^2 + k_z^2) c^2 \omega_c^2 - \beta_y^2 (k_x^2 + k_z^2) c^2 \omega_p^2 - i \beta_y^3 k_x c \omega_p^2 \omega_c] + \beta_y^2 k_z^2 c^2 \omega_c^2 \omega_p^2 = 0. \qquad (15) $$

In this limit, the instability can be shown to result from a coupling between the high-frequency branch of the extraordinary mode with the acoustic mode, as discussed in the following Section 2.2.

We now present numerical solutions of this dispersion relation in the various 2D planes: ($k_x$, $k_z$) in Fig. 3 assuming $k_y = 0$; ($k_x$, $k_y$) in Fig. 4 assuming $k_z = 0$; and ($k_y$, $k_z$) in Fig. 5 assuming $k_x = 0$.

The global trend that emerges from these numerical simulations is, here as well, a maximum growth rate of order $\beta_y \omega_p$ at wavenumbers $\sim \omega_p$, provided the thermal dispersion velocity $\beta_s$ remains much smaller than the drift velocity $\beta_y$.

## 2.2 Interpretation and analytical approximations

The above instability can be best understood in the limit $k_y = 0$, in the non-relativistic regime $\beta_y \ll 1$, which formally corresponds to $\gamma_{sh}\xi_{cr} \ll 1$. In this limit, one can neglect the acceleration of the plasma relative to the far upstream, $\beta_{\mathcal{R}|u} \sim 0$, so that the convective electric field can be neglected; furthermore, $\kappa = \gamma_{\mathcal{R}|u}/u^0 \sim 1$. Although relativistic shock waves should rather lead to $\gamma_{sh}\xi_{cr} \gtrsim 1$, we find little difference in the growth rate between the above approximation and the numerical calculation, suggesting that it remains a good approximation.
<figure><img src="image_6.png" /><figcaption>Figure 4. Same as Fig. 3, in the plane $(k_x, k_y)$, for $k_z = 0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Same as Fig. 3, in the plane $(k_y, k_z)$, for $k_x = 0$. The growth rate vanishes at large values of $k_y$.</figcaption></figure>

In this $\beta_y \ll 1$ regime, the instability involves only velocity fluctuations $\delta u_x$, $\delta u_z$, a density fluctuation $\delta n$, and electromagnetic perturbations $\delta B_x$, $\delta B_z$ and $\delta E_y$. One then finds that a combination of the acoustic mode along $B_x$ and the high-frequency superluminal branch of the extraordinary mode is destabilized by the drift motion that results from the compensation of the current $j_y$.

To see this, we use the perturbed y-component of the electromagnetic vector potential $\delta A_y$ and the displacement $\xi$ of the plasma. Maxwell equations then imply

$$c^2 \nabla \cdot \nabla \delta A_y - \partial_i^2 \delta A_y + 4\pi c \delta j_y = 0, \quad (16)$$

with $\delta j_y = \delta j_y^{(r)} + \delta j_y^{(c)}$, with the following notations:

$$\delta j_y^{(r)} \equiv 2ecn \Delta u_y \quad (17)$$

and

$$\delta j_y^{(c)} \equiv 2ec \beta_y \delta n. \quad (18)$$

Note the difference between $\Delta u_y$ and $\delta u_y$, defined in equation (7). Note also that $u^0 \sim 1$ because we work here in the rest frame of the ambient plasma under the approximation $\gamma_{sh} \xi_{cr} \ll 1$.

The response current $\delta j_y^{(r)}$ evolves according to the dynamical equation:

$$\partial_t \delta j_y^{(r)} = \frac{\omega_p^2}{4\pi} (\delta E_y - \delta u_x B_u). \quad (19)$$

The perturbed bulk velocity can be written: $\delta u = \partial_t \xi$, and $\delta E_y = -\partial_t \delta A_y / c$. Thus we obtain the simple relation

$$\delta j_y^{(r)} = - \frac{\omega_p^2}{4\pi c} (\delta A_y + B_u \xi_x). \quad (20)$$

The dynamics of the centre of mass is governed by a MHD-type equation (with $\rho = 2 \, \text{nm}$):

$$\rho \partial_t \delta u + \rho c_s^2 \nabla \frac{\delta n}{n} = \frac{1}{c} \mathbf{j} \times \delta B + \frac{1}{c} \delta j_y^{(r)} \times \mathbf{B}_u, \quad (21)$$

with of course, $\mathbf{j} \equiv 2 nec_u \mathbf{y}$. Note that $\delta j^{(c)}$ does not contribute to the Lorentz force because the term in $\delta n$ cancels out with the equilibrium condition. Note also that $\delta B_x = -\partial_z \delta A_y$ and $\delta B_z = \partial_x \delta A_y$. In particular the x-component reads:

$$\rho \partial_t^2 \xi_x + \rho c_s^2 \partial_x \frac{\delta n}{n} = \frac{1}{c} j_y \delta B_z + \frac{1}{c} \delta j_y^{(r)} B_u, \quad (22)$$

which can be rewritten as (introducing $\delta \tilde{A}_y \equiv \delta A_y / B_u$):

$$\partial_t^2 \xi_x + c_s^2 \partial_x \frac{\delta n}{n} = \omega_c \beta_y c \partial_x \delta \tilde{A}_y - \omega_c^2 (\delta \tilde{A}_y + \xi_x). \quad (23)$$

One can use also the z-component; however it turns out that the equation for sound evolution is more convenient; we obtain it by taking the divergence of the dynamical equation:

$$(\partial_t^2 - c_s^2 \Delta) \frac{\delta n}{n} = \omega_c^2 \partial_x (\delta \tilde{A}_y + \xi_x) - \omega_c \beta_y c \Delta \delta \tilde{A}_y. \quad (24)$$

Therefore we have obtained three dynamical equations of second order in time derivative that couple $\delta A_y$, $\delta n$ and $\xi_x$. Equation (16) for $\delta A_y$ can be rewritten as

$$c^2 \Delta \delta \tilde{A}_y - \partial_t^2 \delta \tilde{A}_y - \omega_p^2 \delta \tilde{A}_y - \omega_p^2 \xi_x + \frac{\omega_p^2}{\omega_c} \beta_y \frac{\delta n}{n} = 0. \quad (25)$$

This system leads to the following dispersion relation:

$$[P_X(\omega^2) - \beta_y^2 \omega_p^2 k^2 c^2] \omega^2 + \beta_y^2 \omega_p^2 \omega_c^2 k_z^2 c^2 - k_s^2 c_s^2 [P_X(\omega^2) + \frac{k_x^2}{k^2} \omega_c^2 (\omega^2 - k^2 c^2)] = 0 \quad (26)$$

with

$$P_X(\omega^2) \equiv \omega^4 - (\omega_p^2 + \omega_c^2 + k^2 c^2) \omega^2 + \omega_c^2 k^2 c^2, \quad (27)$$

and $k^2 = k_x^2 + k_z^2$; $P_X(\omega^2) = 0$ gives the dispersion relation of the extraordinary mode in the cold plasma limit.

This dispersion relation matches well equation (15) up to relativistic corrections in $\beta_y$. Let us discuss equation (26) in several limits of interest.
### 2.2.1 Cold and weakly magnetized: $\beta_s^2 \ll \sigma \ll 1$

Let us analyse the instability in the cold plasma limit, and at small values of $\sigma$, however not necessarily smaller than $\beta_y^2$ when this parameter is small. The dispersion relation reduces to:

$$P_X(\omega^2) - \omega_p^2 k^2 \beta_y^2 c^2 = 0. \quad (28)$$

This leads to a negative root in $\omega^2$:

$$\omega_{\perp}^2 = -(\beta_y^2 - \sigma) \omega_p^2 F(k^2 \delta^2), \quad (29)$$

with $\delta \equiv c/\omega_p$, and

$$F(k^2 \delta^2) = \frac{2k^2 \delta^2}{1 + k^2 \delta^2 + [(1 + k^2 \delta^2)^2 + 4(\beta_y^2 - \sigma)k^2 \delta^2]^{1/2}} \quad (30)$$

In the latter expression, the contribution of $\sigma$ must be kept when it is no longer negligible compared to $\beta_y^2$. For $k^2 \delta^2 \ll 1$, $F(k^2 \delta^2) \simeq k^2 \delta^2$ and

$$\omega_{\perp}^2 \approx -\omega_p^2 (\beta_y^2 - \sigma) k^2 \delta^2. \quad (31)$$

For $k^2 \delta^2 \gg 1$, $F(k^2 \delta^2) \simeq 1$ and

$$\omega_{\perp}^2 \approx -\omega_p^2 (\beta_y^2 - \sigma); \quad (32)$$

which gives the maximum growth rate. Clearly the instability occurs at low magnetization, precisely when $\sigma < \beta_y^2$, in very good agreement with the analysis of the previous section.

### 2.2.2 Long wavelength modes, $k^2 \delta^2 \ll 1$ and finite $\sigma$

In this limit $k^2 \delta^2 \ll 1$, we find

$$\omega_{\perp}^2 \approx -\beta_y \omega_p^2 \sqrt{\frac{\sigma}{1+\sigma}} k_z \delta, \quad (33)$$

i.e. a growth rate for small $\sigma$

$$\Im \omega \approx \sqrt{\omega_c k_z \beta_y c}, \quad (34)$$

which extends previous results obtained in Pelletier et al. (2009) and in Casse, Marcowith & Keppens (2013) in the MHD regime for similar configurations (see also Riquelme & Spitkovsky 2010 and Nekrasov 2013 for similar configurations in the non-relativistic limit). It thus indicates that this instability has a kinetic origin and that the MHD solution describes its long wavelength behaviour.

### 2.2.3 Warm plasma with $\sigma \ll \beta_s^2$

From the general dispersion relation we find:

$$\omega_{\perp}^2 = -\omega_p^2 [\beta_y^2 - \beta_s^2(1 + k^2 \delta^2)] F_s(k^2 \delta^2, \beta_s), \quad (35)$$

where

$$F_s(k^2 \delta^2, \beta_s) = 2k^2 \delta^2 \left\{ 1 + (1+\beta_s^2) k^2 \delta^2 + \left[ [1 + (1+\beta_s^2) k^2 \delta^2]^2 + 4 (\beta_y^2 - \beta_s^2(1+k^2\delta^2)k^2\delta^2)^{1/2} \right]^{-1} \right\}, \quad (36)$$

which can be well approximated by

$$F_s(k^2 \delta^2) \approx \frac{k^2 \delta^2}{1 + k^2 \delta^2}. \quad (37)$$

The main conclusion is that temperature effects quench the instability when $\beta_s \gtrsim \beta_y$, as reported in the previous section.

## 2.3 Description and evolution

The instability presents the character of a common current instability that is triggered when the drift velocity is larger than the sound velocity and also the character of a Weibel type electromagnetic instability when the threshold is strongly overstepped. The growth rate can reach values as large as $\omega_p$ and makes the instability faster than all instabilities previously studied, including the filamentation instability triggered by the reflected particles ($\Im \omega \simeq \sqrt{\xi_{cr}}\omega_p$), the oblique two-stream instability ($\Im \omega \simeq \xi_{cr}^{1/3}\omega_p$), etc.

We find that this instability is quenched at high temperatures, when $\beta_y \lesssim \beta_s \sim \sqrt{kT/mc^2}$. However, in the precursor of relativistic shocks, one expects $\beta_y \sim 1$ and for pair shocks, the preheating inside the precursor remains moderate. Therefore, such temperature effects are not expected to contribute strongly.

In the 2D setting $k_y = 0$, this instability leads to filamentation of the plasma in a way similar to the standard Weibel-filamentation instability, with some noticeable differences. In particular, the current perturbation is here produced by a global charge neutral density variation, $\delta j_y^{(c)} = \beta_y \delta nec$, not by a charge perturbation as in the Weibel/filamentation instability. This density variation is itself produced by the compression effect associated with the Lorentz force, derived from the drift $\beta_y$. In contrast, the perturbed current in the Weibel/filamentation instability $\delta j^{(w)} = \beta_w \delta \rho e c$, with $\beta_w$ the drift velocity of two $e^-$ counterstreaming beams (assuming that charge neutralization is ensured, e.g. by ions) and $\delta \rho e$ the charge perturbation (as before). The Lorentz force then couples this charge perturbation to the electromagnetic potential through

$$\partial_t^2 \delta \rho = - \beta_w \frac{\omega_p^2}{4\pi e c} \Delta \delta A_y. \quad (38)$$

In this counterstreaming (symmetric) situation, the Weibel instability gives rise to small-scale magnetic perturbations with a growth rate similar to that of the current-driven filamentation instability. The difference pointed above, namely charge perturbation versus density perturbation, brings in a major difference between these two instabilities, which is related to the polarity of the current filaments. While in the Weibel instability, the counterstreaming beams contain particles of similar charge, which thus deviate in a perturbed magnetic field in different directions to form filaments of opposite current, in the current-driven filamentation instability, the beams contain particles of opposite charge, which thus deviate in the same direction and create filaments with a current oriented in the same direction, i.e. so as to compensate the current of the suprathmal particles. This picture is sketched in Fig. 6. Current-driven filamentation is thus subject to coalescence and reconnection. The non-linear evolution of this instability will be addressed in a forthcoming study (Plotnikov et al., in preparation).

## 3 DISCUSSION

In our treatment of the instability, we have neglected the response of the plasma of suprathmal particles. This choice is dictated by simplicity, as including the response involves doubling the number of fluid variables, which renders the problem intractable. However, one should expect this approximation to be valid at maximal growth rate, since $\Im \omega$ then becomes larger than the plasma frequency of the suprathmal particles, $\omega_{pb} = \xi_{cr}^{1/2}\omega_p$, with $\omega_p$ the plasma frequency of the ambient (upstream) plasma. As the instability develops and turbulence grows, one should of course expect the orbits of these suprathmal particles to deviate from their zeroth order form given
<figure><img src="image_8.png" /><figcaption>Figure 6. Sketch of the development of the current-driven filamentation instability in 1D, with a perturbation along x. The current carried by the suprathermal particles, oriented along −y, is not indicated here. Electrons and positrons of the ambient plasma flow in opposite directions to compensate this current; if present, a magnetic fluctuation along z leads to density enhancements along filaments, thereby creating a current perturbation which feeds back positively on the magnetic fluctuation.</figcaption></figure>

in Appendix A; this influence will be made more precise in the following (Section 3.2).

## 3.1 Relevance to relativistic shocks

Let us now discuss why the current-driven filamentation instability is likely to play a central role in shaping the precursor, the shock and the acceleration process in the relativistic mildly magnetized regime.

Advection through the shock front provides a crucial limitation for the growth of instabilities upstream of a relativistic shock front. In the upstream plasma rest frame, this can be understood as follows: the precursor extends at most to a distance $c/(\gamma_{sh}\omega_c)$ (e.g. Milosavljević & Nakar 2006; Pelletier et al. 2009), because the suprathermal particles only rotate by an amount $1/\gamma_{sh}$ before being caught back by the shock front; this takes a time $t_u \sim \gamma_{sh}\omega_c^{-1}$, but the distance between the shock front and the tip of this precursor does not exceed $t_u(1 - \beta_{sh}^2) \sim t_u/(2\gamma_{sh}^2)$. Therefore, as measured in the upstream plasma rest frame (indicated by $|u$), any instability whose growth rate $\Im\omega_{|u} \lesssim \gamma_{sh}\omega_c$ cannot grow on the crossing time of the precursor. For the filamentation instability, $\Im\omega_{|u} \sim \xi_{cr}^{1/2}\omega_p$, therefore the instability can grow only if $\gamma_{sh}^2\sigma\xi_{cr}^{-1} \lesssim 1$ (Lemoine & Pelletier 2010, 2011). This indicates that mildly magnetized and/or large Lorentz factor shock waves cannot be mediated by the Weibel-filamentation instability, as mentioned in the Introduction.

The present current-driven filamentation instability modifies this picture, because it grows faster than any of the other instabilities discussed in the context of relativistic shocks, and mostly because of the impact of the current on the incoming plasma in the shock front frame: as discussed in Appendix A and Section 2, if $\gamma_{sh}\xi_{cr} \gtrsim 1$, the upstream cannot compensate the current at rest; it is therefore accelerated along x to a Lorentz factor $\gamma_{R|u} \sim \gamma_{sh}\xi_{cr}/2$ in the upstream rest frame, and its apparent density increases by a similar amount. In the shock front frame, the incoming plasma is slowed down to velocities $\beta_{x,in} \sim -(1 - \xi_{cr}^2/2)$, which means that the rest frame of the plasma effectively moves with a Lorentz factor (along $-x$): $\gamma_{R|sh} \sim 1/\xi_{cr}$. This change of rest frame, relative to far infinity, strongly modifies the criterion under which the instability has or does not have time to grow. In the $\mathcal{R}$ frame, which defines the rest frame of the background plasma after its acceleration phase, the shock front moves with a Lorentz factor $\gamma_{R|sh}$, therefore the

precursor size extends to $c/(\gamma_{R|sh}\omega_c)$ and the time-scale for a plasma mode to cross this precursor now reads

$$t_{x|\mathcal{R}} \approx \frac{1}{\gamma_{R|sh}|\beta_{R|sh}|\omega_c}, \qquad (39)$$

so that the instability can grow whenever $\Im\omega t_{x|\mathcal{R}} \gtrsim 1$, or

$$\sigma \lesssim \xi_{cr}^2. \qquad (40)$$

For typical values $\xi_{cr} \sim 0.1$, this implies that growth is possible up to magnetization levels $\sigma \sim 10^{-2}$, irrespective of the Lorentz factor of the shock. The latter point is of importance, because it guarantees the growth of instabilities at large $\gamma_{sh}$, for which the precursor becomes very short in the upstream rest frame. This result appears compatible with recent PIC simulations, as we argue in Section 3.3.

Once micro-turbulence grows upstream of a relativistic collisionless shock, one may expect the Fermi process to develop (e.g. Lemoine et al. 2006; Niemiec et al. 2006) although how well it develops depends on the relative efficiency of scattering in the micro-turbulence relatively to advection in the large-scale field (Pelletier et al. 2009; Lemoine & Pelletier 2010). To discuss this on quantitative grounds, we write the scattering frequency in the downstream rest frame

$$v_s \sim c\lambda_{\delta B}/r_g^2 \sim \epsilon_{B,d} (\lambda_{\delta B}\omega_p/c) \omega_p, \qquad (41)$$

$\epsilon_{B,d}$ denoting an average value of the equipartition fraction of the magnetic field downstream of the shock, $\lambda_{\delta B}$ representing the coherence length of the field; the above equation holds for typical supra-thermal particles of Lorentz factor $\gamma_{sh}$ in the downstream frame. As discussed in Lemoine & Pelletier (2010), scattering beats advection, hence the Fermi process develops, when

$$v_s \gg \omega_c \Leftrightarrow \sigma \ll \epsilon_{B,d}^2 (\lambda_{\delta B} \omega_p / c)^2. \qquad (42)$$

PIC simulations suggest $\epsilon_{B,d} \sim 0.01$ and $\lambda_{\delta B} \sim 1 - 10c/\omega_p$ with some degree of uncertainty. Nevertheless, the above result indicates that the current-driven instability that we are discussing here must also play a key role in the switch-on of the Fermi process, by building up the micro-turbulence for any value of the shock Lorentz factor, up to magnetization levels as high as $\sigma \lesssim 10^{-2}$.

## 3.2 Current-driven instability versus Weibel/filamentation

At very low magnetization levels, one must expect this current-driven filamentation to gradually disappear, once the other more standard (Weibel-filamentation, two stream, etc.) instabilities can grow. To see this, consider the extreme $\sigma \to 0$ limit: the Weibel/filamentation instability then grows, excites turbulence which scatters the suprathermal particles; since this turbulence has no preferred direction in the transverse plane (y, z), no net perpendicular current arises and current-driven filamentation does not take place.

At finite magnetization, the average current does not vanish, but it may be randomized by the micro-turbulence. This effect has not been taken into in the present calculations, which work at linear order and which neglect the response of the cosmic rays. In order to quantify the magnitude of the back-reaction of the turbulence on the particle trajectories, one must compare the upstream residence time derived under the assumption that microturbulence controls the scattering process with that derived assuming a coherent gyration in the background field. Furthermore, this comparison must be made upstream, in the proper frame of the micro-turbulence. In what follows, we assume that this frame corresponds to $\mathcal{R}$. In this $\mathcal{R}$ frame, the turbulent magnetic field strength $\delta B_{\mathcal{R}} \simeq \delta B/\gamma_{R|sh}$ relatively
to that measured in the shock front; similarly, the typical Lorentz factor of a supra-thermal particle can be written $\bar{\gamma}_R \simeq \gamma_{R|sh} \gamma_{sh}$; the background field has a strength $B_{z|R} \simeq \gamma_{R|u} B_u$ (see equation 2). In this $R$ frame, return to the shock takes place once the particle has been scattered by an angle $\delta\theta_R \sim 1/\gamma_{R|sh}$ (see the discussion in Milosavljević & Nakar 2006; Plotnikov et al. 2013). If the supra-thermal particles gyrate coherently in the background electromagnetic field, return occurs on a time-scale $t_{\omega|R} \sim \delta\theta_R/\omega_R$, with $\omega_R \simeq (\gamma_{R|u}/\bar{\gamma}_R)\omega_c$. If micro-turbulence controls the scattering with scattering frequency $v_{s|R} \sim c\lambda_{\delta B|R}/r_g^2$, return takes place on a time-scale $t_{rv|R} \sim \delta\theta_R^2/v_{s|R}$. Comparing the two time-scales leads to a critical magnetization level:

$$ \sigma_c \sim \xi_{cr}^2 \epsilon_{B,u}^2 (\lambda_{\delta B|R} \omega_p / c)^2 . \quad (43) $$

The quantity $\epsilon_{B,u}$ denotes the typical level of micro-turbulence, upstream of the shock, as measured in the shock front frame. The factor $\xi_{cr}$ appears in this formula because the comparison has been made in the $R$ frame. If the upstream magnetization $\sigma \lesssim \sigma_c$, then micro-turbulent scattering efficiently randomizes the trajectories in the shock front plane, hence the perpendicular current as well. Conversely, if $\sigma \gtrsim \sigma_c$, the return trajectories maintain their coherence, hence the current-driven instability develops efficiently.

An interesting question is what happens at large Lorentz factors and low magnetization levels $\sigma \ll \sigma_c$, where feedback from the turbulence should not be neglected, but where the Weibel-filamentation instability does not have time to grow (in the absence of slow-down of the plasma; see below). This area of parameter space corresponds to $\sigma\gamma_{sh}^2\xi_{cr}^{-1} \gtrsim 1$ and $\sigma \lesssim \sigma_c$. Our analysis suggests that the current-driven instability must develop at the tip of the precursor, where the turbulence is sufficiently weak that its back-reaction can be neglected. Furthermore, the deceleration of the plasma, which results from current compensation, now allows the Weibel/filamentation instability to grow: equation (39) indicates that growth becomes possible in the $R$ frame whenever $\sigma \lesssim \xi_{cr}^3$. This instability may then step over closer to the shock front, where the back-reaction of the turbulence strongly randomizes the return trajectories of the supra-thermal particles.

Nevertheless, one expect the precursor to be shaped by the size $c/\omega_c$ if the current-driven instability shapes the precursor, or even the tip of the precursor: beyond that length scale, the turbulence must die away quickly, because the plasma has not yet slowed down and instabilities cannot grow there; inside the precursor, one may expect some form of equilibrium to be reached between the level of the turbulence, the slow-down of the plasma and the growth rate of the instabilities. Its detailed study lies beyond the present work.

This description contrasts with what one expects in the region of parameter space in which the Weibel-filamentation instability can grow without the slow-down imposed by the current, i.e. $\sigma\gamma_{sh}^2\xi_{cr}^{-1} \lesssim 1$ and $\sigma \ll 10^{-5}$. There, as discussed above, the current is mostly randomized by the near isotropicity of the trajectories of suprathermal particles in the shock front plane. In this limit, the precursor extends to a scale $\epsilon_B^{-1}(\lambda_{\delta B}\omega_p/c)^{-1}c/\omega_p$, smaller than $c/\omega_c$, since the return of suprathermal particles is controlled by the scattering in the small-scale turbulence (Milosavljević & Nakar 2006; Pelletier et al. 2009). This situation actually matches the unmagnetized shock limit; hence, one may expect to find a universal precursor profile, independent of the magnetization parameter. The detailed discussion of the profile in this regime is also left open for further study.

### 3.3 Comparison to PIC simulations

PIC simulations offer valuable tools to probe the physics of relativistic collisionless shock waves. So far, most studies have discussed the unmagnetized or strongly magnetized limit and few have addressed the mild magnetization regime, of interest here. We thus confront our findings to the recent simulations of Sironi et al. (2013), which have explored the regime of moderate magnetizations $\sigma = 10^{-4} \rightarrow 10^{-2}$ at various shock Lorentz factors $\gamma_{sh} = 5 \rightarrow 200$. Such simulations are performed in the downstream plasma rest frame, which does not differ much from the shock rest frame. In this rest frame, the slow-down of the plasma along $x$ is difficult to measure, because the relative modification of $u_{x,in}$ is only of an order of $\xi_{cr}$ (see Appendix A).

However, their fig. 7 is particularly interesting, because it reveals a precursor whose profile does not depend on $\sigma$, provided one rescales the distances by $\sigma^{1/2} = \omega_c/\omega_p$, i.e. provided the distances are expressed in units of $c/\omega_c$. It is actually possible to infer directly from their figure the typical scale height of the precursor, $\sim 2c/\omega_c$, with a rough exponential dependence. For the parameters probed in this figure, $\gamma_{sh} = 21 (\simeq \sqrt{2}\gamma_0$ with their $\gamma_0 = 15)$ and $\sigma = 10^{-4} \rightarrow 10^{-3}$, the Weibel-filamentation instability cannot grow without the slow-down of the plasma imparted by the current-driven filamentation. Therefore, these simulations directly probe the region of parameter space discussed above, in which the current-driven filamentation instability plays the central role. The structure of the precursor conforms well to the expectations, with a size $\sim c/\omega_c$.

In their fig. 5, these authors show the magnetic structure of the precursor in 3D simulations for similar parameters; the magnetic field appears to be structured in sheets parallel to the $x-y$ plane rather than filaments oriented along $x$, which would be expected for a standard Weibel/filamentation instability. Finally, they report no dependence on the shock Lorentz factor, whereas a rather strong dependence is expected if the Weibel-filamentation instability alone shapes the precursor: as the line $\sigma\gamma_{sh}^2\xi_{cr}^{-1} = 1$ is crossed, one expects to transit in a region in which the Weibel-filamentation instability can no longer grow. This independence relative to the Lorentz factor directly results from the slow-down imposed by the current compensation in the $\gamma_{sh}\xi_{cr} \gg 1$ limit: inside the precursor, everything happens as if the shock were moving relative to upstream with the Lorentz factor $\gamma_{R|sh} \sim 1/\xi_{cr}$, so that all memory of the initial $\gamma_{sh}$ is lost.

These trends strongly suggest that the present current-driven filamentation instability shapes the precursor and the shock of weakly magnetized ($\sigma \ll 1$) relativistic shock waves.

Finally, the picture that we have elaborated in Section 3.1 also allows us to understand, at least qualitatively, the results of Sironi et al. (2013) concerning the development of Fermi acceleration. Their simulations indicate that Fermi acceleration develops for any value of the shock Lorentz factor, for magnetization levels $\sigma \lesssim 10^{-5}$. This conforms well with equation (42) and the discussion in Lemoine & Pelletier (2010). There, current-driven filamentation can grow, irrespectively of the shock Lorentz factor; it builds up turbulence and, because $\sigma \lesssim 10^{-5}$, scattering in the micro-turbulence downstream of the shock front beats advection, hence the Fermi process develops. At larger values of $\sigma$, the same simulations indicate that Fermi acceleration develops in a restricted dynamic range, with a maximum energy scaling as $\sigma^{-1/4}$. Equation (42), taken at face value, would indicate that Fermi acceleration should not develop in this limit. However, this argument assumes a homogeneous micro-turbulence downstream of the shock, of strength $\epsilon_{B,d}$, whereas the micro-turbulence seen in PIC simulations actually decreases away
from the shock front. If the law of evolution of $\epsilon_B$ were known, one could improve on equation (42) by comparing the scattering time in this evolving micro-turbulence and the gyration time in the background field. In the absence of such a well-defined law, one can nevertheless understand on a qualitative level the scaling of the maximal energy: as the magnetization increases beyond $10^{-5}$, the condition $\sigma \lesssim \epsilon_{B,d}^2$ remains true only in a finite layer close to the shock front; since the scattering length-scale evolves as the square of the particle energy, the restricted size of this layer leads to the existence of a maximal energy. Let us note that if this layer were of infinite extent, there would nevertheless be a maximal energy, scaling as $\sigma^{-1/2}$, as discussed in Pelletier et al. (2009).

### 3.4 Consequences

The above discussion directly impacts our understanding of shock structuration and of particle acceleration. For instance, Sironi et al. (2013) argue that in front of the shock, there exists a layer of size $\sim c/\omega_c$ filled with Weibel turbulence at a level $\epsilon_B \sim 10^{-2}$; this observation is based on the simulations reported above, in the range $\sigma = 10^{-4} \rightarrow 10^{-2}$. According to the above discussion, this layer actually reflects the constrained growth of current-driven filamentation and Weibel-filamentation instabilities in the precursor, whose size is set by the current profile, which extends on $c/\omega_c$, and the turbulence is not of Weibel origin.

These authors then extrapolate their results to the regime of low magnetization $\sigma \ll 10^{-5}$ to discuss the maximal energy of particles accelerated at relativistic shocks. The above arguments indicate that such an extrapolation is not justified, because the physics of the precursor are likely to change as one transits from the region controlled by the current-driven filamentation instability to that controlled solely by the Weibel-filamentation mode. In particular, as $\sigma \rightarrow 0$, the diverging scale $c/\omega_c$ must decouple and one expects the precursor profile to be entirely controlled by the micro-turbulence, as in the unmagnetized limit. The above discussion indicates that this transition takes place close to the line $\sigma\gamma_{sh}^2\xi_{cr}^{-1} \sim 1$ and $\sigma \sim 10^{-5}$ to $10^{-4}$; for $\gamma_{sh} = 21$ as used in these simulations, both limits reduce to the latter $\sigma \sim 10^{-5}$ to $10^{-4}$.

## 4 CONCLUSIONS

This work reports on a new current-driven filamentation instability upstream of a magnetized relativistic collisionless shock front. As viewed in the shock front frame, the suprathermal particles, which are reflected on the shock front, or accelerated at the shock, gyrate around the perpendicular magnetic field in the shock precursor, thereby depositing a strong current $j_{cr} \sim \xi_{cr}\gamma_{sh} n_u e c$, which is perpendicular to both the magnetic field and the shock normal. As the incoming plasma enters the precursor, it seeks to compensate this current within a few skin depth scales. If $\xi_{cr}\gamma_{sh} \gtrsim 1$, which is a likely situation for highly relativistic shocks, the incoming plasma cannot compensate this current in the upstream rest frame; it is thus accelerated to a large Lorentz factor $\sim \xi_{cr}\gamma_{sh}/2$ (relative to far upstream), which increases the apparent density of the plasma by a similar factor; particles then drift at relativistic velocities in the perpendicular direction to achieve current compensation, electrons and positrons drifting in opposite directions. In the shock front rest frame, the incoming plasma is decelerated along the shock normal at the same time as it is accelerated in this perpendicular direction.

As we have argued, this current destabilizes a combination of the high-frequency branch of the extraordinary mode and of the acoustic mode along the magnetic field. In a 2D configuration, in

which one neglects perturbations along the direction of the current, this instability bears some resemblance to the Weibel-filamentation instability. However, in the present case, the electromagnetic perturbation couples to a density fluctuation, not to a charge fluctuation, because the counterstreaming electrons and positrons carry opposite charges. This leads to the formation of current filaments of a same polarity, all currents being oriented so as to compensate the cosmic ray current induced in the precursor. We find that this instability has a very fast growth rate, of order $\gamma\omega \sim \beta_y\omega_p$ on skin depth scales, with $\beta_y \sim 1$ the drift velocity. This instability is likely to play a key role in shaping the precursor of weakly magnetized relativistic collisionless shocks, in which the growth of other instabilities is very often impeded by the fast transit across the precursor.

In particular, we have shown that this instability can grow at any value of the Lorentz factor, provided the magnetization parameter $\sigma \lesssim \xi_{cr}^2 \sim 10^{-2}$. The relative independence to the Lorentz factor of the shock, which controls the size of the precursor $c/(\gamma_{sh}\omega_c)$ (upstream rest frame), stems from the deceleration that the incoming plasma suffers inside the precursor: the relative Lorentz factor between the shock front frame and the rest frame of the plasma now falls to $\gamma_{R|sh} \sim 1/\xi_{cr}$, independent of $\gamma_{sh}$. In this picture, the shock foot plays the role of a buffer that transforms the interaction with the fast incoming flow into a more moderate regime, depending on the parameter $\xi_{cr}$, over a well-defined distance $\xi_{cr}c/\omega_c$ (in the instantaneous rest frame of the incoming plasma).

In previous studies, we have argued that the filamentation, oblique two-stream modes etc., can grow only at small values of $\sigma$ and moderate values of $\gamma_{sh}$, e.g. such that $\sigma\gamma_{sh}^2\xi_{cr}^{-1} \lesssim 1$ for the Weibel-filamentation mode (Lemoine & Pelletier 2010, 2011). Otherwise, the incoming plasma transits faster across the precursor than a growth time of the instability. Therefore, the current-driven filamentation instability emerges as the leading instability outside this region of parameter space. At very low magnetizations, $\sigma \ll 10^{-5}$, and in the region where the standard filamentation mode can grow, the current-driven filamentation instability should gradually disappear, as the turbulent small-scale electromagnetic fields randomize the return trajectories of the suprathermal particles in the shock front plane. In this limit, one transits to the unmagnetized limit, in which the precursor size is no longer controlled by the background magnetic field, but by the profile of the micro-turbulence.

Outside this region, up to $\sigma \sim 10^{-2}$, the current-driven filamentation instability is likely to play a dominant role. The interesting physics of the shock at low magnetizations and at Lorentz factor so large that the standard Weibel-filamentation mode cannot grow deserves close scrutiny. In this region, the current filamentation instability can grow in the absence of strong microturbulence; however the very growth of this instability and of the filamentation mode, thanks to the deceleration of the plasma, builds up the small-scale turbulence, which then back reacts on the current profile. The profile of the precursor in this regime is left open for further study.

Our analysis at linear level indicates that the growth rate of the current-driven filamentation instability is maximal on plasma skin depth scales. This does not affect previous results concerning the maximal energy of accelerated particles, which assume micro-turbulence set on skin depth scales (see e.g. Kirk & Reville 2010; Bykov et al. 2012; Plotnikov et al. 2013).

## ACKNOWLEDGEMENTS

This work has been financially supported by the Programme National Hautes Energies (PNHE).
REFERENCES

Achterberg A., Wiersma J., 2007, A&A, 475, 19

Achterberg A., Wiersma J., Norman C. A., 2007, A&A, 475, 1

Alsop D., Arons J., 1988, Phys. Fluids, 31, 839

Bret A., Gremillet L., Bénisti D., 2010, Phys. Rev. E, 81, 036402

Bykov A., Gehrels N., Krawczynski H., Lemoine M., Pelletier G., Pohl M., 2012, Space Sci. Rev., 173, 309

Casse F., Marcowith A., Keppens R., 2013, MNRAS, 433, 940

Gallant Y. A., Hoshino M., Langdon A. B., Arons J., Max C. E., 1992, ApJ, 391, 73

Haugbølle T., 2011, ApJ, 739, L42

Hoshino M., Arons J., 1991, Phys. Fluids B, 3, 818

Hoshino M., Arons J., Gallant Y. A., Langdon A. B., 1992, ApJ, 390, 454

Keshet U., Katz B., Spitkovsky A., Waxman E., 2009, ApJ, 693, L127

Kirk J., Reville B., 2010, ApJ, 710, 16

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

Lemoine M., Pelletier G., 2011, MNRAS, 417, 1148

Lemoine M., Pelletier G., Revenu B., 2006, ApJ, 645, L129

Lemoine M., Pelletier G., Gremillet L., Plotnikov I., 2014, Europhys. Lett., submitted

Lyubarsky Y., Eichler D., 2006, ApJ, 647, L1250

Martins S. F., Fonseca R. A., Silva L. O., Mori W. B., 2009, ApJ, 695, L189

Medvedev M. V., Loeb A., 1999, ApJ, 526, 697

Milosavljević M., Nakar E., 2006, ApJ, 651, 979

Nekrasov A. K., 2013, Plasma Phys. Control. Fusion, 55, 085007

Niemiec J., Ostrowski M., Pohl M., 2006, ApJ, 650, 1020

Nishikawa K.-I. et al., 2009, ApJ, 698, L10

Pelletier G., Lemoine M., Marcowith A., 2009, MNRAS, 393, 587

Plotnikov I., Pelletier G., Lemoine M., 2013, MNRAS, 430, 1208

Rabinak I., Katz B., Waxman E., 2011, ApJ, 736, 157

Riquelme M., Spitkovsky A., 2010, ApJ, 717, 1054

Shaisultanov R., Lyubarsky Y., Eichler D., 2012, ApJ, 744, 182

Sironi L., Spitkovsky A., 2009, ApJ, 698, 1523

Sironi L., Spitkovsky A., 2011, ApJ, 726, 75

Sironi L., Spitkovsky A., Arons J., 2013, ApJ, 771, 54

Spitkovsky A., 2008a, ApJ, 673, L39

Spitkovsky A., 2008b, ApJ, 682, L5

Wiersma J., Achterberg A., 2004, A&A, 428, 365

APPENDIX A: PROFILE OF THE PRECURSOR

We construct here the profile of the precursor in the cold plasma limit, in the shock rest frame. We seek here a 1D zeroth-order stationary solution of the shock precursor, with $\partial_t = \partial_y = \partial_z = 0$, as dictated by the geometry of the problem. The perturbation of this solution leads to the linear system discussed in Section 2.1 and Appendix B. As discussed in Section 2, the zeroth order solution is characterized by the profile of the magnetic field $B = B_z z$, the convective electric field $E = E_y y$, and the fluid four-velocities of the various species.

Note that the x-component of the current density vanishes for both incoming particles and suprathelial particles, as a consequence of the stationary state: current conservation $\partial_\mu j^\mu_\alpha = 0$ for any species $\alpha$ implies the conservation law $\partial_x j^x_\alpha = 0$; since the x-component of the current density of incoming particles, summed over electrons and positrons, vanishes as $x \to +\infty$, it also vanishes in the precursor, and similarly for the suprathelial particles. As particles gyrate in the $(x, y)$ plane, we set $u^z_\alpha = 0$, hence $j_{z,\alpha} = 0$ for all species.

Furthermore, we do not expect any non-zero $E_x$ component to emerge inside the precursor because of the charge symmetry of the

pair plasma. One can check that the above solution is self-consistent. In particular, the magnetic field does not possess other components, as a result of $\nabla \cdot B = 0$, $\partial_y = \partial_z = 0$ and $j_z = 0$.

A1 Simplified MHD model

In Section 2, we provide a relativistic two-fluid description of the instability, the term two-fluid referring to the electrons and positrons of the incoming background plasma. This description thus extends beyond any MHD picture of the instability, up to the inertial scale of the pair plasma. Nevertheless, it is instructive to describe briefly the structure of the precursor in an ideal MHD picture, in which one assumes that the magnetic field remains frozen in the plasma all throughout the precursor.

Treating the suprathelial particle component as a tenuous fluid carrying a current density $j_{cr} = j_{y,cr} y$, with $j_{y,cr} = -\gamma_{sh}\xi_{cr}n_u e c$, the electric field is fixed through the frozen-in condition:

$$E_y = \bar{\beta}_{x,in} B_z, \qquad (A1)$$

with $\bar{\beta}_{x,in}$ denoting the centre-of-mass three-velocity x-component of the incoming plasma. Then $\nabla \times E = 0$ imposes $\partial_x(\bar{\beta}_{x,in} B_z) = 0$, or

$$\bar{\beta}_{x,in} B_z = \beta_{sh} B_\infty, \qquad (A2)$$

with $B_\infty = \gamma_{sh} B_u$. To keep the analysis brief, here, we assume $\gamma_{sh}\xi_{cr} \ll 1$, meaning that the velocity of the electrons/positrons of the background plasma along the $y$ direction is much smaller than $c$. This allows us to set $\bar{\beta}_{x,in} \simeq \bar{u}_{x,in}/(1 + \bar{u}^2_{x,in})^{1/2}$ in the above equations, with $\bar{u}_{x,in}$ the x-component of the centre-of-mass four-velocity.

The current density flowing in the incoming background plasma is itself fixed through

$$j_{y,in} = -j_{y,cr} - \frac{c}{4\pi} \partial_x B_z. \qquad (A3)$$

Particle number conservation $\partial_x(n_u u_x) = 0$ and energy-momentum conservation in the cold plasma limit then lead to the equation:

$$n_u m_e c^2 \bar{u}_{x,in} \partial_x \bar{u}_{x,in} = \frac{1}{c} j_{y,in} B_z. \qquad (A4)$$

This equation of motion becomes an equation for $\bar{u}_{x,in}$, once equations (A2) and (A3) have been taken into account. This equation can be rewritten in the following compact form:

$$\left[ 1 - \frac{\beta_{sh}^2 \gamma_{sh}^2 \bar{\beta}_{x,in}^2}{\bar{u}_{x,in}^4} \sigma \right] \bar{\beta}_{x,in} \bar{u}_{x,in} \partial_x \bar{u}_{x,in} = \beta_{sh} \gamma_{sh} \xi_{cr} \frac{\omega_c}{c}. \qquad (A5)$$

In this equation, we have used the definition of the magnetization parameter, equation (13) and $\omega_c = eB_u/(m_e c) = eB_\infty/(\gamma_{sh}m_e c)$. Equation (A5) is particularly useful, because it allows us to obtain a quick estimate of the slow-down of the plasma due to the Lorentz force: one first notes that the second term in the brackets, which originates from the uncompensated part of the current in the precursor, is much smaller than unity, and can be safely neglected; then, one finds that between entry into the precursor and shock crossing, the variation of $\bar{u}_{x,in}$ reads

$$\Delta \bar{u}_{x,in} \approx -u_{\infty} \xi_{cr}, \qquad (A6)$$

with $u_\infty = \gamma_{sh} \beta_{sh}$; note that the scale of variation is set by the precursor size $c/\omega_c$. Assuming now that the transverse three-velocity
of electrons and positrons, along y, is of order $\pm\gamma_{sh}\xi_{cr}$, one can check that the three-velocity $\bar{\beta}_{x,in,0}$ close to the shock front is of the order of

$$ \bar{\beta}_{x,in,0} \approx \beta_{sh} \left( 1 - \frac{\xi_{cr}^2}{2} \right). \qquad (A7) $$

These results will remain true in the following multi-fluid description, even at large values of the quantity $\gamma_{sh}\xi_{cr}$. In the above MHD model, $\partial_x B_z/B_z = -\partial_x \bar{\beta}_{x,in}/\bar{\beta}_{x,in}$, therefore the above scalings allow us to derive an estimate of $(c/4\pi)\partial_x B_z$, which characterizes the departure from current compensation in the precursor:

$$ \left| \frac{4\pi j_{y,cr}}{c \partial_x B_z} \right| \approx \frac{1}{\xi_{cr} \sigma} \gg 1, \qquad (A8) $$

indicating that the current is indeed compensated to very high accuracy in the precursor.

## A2 Multi-fluid model

We now turn to a more exhaustive multi-fluid model of the precursor, which is necessary to construct the steady state on which the linear analysis of Section 2 relies. In particular, we relax the frozen-in condition of the magnetic field inside the precursor and we follow the kinematics of the various particle populations along the y direction. Of course, well outside the precursor, one still assumes $E_y = \beta_{sh}B_z$, corresponding to the assumption of zero electric field in the rest frame of the background plasma as $x \to +\infty$.

We consider the following populations of particles: the incoming particles, denoted by the subscript in, and the suprathermal particle population, which we divide into two sub-populations, those moving towards $+x$ from the shock front up to the tip of the precursor (subscript r+) and those moving towards $-x$ from the tip of the precursor towards the shock front (subscript r−). We set the shock front at $x=0$ and the tip of the precursor at $x_1$. All throughout this section, we denote by $u_\alpha^\mu$ the four-velocity of the positron component of species $\alpha$, with $\alpha \in (\text{in}, r+, r−)$. As discussed in Section 2, the x-components of the four-velocities of the electrons match those of the positrons, while the y-components are opposite.

Alsop & Arons (1988) have described the structure of the precursor of a strongly magnetized relativistic shock; they do so by solving the fluid and Maxwell equations with one population of incoming particles, which gyrate in the compressed magnetic field. The present description is slightly different: we set a boundary at $x=0$, corresponding to the shock transition, into which the incoming population flows and out of which the suprathermal particle population emerges, with no specific relation between these two populations.

In the cold plasma limit, the coherent rotation of the suprathermal particles at the tip of the precursor implies $\beta_{x,r\pm}(x_1) = 0$, therefore $n_{r\pm}(x_1) \to +\infty$, and consequently $|j_y(x_1)| \to \infty$. This singular behaviour disappears of course when warm plasma effects are introduced. Indeed, the suprathermal particle population should be described in the present shock rest frame as a relativistically hot plasma with mean Lorentz factor $\sim\gamma_{sh}$ and roughly isotropic distribution function. Such effects are discussed in the next Appendix A3. The cold plasma approximation, which we use here, has the advantage of providing quantitative estimates for the various quantities used in the manuscript.

The electromagnetic profile is thus determined by $B_z \equiv \gamma_{sh}B_u(1+b)$, by the current $j_y$ and the four-velocities of the respective fluids.

This profile of the precursor can be solved as a shooting problem, with three parameters to be determined by the boundary conditions: $b_1, \gamma_{r1}$, corresponding respectively to the deviation from $\gamma_{sh}B_u$, the Lorentz factor of suprathermal particles at the tip of the precursor, and $x_1$. The boundary conditions are:

$$ \begin{aligned} u_{x,r+,0} &= u_{sh}, & u_{y,r+,0} &= 0, \\ n_{r-,0} u_{y,r-,0} &= -n_{in,0} u_{y,in,0}. \end{aligned} \qquad (A9) $$

The first two conditions specify the initial data for the suprathermal particle population: we have chosen here a normal incidence to the shock front and a Lorentz factor $\gamma_{sh}$, as expected at relativistic shocks. The third condition imposes a vanishing net flux of particles along the shock front in the y direction.

In the cold plasma limit, and under the stationary state approximation $\partial_t = 0$, the fluid equations $\partial_\mu(n_\alpha u_\alpha^\mu) = 0$ and $\partial_\mu T_\alpha^{\mu\nu} = +e n_\alpha u_\alpha^\mu F_\mu^\nu$ (for the positron components) read:

$$ \begin{aligned} \partial_x (n_\alpha u_\alpha^x) &= 0, \\ \beta_{x,\alpha} \partial_x u_\alpha^x &= \frac{e}{m_e} \beta_{\alpha,y} B_z, \\ \beta_{x,\alpha} \partial_x u_\alpha^y &= \frac{e}{m_e} (E_y - \beta_{\alpha,x} B_z), \\ \beta_{x,\alpha} \partial_x \gamma_\alpha &= \frac{e}{m_e} \beta_{y,\alpha} E_y. \end{aligned} \qquad (A10) $$

Here, $\gamma_\alpha \equiv u_\alpha^0$. For the various species, the continuity equations imply that at each point: $n_{in}u_x = n_\infty u_\infty$ with $u_\infty = \gamma_{sh}\beta_{sh} < 0$, $n_{r+}u_{x,r+} = \xi_{cr}n_\infty u_{x,r+,0}$, $n_{r-}u_{x,r-} = -\xi_{cr}n_\infty u_{x,r+,0}$. The quantity $n_\infty$ represents the proper particle density as $x \to +\infty$, while $u_{x,r \pm 0}$ represents the x-component of the four-velocity of species $r \pm$ at the shock front.

Complemented with Ampère's law $\partial_x B_z = -4\pi j_y/c$, the system equation (A10) may then be rewritten:

$$ \begin{aligned} \beta_{x,\alpha} \partial_x \beta_{x,\alpha} &= \frac{\omega_{L,\alpha}}{c} [(1+b) - \beta_{x,\alpha} \beta_{sh}] \beta_{y,\alpha}, \\ \beta_{x,\alpha} \partial_x \beta_{y,\alpha} &= \frac{\omega_{L,\alpha}}{c} [\beta_{sh}(1 - \beta_{y,\alpha}^2) - \beta_{x,\alpha}(1+b)], \\ \beta_{x,\alpha} \partial_x \gamma &= \frac{\omega_{L,\alpha}}{c} \gamma \beta_{y,\alpha} \beta_{sh}, \\ \partial_x b &= -\frac{\omega_c}{\sigma \gamma_{sh} n_\infty c} (n_{in} u_{y,in} + n_{r+} u_{y,r+} + n_{r-} u_{y,r-}), \end{aligned} \qquad (A11) $$

with $\omega_{L,\alpha} \equiv e\gamma_{sh}B_u/(\gamma_\alpha mc) = \omega_c\gamma_{sh}/\gamma_\alpha$, in terms of $\omega_c \equiv eB_u/(mc)$ the upstream cyclotron frequency, which sets the spatial scale $c/\omega_c$ of the precursor. As discussed above, $u_{y,in}, u_{y,r+}$ and $u_{y,r-}$ represent the y-components of the four-velocities of the incoming, suprathermal r+ and r− positron components, respectively. The last equation for b implicitly uses the fact that the y-velocities of electrons are opposite to those of the positrons, for both incoming and suprathermal particles, hence their y-current densities add up; the magnetization $\sigma$ is defined in equation (13). This last equation holds in the shock precursor where the various populations mix.

Given the above three parameters $b_1, \gamma_{r1}$ and $x_1$, these fluid equations must then be matched to the boundary conditions; this determines the profile of the precursor.
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Structure of the precursor for $\sigma = 0.01$, $\xi_{cr} = 0.1$ and $\gamma_{sh} = 1000$ ($r_L \equiv c/\omega_c$). Top panel: spatial profiles of the y-current carried by the suprathelial particle population in units of $en_u c$ (in blue), and of the compensating current carried by the inflowing ambient plasma (in red); in green, the spatial profile of the perturbed magnetic field $b = B_z/B_{z,0} - 1$. Bottom panel: spatial profiles of the x-velocities of the suprathelial particle population (in blue), of the inflowing background plasma (in red), and of the y-velocity of the background plasma positrons (in green).</figcaption></figure>

Numerical examples of the profile are represented in Fig. A1. We have set $\sigma = 0.01$, $\xi_{cr} = 0.1$ and $\gamma_{sh} = 10^3$, but the profile does not depend on $\gamma_{sh}$ in the ultra-relativistic limit; it is entirely controlled by $\sigma$ and $\xi_{cr}$.

One can obtain an approximation to the above profile as follows. In the vicinity of $x_1, n_{r\pm}(x) \gg n_{in}(x)$, therefore the incoming particle contribution to Ampère's law can be neglected. Furthermore, one can approximate the motion of r+ particles close to $x_1$ as uniform deceleration, implying $\beta_x \simeq [2|\dot{\beta}_x(x_1)|(x-x_1)]^{1/2}$, with $|\dot{\beta}_x(x_1)| = \omega_{L1}(1+b_1)\beta_{y,r+}(x_1)$ given that $\beta_x(x_1) = 0$, and $\omega_{L1} = \omega_{L,r+}(x_1)$. This allows us to determine the singular profile of the density close to $x_1$, using the continuity equation. Plugging this result and the similar estimate for r- particles into Ampère's law, one derives

$$b \approx b_1 \left[ 1 + \frac{\xi_{cr}}{\sigma} \frac{|\beta_{x,r+}(0)| \sqrt{\dot{\beta}_{y,r+}(x_1)}}{\sqrt{2}(1+b_1)^{1/2}b_1} \omega_{L1}^{1/2} (x-x_1)^{1/2} \right] . \quad (A12)$$

The term in the brackets determines the scale over which $b$ varies close to $x_1$, $\Delta x \sim \sigma^2 \xi_{cr}^{-2} b_1^3 c / (\sqrt{2}\omega_{L1})$. Using Ampère's law with $\partial_x b \sim b_1 / \Delta x$, $\omega_{L1} \sim \omega_c$ and assuming $b \gg 1$ leads to

$$b_1 \sim \left( \frac{\xi_{cr}}{\sigma} \right)^{1/2} . \quad (A13)$$

The above turns out to provide the correct scaling seen in the numerical calculations. In turn, this leads to $\Delta x \sim \sigma^{1/2} \xi_{cr}^{-1/2} c / \omega_c \sim \xi_{cr}^{-1/2} c / \omega_p$: current compensation takes place on skin depth scales, as anticipated in Lemoine & Pelletier (2011).

Outside the precursor, the field goes down to its asymptotic far upstream value on skin depth scales as well. Equations (A11) can be used in this region, with $n_{r\pm} \to 0$ in Ampère's law. As discussed in Alsop & Arons (1988), the system then admits the two integrals of motion

$$\gamma_{in} = \gamma_{sh} (1 - \sigma b),$$

$$u_{x,in} = u_{\infty} \left[ 1 - \frac{\sigma}{2\beta_{sh}^2} b(b+2) \right] . \quad (A14)$$

These two integrals, combined with equations (A11), allow us to derive the following equation for the profile of $b$:

$$\partial_x b = - \frac{\omega_c}{\sqrt{\sigma c}} \frac{[b^2 - \sigma b^2 / (\gamma_{sh}^2 \beta_{sh}^2) - \sigma b^3 / \beta_{sh}^2 - \sigma b^4 / (4\beta_{sh}^2)]^{1/2}}{1 - \sigma b(b+2) / (2\beta_{sh}^2)} . \quad (A15)$$

This equation reveals the length scale of the profile: $c/\omega_p$, and allows us to solve for $b$, by integrating from $b_1$ up to $+\infty$, then for $u_{in}$.

Using the integrals of motion, one computes the typical change in Lorentz factor at the entrance into the precursor,

$$\gamma(x_1) = \gamma_{sh} \left(1 - \sqrt{\sigma \xi_{cr}}\right), \\ u_{x,in}(x_1) \approx u_{\infty} \left(1 - \frac{\xi_{cr}}{2}\right), \quad (A16)$$

$$|u_{y,in}(x_1)| \approx -u_{\infty} \sqrt{\xi_{cr}}. \quad (A17)$$

The variation in Lorentz factor is small compared to that of $u_x$ and $u_y$, but the slow-down along $x$ is substantial: at $x_1$, the particles move at velocity $\beta_{x,in}(x_1) \simeq 1 - \xi_{cr}/2$ in the shock front frame. This slow-down is obvious in Fig. A1.

Well inside the precursor, current compensation implies

$$|u_{y,in}| \approx \xi_{cr} \gamma_{sh}. \quad (A18)$$

In order to derive the slow-down imparted to incoming particles, one first notes that $b \ll 1$ outside the peak at the tip of the precursor, as indicated by equations (A12) and (A13). The dynamics of incoming particles is then given by equation (A11) with $b \ll 1$, which implies that the flow is slowed by an amount

$$\Delta u_{x,in} \approx \gamma_{sh} \xi_{cr}, \quad (A19)$$

between the far upstream value and the value of $u_{x,in}$ well inside the precursor. This value matches that at entry into the precursor, equation (A16), and it also matches the value obtained in the simplified MHD model, equation (A6). This slow-down appears as a direct consequence of current compensation, which imposes a Lorentz force directed in the $+x$ direction. In a similar way, one derives $\Delta\gamma_{in} \sim -\gamma_{sh}\xi_{cr}$. Thus the Lorentz factor of both flows remains large after its modification by the Lorentz force. In terms of three-velocity, this implies that $\beta^2$ remains close to unity, up to $1/(2\gamma_{sh}^2)$.
Using equations (A18) and (A19), one derives the three-velocities well inside the precursor:

$$|\beta_y| \approx \xi_{cr}, \qquad (A20)$$

and, at large values of $\gamma_{sh}$,

$$\beta_{x,in} \approx \beta_{sh} (1 - \xi_{cr}^2/2). \qquad (A21)$$

Therefore, if $\gamma_{sh}\xi_{cr} \gg 1$, then $\beta_{x,in} \simeq -(1 - \xi_{cr}^2/2)$, while $\beta_{x,in} \simeq \beta_{sh}$ in the opposite limit, which corresponds to negligible, sub-relativistic deceleration.

Assuming that $\gamma_{sh}\xi_{cr} \gg 1$, the relative Lorentz factor between the shock front frame and the frame $\mathcal{R}$ in which the incoming is at rest along $\mathbf{x}$, i.e. $u_{x,in}|_{\mathcal{R}} \equiv 0$, has fallen from $\gamma_{sh}$ outside the precursor down to

$$\gamma_{\mathcal{R}|sh} \approx \frac{1}{\xi_{cr}}. \qquad (A22)$$

As viewed in the upstream rest frame, the ambient plasma has been picked up by the current layer and accelerated towards $+\mathbf{x}$ to a Lorentz factor

$$\gamma_{\mathcal{R}|u} \approx \gamma_{sh} \frac{\xi_{cr}}{2} \quad (\gamma_{sh} \xi_{cr} \gg 1). \qquad (A23)$$

Finally, in the $\mathcal{R}$ frame in which the ambient plasma is at rest, the particles move with velocity $|\beta_{y,in}| \sim 1$ with bulk Lorentz factor $\sim \gamma_{sh}\xi_{cr}/2$, provided of course that $\gamma_{sh}\xi_{cr} \gg 1$. In the opposite (weak current) limit, $\gamma_{sh}\xi_{cr} \ll 1$, one finds $|\beta_{y,in}| \sim \gamma_{sh}\xi_{cr}$, $\beta_{\mathcal{R}|u} \sim 0$ and $\gamma_{\mathcal{R}|u} \sim 1$; similarly, $\gamma_{\mathcal{R}|sh} \sim \gamma_{sh}$.

### A3 Warm plasma limit

The above discussion assumed a cold plasma of returning particles, with initial momentum (on the shock surface) directed along the shock normal. Here, we introduce the effects of angular dispersion of the beam of returning particles. For simplicity, we neglect the dispersion in Lorentz factor of the returning particles; this dispersion can be taken into account but it should not modify strongly the overall shape of the current profile.

The number density of returning particles at the shock front (considering $e^+/e^-$ species altogether), with momentum oriented within a solid angle element $d\Omega_i$, is written as $dn_{r+,i}(\Omega_i)$. The magnitude of the current deposited by those particles in the precursor can be written as:

$$d j_{y,r+}(x) = |\beta_y(x)| dn_{r+}(x, \Omega_i)ec. \qquad (A24)$$

Assuming that the particle population r-deposits the same amount of current as r+, and using the equation of conservation for the number density of r+ particles, the total current element deposited by supra-thermal particles emitted in the $\Omega_i$ direction reads:

$$d j_y(x) \approx 2 \frac{|\beta_{y,r+}(x)|}{\beta_{x,r+}(x)} \beta_{x,r+,i} dn_{r+,i}(\Omega_i)ec, \qquad (A25)$$

$\beta_{x,r+,i}$ denoting the initial $x$-component of the three-velocity of r+ particles. This equation can be simplified using the result of the previous section, which indicate that $\beta_{x,r+}(x) \simeq |2\dot{\beta}_{x,r+}(x_1)(x - x_1)|^{1/2}$ in the vicinity of the turning point $x_1$, so that most of the current $d j_y(x)$ is deposited at $x_1$. Note that $x_1$ depends on the initial direction $\Omega_i$. We then approximate the spatial profile of the current element equation (A25) with a delta function in $x$:

$$d j_y(x) \approx \mathcal{A} \delta(x - x_1) d F_{r+,i}(\Omega_i) e \qquad (A26)$$

<figure><img src="image_10.png" /><figcaption>Figure A2. Upper panel: profile of $j_y$ (in units of $\gamma_{sh}\xi_{cr}n_{\infty}ec$) carried by returning particles as a function of $x/r_L$ ($r_L = c/\omega_c$) in the limit $\gamma_{sh} \gg 1$, including the effects of angular dispersion at the shock. Lower panel: modulus of the four-velocity components $|u_{x,in}|$ (solid red), $|u_{y,in}|$ (dashed blue) and $\gamma_{\mathcal{R}|sh}$ (dotted orange); the lower panel assumes $\gamma_{sh} = 100$ and $\xi_{cr} = 0.1$.</figcaption></figure>

with $d F_{r+,i}(\Omega_i) = \beta_{x,r+,i} c dn_{r+,i}(\Omega_i)$ the initial flux element. The prefactor is calculated by normalizing the integrated current element along $x$ in equation (A26) to that obtained in equation (A25).

In order to express $\mathcal{A}$ as a function of the initial velocities $\beta_{x,r+,i}$ and $\beta_{y,r+,i}$, one needs to express the quantity $|\beta_{y,r+}(x)/\beta_{x,r+}(x)|$ in the vicinity of $x_1$ using the equations of motion. These equations of motion must be written in the upstream rest frame then Lorentz transformed to the shock frame. We compute the trajectories of the returning particles in the background electromagnetic field, neglecting in particular the perturbed component of the magnetic field; this should remain a good approximation, given that the overall effect of the angular dispersion of the beam is to spread out over the precursor length scale the current profile. One then obtains first the turning point:

$$x_1(\Omega_i) = \frac{c}{\omega_{c,0}} \gamma_{sh}^3 (1 - \beta_{sh} \beta_{x,r+,i}) \times [\beta_{x,i|u} \sin \omega_1 + \beta_{y,i|u} (1 - \cos \omega_1) + \beta_{sh} \omega_1], \qquad (A27)$$

as a function of the upstream-frame initial velocities

$$\beta_{x,i|u} = \frac{\beta_{x,r+,i} - \beta_{sh}}{1 - \beta_{x,r+,i}\beta_{sh}}, \quad \beta_{y,i|u} = \frac{\beta_{y,r+,i}}{\gamma_{sh}(1 - \beta_{x,r+,i}\beta_{sh})}, \qquad (A28)$$

and the quantity $\omega_1$, which is defined implicitly by:

$$\beta_{x,i|u} \cos \omega_1 + \beta_{y,i|u} \sin \omega_1 = -\beta_{sh}. \qquad (A29)$$

Recall that $\beta_{sh} < 0$ in our present notations. The initial cyclotron frequency of the returning particles reads $\omega_{c,0} = e\gamma_{sh}B_u/(\gamma_{0,r+}mc)$, with $\gamma_{0,r+} \simeq \gamma_{sh}$ their initial Lorentz factor. One derives eventually:

$$\mathcal{A} = 2\sqrt{2} \frac{c}{\omega_{c,0}} \left[ \gamma_{sh} (1 - \beta_{sh} \beta_{x,r+,i}) \right]^{1/2} \times \left| \left( \beta_{x,r+,i} \sin \omega_1 - \beta_{y,r+,i} \cos \omega_1 \right) \right|^{1/2}. \qquad (A30)$$
Finally, the flux is normalized through $\int dF_{r+,i}(a\Omega_i) = \gamma_{sh}\xi_{cr}n_{\infty}c$.

In the limit $\gamma_{sh} \to +\infty$, all above quantities reach finite asymptotes, as it should; we use these asymptotic values in the numerical calculation of the integral over the angular variables. One finally obtains the current profile depicted in Fig. A2.

This profile allows us to estimate the velocity profile of the incoming plasma inside the foot. As in the cold plasma limit, current compensation imposes the following scalings inside the precursor

$$|u_{y,in}| \sim \xi_{cr}\gamma_{sh}, \quad u_{x,in} \sim -(1 - \xi_{cr})\gamma_{sh}, \qquad (A31)$$

so that the relative Lorentz factor between the shock front frame and the frame $\mathcal{R}$ in which the incoming is at rest along $+x$ is, as before, $\gamma_{\mathcal{R}|sh} \simeq 1/\xi_{cr}$ if $\gamma_{sh}\xi_{cr} \gg 1$. Fig. A2 shows a numerical calculation of the evolution of $u_{x,in}$, $|u_{y,in}|$ and $\gamma_{\mathcal{R}|sh}$ inside the precursor (assuming $|b| \ll 1$) for $\gamma_{sh} = 100$ and $\xi_{cr} = 0.1$, which confirms the above scalings.

## APPENDIX B: LINEAR SYSTEM

We explicit here the linear system used to compute the dispersion relation, for reference. We rescale the time and space derivatives by $\omega_c$ (cyclotron frequency in the upstream rest frame): $\partial_{\tilde{t}} \equiv \omega_c^{-1}\partial_t$, $\partial_{\tilde{x}} \equiv c\omega_c^{-1}\partial_x$ etc. We rescale all electromagnetic fields by the background value $B_{z\mathcal{R}}$ (e.g. $\delta\tilde{B}_x \equiv \delta B_x/B_{z\mathcal{R}}$) and we introduce the notations: $\kappa \equiv \gamma_{\mathcal{R}|sh}/u^0$, $\tilde{\beta}_s^2 \equiv \beta_s^2/u^0$, $\delta_n \equiv \delta n/n$, $\delta_\rho \equiv \delta\rho/n$, and we rescale $\delta u^\mu$ and $\Delta u^\mu$ by $u^0$, e.g. $\delta\tilde{u}^\mu \equiv \delta u^\mu/u^0$. This leads

to the following adimensional system

$$\partial_{\tilde{t}}\delta_n + \beta_y\partial_{\tilde{y}}\delta_\rho + \beta_y\partial_{\tilde{t}}\Delta\tilde{u}_y + \partial_{\tilde{\bar{s}}}\delta\tilde{u}_x + \partial_{\tilde{y}}\delta\tilde{u}_y + \partial_{\tilde{z}}\delta\tilde{u}_z = 0$$

$$\partial_{\tilde{t}}\delta_\rho + \beta_y\partial_{\tilde{y}}\delta_n + \beta_y\partial_{\tilde{t}}\delta\tilde{u}_y + \partial_{\tilde{x}}\Delta\tilde{u}_x + \partial_{\tilde{y}}\Delta\tilde{u}_y + \partial_{\tilde{z}}\Delta\tilde{u}_z = 0$$

$$\partial_{\tilde{t}}\delta\tilde{u}_x + \beta_y\partial_{\tilde{y}}\Delta\tilde{u}_x + \tilde{\beta}_s^2\partial_{\tilde{x}}\delta_n - \kappa\Delta\tilde{u}_y - \beta_y\kappa\delta\tilde{B}_z = 0$$

$$\partial_{\tilde{t}}\Delta\tilde{u}_x + \beta_y\partial_{\tilde{y}}\delta\tilde{u}_x + \tilde{\beta}_s^2\partial_{\tilde{x}}\delta_\rho - \kappa\delta\tilde{u}_y - \kappa\delta\tilde{E}_x = 0$$

$$\partial_{\tilde{t}}\delta\tilde{u}_y + \beta_y\partial_{\tilde{y}}\Delta\tilde{u}_y + \tilde{\beta}_s^2\partial_{\tilde{y}}\delta_n + \beta_{\mathcal{R}|u}\kappa\beta_y\delta\tilde{u}_y + \kappa\Delta\tilde{u}_x = 0$$

$$\partial_{\tilde{t}}\Delta\tilde{u}_y + \beta_y\partial_{\tilde{y}}\delta\tilde{u}_y + \tilde{\beta}_s^2\partial_{\tilde{y}}\delta_\rho + \beta_{\mathcal{R}|u}\kappa\beta_y\Delta\tilde{u}_y - \kappa\delta\tilde{E}_y + \kappa\delta\tilde{u}_x = 0$$

$$\partial_{\tilde{t}}\delta\tilde{u}_z + \beta_y\partial_{\tilde{y}}\Delta\tilde{u}_z + \tilde{\beta}_s^2\partial_{\tilde{z}}\delta_n + \beta_y\kappa\delta\tilde{B}_x = 0$$

$$\partial_{\tilde{t}}\Delta\tilde{u}_z + \beta_y\partial_{\tilde{y}}\delta\tilde{u}_z + \tilde{\beta}_s^2\partial_{\tilde{z}}\delta_\rho - \kappa\delta\tilde{E}_z = 0$$

$$\partial_{\tilde{t}}\delta\tilde{B}_x + \partial_{\tilde{y}}\delta\tilde{E}_x - \partial_{\tilde{z}}\delta\tilde{E}_y = 0$$

$$\partial_{\tilde{t}}\delta\tilde{B}_y + \partial_{\tilde{z}}\delta\tilde{E}_x - \partial_{\tilde{x}}\delta\tilde{E}_z = 0$$

$$\partial_{\tilde{t}}\delta\tilde{B}_z + \partial_{\tilde{x}}\delta\tilde{E}_y - \partial_{\tilde{y}}\delta\tilde{E}_x = 0$$

$$\partial_{\tilde{t}}\delta\tilde{E}_x - \partial_{\tilde{y}}\delta\tilde{B}_z + \partial_{\tilde{z}}\delta\tilde{B}_y + \frac{1}{\kappa\sigma}\Delta\tilde{u}_x = 0$$

$$\partial_{\tilde{t}}\delta\tilde{E}_y - \partial_{\tilde{z}}\delta\tilde{B}_x + \partial_{\tilde{x}}\delta\tilde{B}_z + \frac{1}{\sigma}\Delta\tilde{u}_y + \frac{1}{\kappa\sigma}\beta_y\delta_n = 0$$

$$\partial_{\tilde{t}}\delta\tilde{E}_z - \partial_{\tilde{x}}\delta\tilde{B}_y + \partial_{\tilde{y}}\delta\tilde{B}_x + \frac{1}{\kappa\sigma}\Delta\tilde{u}_z = 0. \qquad (B1)$$

This paper has been typeset from a TeX/LATeX file prepared by the author.