# Saturation of the asymmetric current filamentation instability under
conditions relevant to relativistic shock precursors

Virginia Bresci  Institut d’Astrophysique de Paris, CNRS – Sorbonne
Université, 98 bis boulevard Arago, F-75014 Paris  CEA, DAM, DIF, F-91297
Arpajon, France  Laurent Gremillet  CEA, DAM, DIF, F-91297 Arpajon, France
Université Paris-Saclay, CEA, LMCE, 91680 Bruyères-le-Châtel, France  Martin
Lemoine  Institut d’Astrophysique de Paris, CNRS – Sorbonne Université, 98 bis
boulevard Arago, F-75014 Paris

###### Abstract

The current filamentation instability, which generically arises in the
counterstreaming of supersonic plasma flows, is known for its ability to
convert the free energy associated with anisotropic momentum distributions
into kinetic-scale magnetic fields. The saturation of this instability has
been extensively studied in symmetric configurations where the
interpenetrating plasmas share the same properties (velocity, density,
temperature). In many physical settings, however, the most common
configuration is that of asymmetric plasma flows. For instance, the precursor
of relativistic collisionless shock waves involves a hot, dilute beam of
accelerated particles reflected at the shock front and a cold, dense inflowing
background plasma. To determine the appropriate criterion for saturation in
this case, we have performed large-scale 2D particle-in-cell simulations of
counterstreaming electron-positron pair and electron-ion plasmas. We show
that, in interpenetrating pair plasmas, the relevant criterion is that of
magnetic trapping as applied to the component (beam or plasma) that carries
the larger inertia of the two; namely, the instability growth suddenly slows
down once the quiver frequency of those particles equals or exceeds the
instability growth rate. We present theoretical approximations for the
saturation level. These findings remain valid for electron-ion plasmas
provided that electrons and ions are close to equipartition in the plasma flow
of larger inertia. Our results can be directly applied to the physics of
relativistic, weakly magnetized shock waves, but they can also be generalized
to other cases of study.

##  I Introduction

The interpenetration of fast charged particle beams or plasmas gives rise to
the current filamentation instability (CFI), often referred to as the Weibel
instability Weibel (1959); Fried (1959). One of its remarkable features is to
convert part of the free energy associated with the anisotropic momentum
distribution into intense magnetic fields on skin-depth scales, even in the
absence of pre-existing, coherent magnetization Califano  _et al._ (1997);
Bret  _et al._ (2004, 2005, 2007, 2008, 2010a, 2010b); Wiersma and Achterberg
(2004); Achterberg and Wiersma (2007); Achterberg  _et al._ (2007). Its
generic nature, its robustness and its physical implications have thus given
it a compelling role in many fields of research, e.g., high-energy density
physics Lee and Lampe (1973); Molvig (1975); Silva  _et al._ (2002); Adam  _et
al._ (2006); Debayle  _et al._ (2010), laboratory astrophysics Fiuza  _et al._
(2012); Huntington  _et al._ (2015); Lobet  _et al._ (2015); Warwick  _et al._
(2017); Li  _et al._ (2019); Fiuza  _et al._ (2020), cosmology Gruzinov
(2001); Schlickeiser (2005), and high-energy astrophysics, where it is thought
to shape much of the nonthermal electromagnetic radiation from powerful
explosive transients Medvedev and Loeb (1999); Gruzinov and Waxman (1999);
Silva  _et al._ (2003); Frederiksen  _et al._ (2004); Kirk and Reville (2010);
Bykov and Treumann (2011); Medvedev  _et al._ (2011); Plotnikov  _et al._
(2013); Lemoine (2013); Kelner  _et al._ (2013); Spitkovsky (2008a, b).

Specifically, the CFI regulates the structure of weakly magnetized,
collisionless shocks, which form through the counterstreaming of plasma shells
at supersonic speeds Moiseev and Sagdeev (1963). As such, it controls the
production of high-energy particles and radiation in such environments,
whether in the near Sulaiman  _et al._ (2015) or the remote Universe Medvedev
and Loeb (1999). Because of their paramount consequences, those “Weibel-
mediated” collisionless shocks have inspired a large body of literature, both
in the subrelativistic Kato and Takabe (2008, 2010); Bret  _et al._ (2013);
Matsumoto  _et al._ (2015); Bohdan  _et al._ (2020a, b); Ruyer  _et al._
(2015a, b, 2016, 2017) and relativistic regimes Lyubarsky and Eichler (2006);
Lemoine  _et al._ (2006); Kato (2007); Spitkovsky (2008a); Martins  _et al._
(2009); Nishikawa  _et al._ (2009); Keshet  _et al._ (2009); Sironi  _et al._
(2013); Plotnikov  _et al._ (2018); Tomita  _et al._ (2019); Lemoine  _et al._
(2019a, b, c); Pelletier  _et al._ (2019); Vanthieghem  _et al._ (2020). Long-
standing issues, with obvious phenomenological implications in the
aforementioned domains of research, are the level of saturation of the CFI
Davidson  _et al._ (1972); Kato (2005); Califano  _et al._ (1998); Okada and
Ogawa (2007); Shvets  _et al._ (2009); Cagas  _et al._ (2017); Grassi  _et
al._ (2017); Achterberg  _et al._ (2007); Takamoto  _et al._ (2018, 2019) and
the long-term evolution of the self-generated magnetic turbulence Jaroschek
_et al._ (2004); Medvedev  _et al._ (2004); Milosavljevic and Nakar (2006);
Chang  _et al._ (2008); Lemoine (2015); Ruyer  _et al._ (2015a); Vanthieghem
_et al._ (2018). Saturation is thought to occur through either transverse
trapping of the particles in the magnetic filaments Davidson  _et al._ (1972);
Yang  _et al._ (1994); Istomin (2011), cyclotron gyration of the particles
around the magnetic-field extrema Moiseev and Sagdeev (1963); Califano  _et
al._ (1998); Kato (2005); Istomin (2011), or exhaustion of the available
particle current Kato (2005). Most previous studies on the saturation of the
CFI, and, to our knowledge, all those related to relativistic astrophysical
systems, have considered symmetric configurations in which the
interpenetrating plasmas share similar characteristics (i.e., identical
temperatures, densities and drift velocities). In the precursor of
collisionless shock waves, however, the interaction is strongly asymmetric, as
it involves a hot, dilute beam of accelerated particles interacting with a
cold, dense background plasma.

The main objective of the present work is therefore to examine the saturation
processes of the CFI in generic asymmetric configurations. In particular, we
seek to determine which criterion holds, and whether this criterion applies to
the beam, or to the background plasma, an ambiguity which obviously does not
arise in symmetric configurations. We do so by confronting analytical
predictions with particle-in-cell (PIC) simulations of initially unmagnetized,
collisionless plasmas. For most of our study, we consider plasmas composed of
equal mass species, interpenetrating each other at a relativistic velocity.
This configuration is typical of the precursor of a relativistic shock
propagating in a pair plasma, but it is also relevant for the study of the CFI
in asymmetric electron-electron or ion-ion flows. We nonetheless extend our
simulations to the case of electron-ion plasmas, in the ultrarelativistic and
mildly relativistic regimes.

Our study is laid out as follows. In Section II, we recall the salient
features of the CFI in asymmetric counterstreaming flows. In particular, we
emphasize the notion of the preferred “Weibel frame”, in which the instability
is of a purely magnetic nature, and which becomes crucial in the asymmetric
interaction regime. We then discuss the main saturation mechanisms and give
the corresponding estimates of the maximum magnetic field energy. In Sec. III,
we present our PIC simulations for pair plasmas and analyze their results in
light of the above mechanisms. We extend our analysis to the case of electron-
ion counterstreaming configurations in Sec. IV, and finally summarize our
results and conclusions in Sec. VI.

##  II The asymmetric current filamentation instability

In this section as well as the next one, our initial setup comprises two
counterstreaming, unmagnetized pair plasmas drifting along the $x-$axis. We
note the beam with a subscript b and the plasma with a subscript p. By
convention, the beam corresponds to the population with the lower relativistic
plasma frequency. The latter is defined as (cgs units are used throughout)

$\Omega_{\rm p\alpha}=\left(\frac{4\pi
n_{\alpha}e^{2}}{\widetilde{w}_{\alpha}/c^{2}}\right)^{1/2}\,,$ (1)

where $e$ is the elementary charge, $c$ the velocity of light, $n_{\alpha}$
the proper number density and $\widetilde{w}_{\alpha}$ the enthalpy per
particle of charged species $\alpha\in\\{b+,b-,p+,p-\\}$ in its initial state.
Note that, in our notations, $n$ refers to a single charged species; it thus
represents half of the initial total number density of the corresponding
component (beam or plasma). Introducing the corresponding particle mass
$m_{\alpha}$, adiabatic index $\widehat{\Gamma}_{\alpha}$ and proper
temperature $T_{\alpha}$, one has
$\widetilde{w}_{\alpha}=m_{\alpha}c^{2}+\widehat{\Gamma}_{\alpha}k_{\rm
B}T_{\alpha}/(\widehat{\Gamma}_{\alpha}-1)$. This implies
$\widetilde{w}_{\alpha}\simeq m_{\alpha}c^{2}$ for a plasma of subrelativistic
temperature ($k_{\rm B}T_{\alpha}/m_{\alpha}c^{2}\ll 1$), and
$\widetilde{w}_{\alpha}\simeq\widehat{\Gamma}_{\alpha}k_{\rm
B}T_{\alpha}/(\widehat{\Gamma}_{\alpha}-1)$ for a relativistically hot plasma
($k_{\rm B}T_{\alpha}/m_{\alpha}c^{2}\gg 1$). Given the inverse normalized
temperature of species $\alpha$, $\mu_{\alpha}\equiv
m_{\alpha}c^{2}/(k_{B}T_{\alpha})$, one has $\widetilde{w}_{\alpha}\simeq
m_{\alpha}c^{2}$ and $\Omega_{\rm p\alpha}\simeq\omega_{\rm p\alpha}$ for a
plasma of nonrelativistic temperature ($\mu_{\alpha}\gg 1$), but
$\widetilde{w}_{\alpha}\simeq\widehat{\Gamma}_{\alpha}k_{\rm
B}T_{\alpha}/(\widehat{\Gamma}_{\alpha}-1)$, and hence $\Omega_{\rm
p\alpha}\simeq\omega_{\rm p\alpha}\sqrt{\mu_{\alpha}}/2$ (taking
$\widehat{\Gamma}_{\alpha}=4/3$) for a relativistically hot plasma
($\mu_{\alpha}\ll 1$). In the following, use will also be made of

$\omega_{\rm p\alpha}=\left(4\pi n_{\alpha}e^{2}/m_{\alpha}\right)^{1/2}\,,$
(2)

the nonrelativistic plasma frequency of species $\alpha$.

In this configuration, the counterstreaming instability can be described in
wavenumber space $(k_{\parallel},\,\bm{k_{\perp}})$, in terms of the
longitudinal $k_{\parallel}=\bm{k}\cdot\bm{\hat{x}}$ and perpendicular
$\bm{k_{\perp}}$ wavenumbers. This instability breaks into two main branches
Bret  _et al._ (2010a): purely transverse modes (the CFI) with
$k_{\parallel}\ll k_{\perp}$, and the so-called oblique two-stream modes, for
which $k_{\parallel}\sim k_{\perp}$. We neglect here the purely parallel
electrostatic branch with $k_{\perp}\ll k_{\parallel}$, which is usually
subdominant in the relativistic limit. The oblique modes are essentially
electrostatic, while the transverse CFI modes are essentially magnetic Bret
_et al._ (2010c), of direct interest to the present study.

###  II.1 The “Weibel frame”

The notion of being “magnetic” or “electrostatic” is a frame-dependent
statement, which must be made precise in a relativistic setting. This is
discussed in detail in Ref. Pelletier  _et al._ (2019), and we recap here the
most salient features. The purely transverse CFI (meaning
$k_{\parallel}\rightarrow 0$) develops through the pinching of the
counterstreaming plasmas into filamentary structures oriented along
$\bm{\hat{x}}$, each endowed with a net current. These structures are
surrounded by toroidal magnetic fields $\bm{\delta B_{\perp}}$ and radial
electric fields $\bm{\delta E_{\perp}}$. The CFI growth also comes with an
inductive electric field component $\bm{\delta E_{\parallel}}$, oriented along
the drift direction. In a first approximation, the latter field can be
neglected, because its magnitude is of the order of $|\Im\omega/k_{\perp}c|\ll
1$ relative to the magnetic field component. The dominance of the magnetic
component means $\delta B_{\perp}^{2}-\delta E_{\perp}^{2}>0$ for each
unstable wavenumber $\mathbf{k_{\perp}}$, and hence that there exists a frame,
moving at velocity (in units of $c$)

$\bm{\beta_{\rm w}}=\frac{\bm{\delta E_{\perp}}\times\bm{\delta
B_{\perp}}}{\delta B_{\perp}^{2}}\,,$ (3)

in which the transverse electric field component vanishes. In this frame,
which we call the “Weibel frame”, the CFI can be regarded as purely magnetic,
up to the weak inductive component which cannot be erased by a Lorentz boost.

In the precursor of relativistic shocks, this frame gains special importance
because the interaction between the beam of accelerated particles and the
background plasma is so asymmetric that $\beta_{\rm w}\simeq 1$, meaning
$\delta E_{\perp}\simeq\delta B_{\perp}$. It is crucial to properly
characterize this frame, as it controls the heating and slowdown of the
background plasma and because it greatly helps evaluate the scattering rate of
suprathermal particles Lemoine  _et al._ (2019a). Hence, the Weibel frame is
connected to acceleration processes and has direct phenomenological
consequences.

Let us consider a set of initial beam $\left\\{n_{\rm b},\,T_{\rm b},\,u_{\rm
b|r}\right\\}$ and plasma $\left\\{n_{\rm p},\,T_{\rm p},\,u_{\rm
p|r}\right\\}$ parameters, where
$u_{\rm\alpha|r}\equiv\gamma_{\rm\alpha|r}\beta_{\rm\alpha|r}$ denotes the
$x-$component of the four-velocity of species $\alpha$, and
$\beta_{\rm\alpha|r}$ and $\gamma_{\rm\alpha|r}$ are the associated normalized
three-velocity (in units of $c$) and Lorentz factor, all defined in some
reference frame (subscript |r). We use proper densities and temperatures
unless explicitly specified otherwise. The corresponding Weibel frame velocity
can be determined in the following two ways.

In the linear phase of the CFI, one can define the Weibel frame as that in
which the electrostatic component of the dispersion relation vanishes. This
has been done in Ref. Ruyer  _et al._ (2016) in the subrelativistic regime,
and in Refs. Pelletier  _et al._ (2019); Lemoine  _et al._ (2019a) in the
relativistic regime. This amounts to setting the $\epsilon_{xy}$ component of
the total dielectric tensor to zero, assuming
$\bm{k_{\perp}}=k_{\perp}\bm{\hat{y}}$. This is not a trivial step, as the
dielectric tensor itself depends on the solution to the dispersion relation,
see Ref. Pelletier  _et al._ (2019) for a discussion of the procedure.

Alternatively, one can describe the nonlinear phase of the instability as a
quasistatic equilibrium between particles and fields, ordered along the
transverse $y-$direction (a reduced 2D $x-y$ geometry is assumed throughout
for simplicity) in a periodic sequence of current filaments. In a four-fluid
(isothermal) description, the density of each component at equilibrium can be
written as a function of the electromagnetic potentials, see Ref. Vanthieghem
_et al._ (2018) for details. Setting the electrostatic contribution to zero
imposes a relationship between the physical characteristics of the fluid, in
the form

$\frac{n_{\rm b}\gamma_{\rm b|w}^{2}\beta_{\rm b|w}}{T_{\rm b}}+\frac{n_{\rm
p}\gamma_{\rm p|w}^{2}\beta_{\rm p|w}}{T_{\rm p}}=0\,,$ (4)

where the normalized $x$-velocities $\beta_{\alpha|\rm w}$ and Lorentz factors
$\gamma_{\alpha|\rm w}$ are here measured in the Weibel frame. The above
equation can be solved to obtain the velocity of the Weibel frame in the
reference frame. As it turns out, both methods give similar expressions for
this velocity under conditions relevant to the precursor of relativistic
shocks. Here, we rely on the latter method and make the result explicit, as
follows.

Writing $\beta_{\alpha|\rm w}$ and $\gamma_{\alpha|\rm w}$ in terms of
$\beta_{\rm\alpha|r}$ and $\gamma_{\rm\alpha|r}$ through standard Lorentz
transforms, one finds that the Weibel frame velocity, relative to the
reference frame, can be expressed as

$\beta_{\rm w|r}=\frac{Q_{\rm w}-\sqrt{Q_{\rm w}^{2}-4}}{2}\,,$ (5)

where

$Q_{\rm w}=\frac{n_{\rm b}\gamma_{\rm b|r}^{2}\left(1+\beta_{\rm
b|r}^{2}\right)/T_{\rm b}+n_{\rm p}\gamma_{\rm p|r}^{2}\left(1+\beta_{\rm
p|r}^{2}\right)/T_{\rm p}}{n_{\rm b}\gamma_{\rm b|r}^{2}\beta_{\rm b|r}/T_{\rm
b}+n_{\rm p}\gamma_{\rm p|r}^{2}\beta_{\rm p|r}/T_{\rm p}}\,.$ (6)

The minus sign in Eq. (5) reflects the fact that $\beta_{\rm
w|r}\simeq\beta_{\rm p|r}$ if the beam component becomes negligible: the
turbulence is then mostly magnetic in the rest frame of the background plasma.
Once $\beta_{\rm w|r}$ is known, the velocity of each species in a given
reference frame can be Lorentz transformed to the Weibel frame.

Henceforth, all velocities or Lorentz factors that do not carry a subscript |r
are understood to be defined in the Weibel frame. The Weibel frame associated
with the _initial_ state of the system is the reference frame in which our
simulations will be conducted. Note that this frame can differ from the
_instantaneous_ Weibel frame that results from the time-evolving properties of
the beam and the plasma as the instability develops. This will be manifest in
our simulations, and we will return to this point in Sec. III.

###  II.2 Linear stage of the CFI growth

Let us first consider the linear properties of the purely transverse CFI
modes, that is, with wave vector $\mathbf{k_{\perp}}=k_{\perp}\bm{\hat{y}}$
and frequency $\omega\equiv i\Gamma_{\rm w}$, where $\Gamma_{\rm w}$ is the
$k_{\perp}$-dependent growth rate. In the linear phase of the CFI, each mode
grows as

$\delta B_{z}(k_{\perp})=\delta B_{0}(k_{\perp})e^{\Gamma_{\rm
w}(k_{\perp})t}\,,$ (7)

where $\delta B_{0}(k_{\perp})$ is the seed magnetic field fluctuation.
Assuming that the magnetic spectrum ends up being dominated by modes of
similar growth rate and seeded by comparable fluctuations, one can infer the
instantaneous growth rate through

$\Gamma_{\rm w}=\frac{1}{2}\frac{{\rm d}}{{\rm
d}t}\ln\left[\frac{\langle\delta B_{z}(t)^{2}\rangle}{\langle\delta
B_{z}(0)^{2}\rangle}\right]\,.$ (8)

The quantity in the rhs can be easily extracted from numerical simulations and
directly compared with analytic estimations of $\Gamma_{\rm w}$. The latter
involve rather heavy calculations of the dielectric tensor contained in the
kinetic dispersion relation of which we will summarize here only the general
key points.

The system is initially charge and current neutral with no equilibrium
electromagnetic fields. Linearizing the Vlasov-Maxwell equations by expressing
every perturbed physical quantity as $\delta\xi\propto
e^{i(\bm{k}_{\perp}\cdot\bm{r}-\omega t)}$ yields the dispersion relation of
the CFI:

$\epsilon_{yy}\left(\epsilon_{xx}-1/\zeta^{2}\right)=\epsilon_{xy}^{2}\,,$ (9)

where $\zeta=\omega/k_{\perp}c$, and the elements of the dielectric tensor are
given by

$\displaystyle\epsilon_{ij}$
$\displaystyle=\delta_{ij}+\sum_{\alpha}\frac{\gamma_{\alpha}\omega^{2}_{p\alpha}}{\zeta^{2}k_{\perp}^{2}c^{2}}\,\int\frac{u_{i}}{\gamma}\frac{\partial
f_{\alpha}^{(0)}}{\partial u^{j}}\,d^{3}u$ (10)
$\displaystyle+\sum_{\alpha}\frac{\gamma_{\alpha}\omega^{2}_{p\alpha}}{\zeta^{2}k_{\perp}^{2}c^{2}}\,\int\frac{u_{i}u_{j}}{\gamma^{2}}\frac{\partial
f_{\alpha}^{(0)}/\partial u_{y}}{\zeta-\beta_{y}}\,d^{3}u\,,$ (11)

where $i,j=(1,2,3)$ and $\mathbf{u}=\gamma\bm{\beta}$. Hence, the
filamentation instability is generally not purely magnetic unless the off-
diagonal term of the dielectric tensor vanishes, which would be the case for
symmetric counterstreaming flows.

In the present work, we consider particle populations characterized by
Maxwell-Jüttner momentum distribution functions. As shown in Ref. Pelletier
_et al._ (2019), approximate growth rates of the CFI can be obtained in two
asymptotic limits that depend on the value of the parameter
$\chi_{\alpha}=\gamma_{\alpha}|\zeta|/\sqrt{1-\zeta^{2}}$. For each plasma
species, we define the hydrodynamic limit in which the thermal velocity spread
of the distribution function is, broadly speaking, smaller than the
(imaginary) “phase velocity” of the waves, and the opposite _kinetic_ limit.
More precisely, recalling that $\mu_{\alpha}=m_{\alpha}c^{2}/k_{\rm
B}T_{\alpha}$, the hydrodynamic (resp. kinetic) limit for the cold plasma
component corresponds to $\tilde{\chi}_{\rm p}\equiv\chi_{\rm p}\sqrt{\mu_{\rm
p}/2}\gg 1$ (resp. $\ll 1$). For the relativistically hot beam component, the
hydrodynamic (resp. kinetic) limit is rather defined as $\chi_{\rm b}\gg 1$
(resp. $\ll 1$), see Pelletier  _et al._ (2019) for details.

We can thus derive two useful approximations of the maximum growth rate and
associated wave number in terms of the nonrelativistic plasma frequencies of
the plasma species, one in the fully kinetic regime – meaning the kinetic
approximation for both species – and one in the combined hydrodynamical (beam)
and kinetic (plasma) regimes, respectively,

$\displaystyle\Gamma_{\rm w,k-k}$ $\displaystyle\simeq\frac{(\omega_{\rm
pb}^{2}\mu_{\rm b})^{3/2}\gamma_{\rm b|p}^{3}\beta_{\rm
b|p}^{3}}{\sqrt{2\pi\mu_{\rm p}}\omega_{\rm pp}^{2}+\frac{3\pi}{2}\omega_{\rm
pb}^{2}\mu_{\rm b}\gamma_{\rm b|p}^{3}}\,,$ $\displaystyle k_{\rm\perp,k-k}$
$\displaystyle\simeq\sqrt{\frac{2}{3}\mu_{\rm b}}\gamma_{\rm b|p}\omega_{\rm
pb}\,,$ (12)

and

$\displaystyle\Gamma_{\rm w,k-h}$ $\displaystyle\simeq\sqrt{\omega_{\rm
pb}^{2}\mu_{\rm b}}\,,$ $\displaystyle k_{\rm\perp,k-h}$
$\displaystyle\simeq(2\pi\omega_{\rm pb}^{2}\mu_{\rm p}\mu_{\rm
b})^{1/6}\omega_{\rm pp}^{2/3}\,.$ (13)

The quantities $\beta_{\rm b|p}$ and $\gamma_{\rm b|p}$ represent the
normalized three-velocity of the beam relative to the plasma and its
corresponding Lorentz factor. It is important to stress that the above
formulae have been derived solving the dispersion relation of the instability
making the approximation of cold plasma ($k_{\rm B}T_{\rm p}\ll
m_{\alpha}c^{2}$) and hot beam ($k_{\rm B}T_{\rm b}\gg m_{\alpha}c^{2}$) in
the respective dielectric tensor. Furthermore, those approximations assume
that the plasma moves at subrelativistic velocities with respect to the Weibel
frame; consequently, it neglects terms of order $\mathcal{O}(\beta_{\rm
p|w})$. Those formulas encompass the majority of the situations addressed in
the following but not all; this will be made explicit.

###  II.3 Saturation criteria for the CFI

In the early linear stage of the instability, the charged particles are
deflected by magnetic field fluctuations with a polarity perpendicular to
their initial drift velocity. As a result, particles of opposite charge from
each component of the system (beam or plasma) are focused in different
regions, forming transverse current modulations or “current filaments”.
Particles of opposite charges from both components concentrate in the same
filaments where their currents add up; this amplifies the initial magnetic
field perturbation, thus leading to the development of the instability.

Eventually, the particle dynamics becomes modified by the fields so that
saturation mechanisms take place. Ultimately, the CFI enters a strongly
nonlinear stage, in which secondary instabilities, such as the merging of
filaments of equal polarity, or the kink of current filaments, can arise, see
Ref. Vanthieghem  _et al._ (2018) for a detailed discussion. The transition
between these two phases, i.e., saturation and the strongly nonlinear stage,
is fraud with ambiguities, as filaments can coalesce while the current
filaments keep building up through the CFI. We define here the saturation as
the point at which the growth of the magnetic energy density is halted, or at
least significantly reduced. This will be made clear in the figures that
follow.

To investigate the saturation of the instability, we will compare the temporal
evolution of the magnetic field as extracted from simulations, with different
criteria of saturation borrowed from the literature, which we summarize below.
We emphasize a key difference with respect to the case of symmetric
counterstreaming plasmas, which are more commonly envisaged. In the asymmetric
configuration, an ambiguity arises as to which component (beam or plasma) is
eventually responsible for the saturation, and through which mechanism. For
this reason, we discuss in the forthcoming paragraphs the saturation criteria
as applied to a generic component. We will then apply each of them to the beam
and to the plasma and compare those to the simulation results in the next
Section.

####  II.3.1 Transverse trapping

The widely used trapping-based saturation criterion, first proposed by
Davidson in the nonrelativistic regime Davidson  _et al._ (1972), and later
generalized to the relativistic regime Yang  _et al._ (1994); Lyubarsky and
Eichler (2006); Achterberg  _et al._ (2007); Kaang  _et al._ (2009), expresses
the fact that, in the weakly nonlinear phase of the CFI, particles quiver
transversely around the center of the filament (i.e., around a magnetic field
node) in which they are focused. Assuming a harmonic $B$-field profile of
amplitude $B$ and wavenumber $k_{\perp}$, a particle of Lorentz factor
$\gamma$ and mass $m$ oscillates at the bounce frequency

$\omega_{\rm B}=\left(\frac{ek_{\perp}\beta_{\parallel}B}{\gamma
m}\right)^{1/2}\,.$ (14)

The onset of saturation can be viewed as when the assumption of zero-order
ballistic particle motion no longer holds. This occurs when $\omega_{\rm B}$
becomes comparable with the instability growth rate, $\Gamma_{\rm w}$.
Introducing $\langle\gamma\rangle$ the typical Lorentz factor of the
considered species, the corresponding saturation magnetic field can thus be
expressed as

$B_{\rm t}=\frac{\Gamma_{\rm w}^{2}}{k_{\perp}}\frac{\langle\gamma\rangle
m}{\beta_{\parallel}e}\,.$ (15)

####  II.3.2 Magnetization limit

In the nonlinear phase of the CFI, the plasma can be modelled as an ensemble
of cylindrical filaments of radius
$r\simeq\lambda_{\perp}/4\simeq\pi/2k_{\perp}$, carrying a current density
$j$. As the $B$-field grows in amplitude, the Larmor radius of the particles,
$r_{\rm L}=\gamma\beta mc^{2}/eB$, shrinks, possibly up to the point where it
becomes smaller than the filament radius. Particles then become _spatially
trapped_ within the filaments in both the longitudinal and transverse
directions, while orbiting around the $B$-field extrema. In the literature,
this limit is often referred to as the “Alfvén limit” Kato (2005). Similarly,
particles gyrating at a Larmor frequency $\omega_{L}=eB/m\gamma$ higher than
the instability growth rate can be regarded as _temporally magnetized_. In
either case, the linear approximation, which assumes rectilinear motion across
the filaments, breaks down. The maximum value of the magnetic field set by
this condition is then given by

$B_{\rm
m}=\max{(B_{\mathrm{m},\,r_{\mathrm{L}}},B_{\mathrm{m},\,\omega_{\mathrm{L}}})}\,,$
(16)

where

$B_{\mathrm{m},\,r_{\mathrm{L}}}=\frac{2}{\pi}k_{\perp}\langle\gamma\beta\rangle\frac{mc^{2}}{e}$
(17)

satisfies the spatial constraint and

$B_{\mathrm{m},\,\omega_{\mathrm{L}}}=\Gamma_{\rm
w}\langle\gamma\rangle\frac{mc}{e}$ (18)

the temporal one. Since the CFI is characterized by $\Gamma_{\rm w}\ll
k_{\perp}c$ in relativistic shock precursors Pelletier  _et al._ (2019), it
follows that usually $B_{\rm m}=B_{\mathrm{m},\,r_{\mathrm{L}}}$ if $\beta\sim
1$. Similar saturation criteria were considered in Moiseev and Sagdeev (1963);
Medvedev and Loeb (1999); Lyubarsky and Eichler (2006); Bret  _et al._ (2013).

####  II.3.3 Particle current limit

The magnetic field is also bounded by above by the maximum current density
that can sustain it Kato (2005). This maximum current density corresponds to
the current carried by one of the two oppositely charged species making up the
component (i.e., beam or plasma) under study. This limit thus tacitly assumes
that, at maximum magnetic field, all the particles of a given component within
a transverse length $\lambda_{\perp}/2$ have undergone complete spatial
separation in two adjacent filaments. Assuming these have a uniform current
density, the $B$-field created by a charged species of initial apparent
density $\gamma n$ (with $\gamma$ characterizing here the drift motion) has a
maximum strength

$B_{p}\simeq 2\pi^{2}\frac{e\gamma
n}{k_{\perp}}\langle\beta_{\parallel}\rangle\,.$ (19)

In the case of complete spatial separation, the contributions of
counterstreaming species of opposite charge should add up within a filament.
Yet in the asymmetric configurations addressed in the following, only the
particle limit associated with the component that carries most, if not all of
the particle current density, matters. We will therefore identify that
component, which will turn out to be the beam component in most cases, and
ignore the particle limit associated with the other (background plasma)
component.

###  II.4 Analytical estimates

The hierarchy among the above saturation criteria depends on the
characteristic wave number of the instability and the growth rate, given that

$\displaystyle\frac{B_{\rm t}}{B_{\rm p}}$
$\displaystyle\sim\bigg{(}\frac{\Gamma_{\rm w}}{\omega_{\rm
p}}\bigg{)}^{2}\,,$ (20) $\displaystyle\frac{B_{\rm m}}{B_{\rm p}}$
$\displaystyle\sim\left(\frac{k_{\perp}c}{\omega_{\rm p}}\right)^{2}\,,$ (21)

where $\omega_{\rm p}$ represents here the nonrelativistic plasma frequency of
the component to which the saturation criterion is applied, and $k_{\perp}$
denotes the dominant transverse wavenumber. Considering first a cold symmetric
counterstreaming configuration, one has $\Gamma_{\rm w}\sim\omega_{\rm pp}$
and $k_{\perp}\gg\omega_{\rm pp}/c$ to leading order, _e.g._ Wiersma and
Achterberg (2004); Bret  _et al._ (2008). As a consequence, $B_{\rm t}\simeq
B_{\rm p}\ll B_{\rm m}$, implying that the trapping and particle limits are
equivalent and determine saturation. For symmetric counterstreaming hot
plasmas, $\Gamma_{\rm w}$ is reduced to values below $\omega_{\rm pp}$,
because it scales with the relativistic plasma frequency $\Omega_{\rm
p}=\omega_{\rm p}\sqrt{\mu}/2$ and $\mu\ll 1$. Consequently, the trapping
criterion is expected to become more stringent than the other two. In
addition, a relativistic temperature likely prevents the oppositely charged
species of a given component from fully segregating from each other within a
filament, further weakening the particle limit in this regime.

In an asymmetric configuration, we identify distinct saturation limits for the
beam and the plasma, using the respective superscripts b and p. As a general
trait of such configurations, we observe that the beam moves at relativistic
velocities in the Weibel frame, while the drift of the background plasma is
most often sub- or mildly relativistic. This can be read off Eq. (4), which
relates the quantities $n_{\rm b}/T_{\rm b}\propto\Omega_{\rm pb}^{2}$ and
$n_{\rm p}/T_{\rm p}\propto\omega_{\rm pp}^{2}\mu_{\rm p}$. The beam is
usually defined as the component with the smaller plasma frequency of the two,
hence $\Omega_{\rm pb}\ll\omega_{\rm pp}$ suggests that $\gamma_{\rm b|\rm
w}^{2}|\beta_{\rm b|w}|\gg\gamma_{\rm p|\rm w}^{2}|\beta_{\rm p|w}|$.
Therefore, one must expect $B_{\rm t}^{\rm b}\gg B_{\rm t}^{\rm p}$. We
anticipate that, for what concerns saturation through trapping, only the
larger of the two values $B_{\rm t}^{\rm b}$ and $B_{\rm t}^{\rm p}$ matters,
and this trend will be confirmed by the simulations.

We also expect, for the same reasons as above, that $B_{\rm t}^{\rm b}<B_{\rm
p}^{\rm b}$ and $B_{\rm t}^{\rm b}\ll B_{\rm m}^{\rm b}$, because of the large
temperature of the beam. Consequently, we may anticipate that the overall
criterion for saturation will be set by the trapping limit of beam particles.

##  III PIC simulations results

We have performed a number of 2D3V (2D in space, 3D in momentum) PIC
simulations of counterstreaming electron-positron pair plasmas, which
initially obey Maxwell-Jüttner distribution functions, using the massively
parallel calder code Lefebvre  _et al._ (2003). To resolve properly the
initial Weibel instability, the cell size is set to $\Delta x=\Delta
y=0.1\,c/\omega_{\rm pp}$ and the simulations are run over $2\times 10^{4}$
time steps of $\Delta t=0.099\,\omega_{\rm pp}^{-1}$ on a 2D $(x,y)$ grid of
$2000\times 2000$ cells. Henceforth, $\omega_{\rm pp}$ represents the
nonrelativistic plasma frequency of each of the two charged species of the
plasma component in its initial state, i.e., $\omega_{\rm pp}=\left(4\pi
n_{\rm p}e^{2}/m_{e}\right)^{1/2}$ ($m_{e}$ is the electron mass). Each cell
contains initially $100$ macro-particles per species, yielding a total number
of about $10^{9}$ macro-particles. Time and length are normalized to the
inverse nonrelativistic plasma frequency $\omega_{\rm pp}^{-1}$ and the plasma
inertial length $c/\omega_{\rm pp}$.

As previously mentioned, we aim to investigate the saturation of the current
filament instability in an asymmetric interaction between a hot dilute beam
and a cold, dense, inflowing plasma as it happens in the precursor of
astrophysical collisionless shock waves in pair plasmas. For this reason, we
initiate our study making use of initial parameters borrowed from a large-
scale shock simulation corresponding to a relative upstream to downstream
Lorentz factor of $10$, as described in Lemoine  _et al._ (2019a). The
parameters of the beam and plasma populations, as measured in the downstream
shock frame, are as follows (temperatures are given in units of
$m_{e}c^{2}/k_{\rm B}$) : $\gamma_{\rm b|d}=1.38$, $\gamma_{\rm p|d}=9.67$,
$T_{\rm b}=45$, $T_{\rm p}=0.2$, and $n_{\rm b}/n_{\rm p}=0.1$. Those values
are extracted from a region deep inside the precursor of the shock, where the
background plasma has been slightly slowed down and heated to mildly
relativistic temperatures. As announced, we then transform those initial
parameters from the downstream shock frame to the Weibel frame. This gives the
set of parameters indicated by (a) in Table 1, hereafter referred to as the
reference run. Note that the plasma moves at subrelativistic velocities in
this Weibel frame, while the beam is now ultrarelativistic. This difference
demonstrates the importance of the Weibel frame regarding the development of
the instability, and more importantly, regarding its saturation, since the
saturation criteria depend on the inertia of the particles, which in turn
depend on the reference frame.

The parameters of subsequent runs have been varied accordingly to fall in the
region of the parameter space dominated by the CFI over electrostatic and
oblique modes. In particular, we investigate a case where the initial beam
proper density is tripled with respect to the reference case [run (b)], one in
which the initial beam proper temperature is reduced by a factor of $1/3$ [run
(c)], one with an initial beam Lorentz factor reduced by a factor of $1/3$
[run (d)]. Finally, we examine two more extreme configurations by reducing the
initial temperature of the beam while increasing its initial Lorentz factor
and the initial plasma temperature by a factor of $10$ each [run (e)] or $30$
each [run (f)]. The latter runs are of particular interest for the present
study, because their parameters are such that the roles of background plasma
and beam are interchanged with respect to other runs.

What we refer to as the beam is set in motion in the positive $x-$direction
and represents the hot cloud reflected by the shock, which encounters the cold
incoming plasma streaming along the negative direction. Correspondingly, the
transverse CFI generates an out-of-plane magnetic field component, $B_{z}$,
aligned with the $\bm{\hat{z}}$ direction, and its associated electrostatic
component $E_{y}$, along $\bm{\hat{y}}$. Since the simulation frame initially
coincides with the Weibel frame, $E_{y}$ remains much smaller than $B_{z}$
during the initial development of the instability. A stronger $E_{y}$ then
emerges gradually, and as time progresses, the physical conditions of the
plasma and/or the beam change, and so does the instantaneous Weibel frame. In
particular, the filamentary structures start to move along $\bm{\hat{x}}$ at
an approximately coherent velocity corresponding to the time-dependent value
of $\beta_{\rm w}$. To discriminate between the various saturation criteria,
the magnetic field is directly extracted from PIC simulations and compared
with the theoretical estimates of the saturated $B$-field given in Sec. II.3.

In what follows, we focus on the linear and saturation phases of the
instability, while the late-time evolution is left aside and treated in Sec.
V.

Run |  $\gamma_{\rm b}$ |  $\gamma_{\rm p}$ |  $T_{\rm b}/\gamma_{\rm b}$ |  $T_{\rm p}/\gamma_{\rm p}$ |  $\gamma_{\rm b}n_{\rm b}/\gamma_{\rm p}n_{\rm p}$
---|---|---|---|---|---
(a) |  $18.9$ |  $1.01$ |  $2.4$ |  $0.2$ |  $1.9$
(b) |  $16.7$ |  $1.05$ |  $2.9$ |  $0.2$ |  $4.6$
(c) |  $16.7$ |  $1.05$ |  $0.95$ |  $0.2$ |  $1.5$
(d) |  $40.8$ |  $1.2$ |  $1.09$ |  $0.2$ |  $3.5$
(e) |  $25.$ |  $5.3$ |  $0.2$ |  $0.4$ |  $0.5$
(f) |  $25.$ |  $15.4$ |  $0.06$ |  $0.4$ |  $0.2$
Table 1:  Summary of simulation parameters for pair plasmas. Run (a) is the
reference simulation. The parameters of the other runs differ from those of
run (a) as follows: (b) $n_{\rm b}\times 3$; (c) $T_{\rm b}/3$; (d)
$\gamma_{\rm b|d}\times 3$; (e) $\gamma_{\rm b|d}/3$,(e) $\gamma_{\rm
b|d}\times 10$, $T_{\rm p}\times 10$, $T_{\rm b}/10$; (f) $\gamma_{\rm
b|d}\times 30$, $T_{\rm p}\times 30$, $T_{\rm b}/30$. The table gives the
simulation parameters once transformed to the Weibel frame. Temperatures are
given in units of $m_{e}c^{2}/k_{\rm B}$.

###  III.1 Reference run

Figure 1: Temporal evolution of the simulated mean $B$-field strength
($B_{z}$, black curves) compared to various saturation criteria for reference
run (a). Top panel: particle ($B_{\rm p}^{\rm b}$, green dashed-dotted curve)
and trapping ($B_{\rm t}^{\rm b}$, blue dotted curve) limits as applied to the
beam particles. Bottom panel: spatial magnetization ($B_{\rm m}^{\rm p}$, red
dashed curve) and trapping ($B_{\rm t}^{\rm p}$, blue dotted curve) limits as
applied to the plasma particles. All curves are in units of $m_{e}c\omega_{\rm
pp}/e$.

The growth of the magnetic field during the linear and saturation phases of
our reference case (a) can be clearly seen in Fig. 1 (thick black line). In
this figure, and subsequent similar ones, the $B$-field strength is expressed
in dimensionless units, $\overline{B}_{z}=eB_{z}/m_{e}c\omega_{\rm
pp}=B_{z}/\sqrt{4\pi n_{p}m_{e}c^{2}}$. The expected maximum growth rate is
$\Gamma_{\rm w}\simeq 0.02\,\omega_{\rm pp}$ at $k_{\perp}\simeq
0.6\,\omega_{\rm pp}/c$, as obtained by solving numerically the dispersion
relation of the CFI Bret  _et al._ (2010b). This computation also yields
$\tilde{\chi}_{\rm p}\simeq 0.006$ and $\chi_{\rm b}\simeq 0.07$, thus showing
that the kinetic limit does apply for both components. For reference, the
approximations of Eq. (12) give $\Gamma_{\rm w,k-k}\simeq 0.01\omega_{\rm pp}$
and $k_{\perp}\simeq 0.7\,\omega_{\rm pp}/c$ in that regime. These predictions
fairly match the simulations results: the growth rate evaluated using Eq. (8)
between $t=200\,\omega_{\rm pp}^{-1}$ and $t=450\,\omega_{\rm pp}^{-1}$ is
$\Gamma_{\rm w}^{\rm PIC}\simeq 8\times 10^{-3}\,\omega_{\rm pp}$, while the
dominant $k_{\perp}$ in the Fourier spectrum of $B_{z}$ at saturation
($t\simeq 500\,\omega_{\rm pp}^{-1}$) is measured to be $k_{\perp}^{\rm
PIC}\simeq 0.8\,\omega_{\rm pp}/c$. Considering that the spectrum of the
instability is rather broad and variable with time, the factor of $\sim 2$
discrepancy between the theoretical and simulation results is not very
significant.

The measured value of $\overline{B}_{z}$ is compared to the saturation limits
$B_{\rm t}^{\rm b}$ and $B_{\rm p}^{\rm b}$ in the upper panel of Fig. 1, and
to $B_{\rm t}^{\rm p}$ and $B_{\rm m}^{\rm p}$ in the lower panel. As
explained earlier, we only plot the maximum of the two “particle limit”
criteria corresponding to either component, since the lower one is not
relevant for determining saturation. In the present case, the current density
carried by the beam largely dominates that of the plasma because $|\beta_{\rm
p}|\ll 1$. We do not plot $B_{\rm m}^{\rm b}$ because it lies far above
$B_{\rm t}^{\rm b}$, as expected from the discussion of Sec. II.4. Recalling
that the limits given in Eqs. (15), (16) and (19) are upper limits, saturation
is expected to occur once the measured $B$ value exceeds one of the
corresponding curves in Fig. 1. All limits shown here are computed from the
instantaneous quantities measured in the simulation, which explains their
evolution in time. A word of caution thus appears necessary regarding $B_{\rm
t}$: as it scales with $\Gamma_{\rm w}^{2}$, which is computed through Eq.
(8), this limit becomes meaningless outside the phase of linear growth of the
CFI. In particular, the fact that $\overline{B}_{\rm t}^{\rm
b}<\overline{B}_{z}$ at early times ($t\lesssim 100\,\omega_{\rm pp}^{-1}$),
does not mean that saturation has occurred. On the other hand, the fact that
$\overline{B}_{\rm t}^{\rm b}$ and $\overline{B}$ cross each other at $t\simeq
400\,\omega_{\rm pp}^{-1}$ is indicative of saturation through trapping.

Figure 2: Comparison of the typical filament size, as extracted from
simulation (black line), with the Larmor radius of plasma particles in the
simulation frame (red dashed line) and of beam particles in the instantaneous
Weibel frame (orange dotted line). Both radii are computed using dynamical
quantities extracted from the reference run (a) as defined in Table 1.

Around $t\simeq 400\,\omega_{\rm pp}^{-1}$ the magnetic field indeed becomes
so strong that the quiver frequency of the beam particles exceeds the growth
rate of the instability. Beam particles can then be regarded as transversely
trapped around the $B$-field nodes (Fig. 1 top panel). To quantify this, we
use the characteristic momentum and Lorentz factor averaged over the Maxwell-
Jüttner distribution as $\langle\gamma_{\rm b}\beta_{\rm
b\parallel}\rangle\simeq\langle\gamma_{\rm b}\rangle\simeq 4\gamma_{\rm
b}/\mu_{\rm b}$. Combining those values with the theoretical estimates of
$\Gamma_{\rm w},k_{\perp}$, and the parameters of Table 1, we derive the
trapping limit as $\overline{B}_{\rm t}^{\rm b}\simeq(\Gamma_{\rm
w}/\omega_{\rm pp})^{2}(\omega_{\rm pp}/k_{\perp}c)\langle\gamma_{\rm
b}\rangle/\beta_{\rm b\parallel}\simeq 0.5$, which matches well the observed
saturation value $\overline{B}_{z}\simeq 1$. $\overline{B}_{\rm t}^{\rm b}$ is
also close to the estimate from the measured values of $\Gamma_{\rm w}$,
$k_{\perp}$, $\langle\gamma_{\rm b}\rangle$ and $\beta_{\rm b\parallel}$, that
is, $\overline{B}_{\rm t}^{\rm b,PIC}\simeq 0.3$. As expected, the particle
limit for the beam lies above those values, $\overline{B}_{\rm p}^{\rm
b}\simeq(\pi/2)(n_{\rm b}/n_{\rm p})(\omega_{\rm
pp}/k_{\perp}c)\langle\gamma_{\rm b}\beta_{\rm b\parallel}\rangle\simeq 50$,
and the magnetization limit well above, $\overline{B}_{\rm m,r_{L}}^{\rm
b}\simeq(2/\pi)(k_{\perp}c/\omega_{\rm pp})\langle\gamma_{\rm b}\beta_{\rm
b\parallel}\rangle\simeq 1.6\times 10^{3}$. As anticipated in Sec. II.4, the
trapping limit for the beam thus appears to provide the relevant criterion for
saturation. Interestingly, $\overline{B}_{\rm t}^{\rm p}\ll\overline{B}_{z}$
at all times, even during linear growth, indicating that the strong quiver
motion of the plasma component does not prevent the CFI from growing, neither
does it matter from the point of view of saturation.

The large value of $\overline{B}_{\rm m}^{\rm b}$ confirms that magnetic
trapping does not act longitudinally, meaning that the Larmor radius of the
beam particles remains much larger than the characteristic radius of a
filament; see in particular Fig. 2 which carries out such a comparison. As
already pointed out in Sec. II.3.2, if the drift velocity is relativistic, as
is the case for the beam particles, the magnetization limit is determined by
the spatial constraint $r_{L}\lesssim r$. We recall that the notion of Larmor
radius implies a constant $B$-field along with a null electric field, and
hence has to be computed in the instantaneous Weibel frame, which departs,
given the development of the instability, from the simulation frame. This
change of frame is relevant for the beam, which moves relativistically in the
simulation frame at $u_{b}\simeq 19\simeq\rm const$, while it can be neglected
for the background plasma, given that its velocity and the Weibel frame
velocity remain sub-relativistic in the simulation frame ($|\beta_{\rm
p|w}|\ll 1$).

To better understand why the particle limit does not provide the relevant
saturation criterion here, we quantify the contribution of the beam to the
total current to this effect. We plot in Fig. 3 the particle current density
($n\langle\gamma\beta_{\parallel}\rangle c$) of each species in a limited
region of the periodic $y$-domain. One can see that the contributions of the
beam and the plasma to the electric current density fluctuations are
comparable in scale, although the beam dominates the total particle current
density, which enters Eq. (19). Importantly, charge separation is not complete
and the filaments are rather diluted than spatially split. For this reason,
the magnetic field associated with the maximum particle limit among the
components remains always greater than the simulated value (compare the green
and black curves in Fig. 1), and therefore does not account for saturation.

Figure 3: Transverse profiles of the beam (blue curves) and plasma (yellow
curves) current densities, at saturation $t\simeq 400\,\omega_{\rm pp}^{-1}$
in the reference run (a) listed, and in a limited region of the periodic
$y$-domain. For each species, the solid and dashed curves correspond to
positrons and electrons, respectively.

Concerning the background plasma, it remains sub-relativistic and
relativistically cold in the Weibel frame, hence $\langle\gamma_{\rm
p}\rangle\sim\gamma_{\rm p}$ and $\langle u_{\rm p}\rangle\sim\gamma_{\rm
p}\beta_{\rm p}$. As previously mentioned, the $B$-field associated with
particle trapping inside the filaments is nearly everywhere much smaller than
measured in the simulation: $\overline{B}_{\rm t}^{\rm p}\simeq(\Gamma_{\rm
w}/\omega_{\rm pp})^{2}(\omega_{\rm pp}/k_{\perp}c)\langle\gamma_{\rm
p}\rangle/\beta_{\rm p\parallel}\simeq 0.001$, as can be verified using the
above theoretical estimates for $\Gamma_{\rm max}$ and $k_{\perp}$. As a
matter of fact, the bottom panel of Fig. 1 shows that the background plasma
particles are rapidly trapped inside the filaments, , both transversely and
longitudinally, since $\overline{B}_{\rm m,r_{L}}^{\rm
p}\simeq(2/\pi)(k_{\perp}c/\omega_{\rm pp})\gamma_{\rm p}\beta_{\rm p}\simeq
0.1<\overline{B}_{z}$. Actually, $\overline{B}_{\rm m,r_{L}}^{\rm p}$ rapidly
approaches $\overline{B}_{z}$ (at $\omega_{\rm pp}t\simeq 200)$ and stays
remarkably close to it at later times. We do not interpret this as a cause for
saturation of the CFI, but rather as a relaxation of the low-inertia
background plasma into the strong magnetic fields driven by the large-inertia
beam particles. In runs (b), (c) and (d), $\overline{B}_{\rm m,r_{\rm L}}^{\rm
p}$ gets even smaller than $\overline{B}_{z}$ during linear growth, indicating
that plasma particles become magnetically trapped inside the filaments without
inhibiting the CFI growth.

Well beyond saturation, the characteristic filament radius $r$ increases,
roughly linearly in time (see Fig. 2), as a consequence of filament
coalescence. However, the $B$-field strength as measured in the Weibel frame,
that is, $(B_{z}^{2}-E_{y}^{2})^{1/2}$, remains approximately constant. The
slow evolution of $\overline{B}$ in the simulation frame results from the slow
evolution of the Weibel frame velocity; it is therefore of kinematic origin.
Interestingly, Fig. 2 shows that the typical Larmor radius of background
plasma particles adjusts at all times to the filament radius, $r_{\rm L,p}\sim
r$, which implies that those particles gain energy inside the growing
filaments. Qualitatively, this process can be related to the chaotic dynamics
of particles trapped in an effective potential characterized by the potential
four-vector $A_{x}\sim rB_{z}$, which tends to bring equipartition between
kinetic $\langle p\rangle$ and potential $eA_{x}/c$ energies, under the
approximate conservation of the canonical momentum $\Pi_{x}=p_{x}+eA_{x}/c$.
Such equipartition indeed corresponds to $r_{\rm L,p}\sim r$.

###  III.2 Scan in parameter space

The parameters of the reference run (a) are such that the beam carries most of
the energy density of the system, and its relativistic plasma frequency is the
lower among the two. The smaller inertia of the background plasma particles,
which remain sub- or mildly relativistic in the Weibel frame, explains why
they relax rapidly in the magnetized filamentary structures while the rigid
beam current keeps driving the instability. For this reference run, we thus
find that the transverse trapping of beam particles provides the relevant
criterion for determining the saturation of the CFI. This general picture
proves robust (_i.e._ , it applies from runs (a) to (e) in Table 1) even if
the initial parameters are pushed to extreme values, though always in the CFI-
dominated regime.

For instance, Fig. 4 compares the saturation criteria for run (e), in which
the initial $\gamma_{\rm b}$ and $T_{\rm p}$ have been multiplied by 10 and
$T_{\rm b}$ divided by 10. The instability grows fast, with a measured growth
rate $\Gamma_{\rm w}^{\rm PIC}\simeq 0.3\,\omega_{\rm pp}$, $k_{\perp}^{\rm
PIC}\simeq 0.9\,\omega_{\rm pp}/c$, saturating at $t\simeq 30\,\omega_{\rm
pp}^{-1}$. Here as well, transverse magnetic trapping of beam particles
appears to control the saturation level, while the Larmor radius of background
plasma particles still adapts to the filaments size. The parameters of this
simulation, though, are such that Eq. (12) cannot be applied because the
plasma is hot, and because it moves at relativistic velocities in the Weibel
frame. Solving numerically the dispersion relation of the CFI, we obtain
$\Gamma_{\rm w}\simeq 0.3\,\omega_{\rm pp}$ at $k_{\perp}\simeq
1.2\,\omega_{\rm pp}/c$, which nicely agrees with the PIC values. We then
obtain $\overline{B}_{\rm t}^{\rm b}\simeq 37$, a factor of a few above the
simulated value $\overline{B}_{z}\simeq 10$, and slightly below the
theoretical particle limit $\overline{B}_{\rm p}^{\rm b}\simeq 44$. The time
evolution of these limits, computed with the instantaneous measured values and
plotted in the top panel of Fig. 4, confirms that saturation results from
transverse trapping of the beam particles. Moreover, the closeness of the PIC
field value and plasma magnetization limit (compare $\overline{B}_{z}$ and
$\overline{B}_{\rm m}^{\rm p}$ in the bottom panel of Fig. 4) indicates that
the plasma particles are fully trapped in the filaments, as before.

Figure 4: Same as Fig. 1 for simulation run (e).

Case (f) of Table 1, where $T_{\rm b}$ is reduced by a factor of $30$ while
$\gamma_{\rm b}$ and $T_{\rm p}$ are increased by the same amount, provides an
exception to that general picture. In this particular configuration, both the
beam and the plasma become relativistically hot, leading to comparable initial
relativistic plasma frequencies, namely, $\Omega_{\rm pb}\simeq
0.14\,\omega_{\rm pp}$ and $\Omega_{\rm pp}\simeq 0.2\,\omega_{\rm pp}$. One
can then hardly discern which plays the role of the beam and which plays the
role of the background plasma. What matters for the (transverse or
longitudinal) trapping limits, however, is the inertia of the particles. Here,
$\langle u_{\rm b}\rangle\simeq 40$ and $\langle u_{\rm p}\rangle\simeq 100$
initially, so that the background plasma particles will be trapped later than
the beam particles.

Figure 5: Same as Fig. 1 for simulation run (f). In the top panel is also
plotted the time evolution of the $B$-field associated with the spatial
magnetization limit as applied to the beam ($B_{\rm m}^{\rm b}$, red dashed
curve).

In detail, we measure $\Gamma_{\rm w}^{\rm PIC}\simeq 0.3\,\omega_{\rm pp}$
and $k_{\perp}^{\rm PIC}\simeq 3\omega_{\rm pp}/c$, in fair agreement with the
numerical solution to the CFI dispersion relation ($\Gamma_{\rm w}\simeq
0.4\,\omega_{\rm pp}$ at $k_{\perp}\simeq 1.6\,\omega_{\rm pp}/c$) and which
translates into a plasma trapping limit, $\overline{B}_{\rm t}^{\rm p}\simeq
38$, exceeding the beam trapping limit, $\overline{B}_{\rm t}^{\rm b}\simeq
15$. Moreover, since the plasma now carries a larger current density than the
beam, it gives a greater particle limit: $\overline{B}_{\rm p}^{\rm p}\simeq
180$ vs. $\overline{B}_{\rm p}^{\rm b}\simeq 38$. Those limits have been
evaluated using the simulation parameters; they qualitatively match (yet
overestimate by a factor of a few) the values obtained using the instantaneous
simulation parameters (as plotted in Fig. 5). We therefore expect saturation
to be determined by transverse plasma trapping as confirmed by Fig. 5.

In summary, we observe that the CFI growth rate is set by the species with the
lower (relativistic) plasma frequency, while the saturation level is
determined by that component with the larger inertia per particle, according
to the transverse trapping criterion. The expected overall $B$-field amplitude
at saturation can thus be approximated as

$\frac{B_{\rm sat.}}{\sqrt{4\pi n_{\rm p}mc^{2}}}\simeq\left(\frac{\Gamma_{\rm
w}}{\omega_{\rm pp}}\right)^{2}\frac{\omega_{\rm pp}}{k_{\perp}c}\,,{\rm
max}\left(\langle\gamma_{\rm b}\rangle,\langle\gamma_{\rm p}\rangle\right)\,.$
(22)

In the forthcoming section, we extend this analysis to electron-ion
compositions.

##  IV The electron–ion case

The presence of ions introduces a new scale in the problem, associated with
the hierarchy $m_{i}/m_{e}$ ($m_{i}$ ion mass). If both ions and electrons are
cold, the ratio of ion to electron plasma frequencies scales in proportion to
$\sqrt{m_{e}/m_{i}}$. If the electrons are heated to such a degree that their
effective inertia becomes similar to that of the ion species, then the above
hierarchy disappears: both species share a similar relativistic plasma
frequency, and hence the electron-ion component effectively behaves as a pair
plasma. Thus, one may expect to obtain results similar to those for the pair
systems examined in the previous section.

In the particular context of relativistic shock physics, it is known that
electrons are efficiently heated up to near equipartition in
ultrarelativistic, weakly magnetized conditions (see _e.g._ Sironi  _et al._
(2015); Vanthieghem  _et al._ (2020) and references therein). By contrast, in
the mildly relativistic and magnetized regime, electron heating appears to be
weak, implying that some hierarchy between the response of electrons and ions
remains preserved. Both situations will be addressed in the following. In
order to be able to capture the physics of the instability for both electron
and ion species, with a sufficient number of macro-particles per cell and
spatial extent, we will adopt an ion-to-electron mass ratio $m_{i}/m_{e}=100$.

Run |  $\gamma_{\rm b}$ |  $\gamma_{\rm p}$ |  $\frac{T_{\rm be}}{\gamma_{\rm b}}$ |  $\frac{T_{\rm bi}}{\gamma_{\rm b}}$ |  $\frac{T_{\rm pe}}{\gamma_{\rm p}}$ |  $\frac{T_{\rm pi}}{\gamma_{\rm p}}$ |  $\frac{\gamma_{\rm b}n_{\rm b}}{\gamma_{\rm p}n_{\rm p}}$
---|---|---|---|---|---|---|---
(i1) |  $18.9$ |  $1.011$ |  $2.4$ |  $0.024$ |  $0.20$ |  $0.0020$ |  $1.9$
(i2) |  $6.7$ |  $1.00075$ |  $4.5$ |  $0.15$ |  $0.2$ |  $0.002$ |  $1.3$
Table 2: Parameters of the electron-ion simulations, once transformed to the
Weibel frame. Run (i1) is analogous to run (a) in which the positrons have
been replaced with ions. Run (i2) treats a mildly relativistic regime.
Electron and ion temperatures are given in units of $m_{e}c^{2}/k_{\rm B}$ and
$m_{i}c^{2}/k_{\rm B}$, respectively.

###  IV.1 Ultrarelativistic regime

Let us first examine the saturation criteria for case (i1) described in Table
2. The parameters of this run are obtained from run (a) by replacing the
positrons with ions of charge $+e$ and mass $m_{i}=100\,m_{e}$. The beam
electrons and ions then have a comparable inertia: $\langle p_{\rm
be}\rangle/m_{e}c\simeq 4\gamma_{\rm b}T_{\rm be}\simeq 3430$ and $\langle
p_{\rm bi}\rangle/m_{e}c\simeq\gamma_{\rm b}\beta_{\rm
b}(m_{i}/m_{e})K_{3}(\mu_{\rm bi})/K_{2}(\mu_{\rm bi})\simeq 4520$ ($K_{n}$ is
the modified Bessel function of the $n$th kind). Note that the beam ions have
a proper temperature $T_{\rm bi}/m_{i}c^{2}\simeq 0.45$, so that they cannot
be considered as fully relativistic.

Accounting for the ion mass modifies the trapping and magnetization limits as

$B_{\rm t}^{\rm i}=\frac{\Gamma_{\rm
w}^{2}}{k_{\perp}}\frac{\langle\gamma\rangle m_{i}}{\beta_{\parallel}e}$ (23)

and

$B_{\mathrm{m},\,r_{\mathrm{L}}}^{\rm
i}=\frac{2}{\pi}k_{\perp}\langle\gamma\beta\rangle\frac{m_{i}c^{2}}{e}\,.$
(24)

The time evolution of the simulated mean $B$-field is plotted in Fig. 6.
Unlike previous studies (_e.g._ Ruyer  _et al._ (2015a)), the system does not
experience an early phase governed by electrons, in which the CFI grows
faster, before moving to a regime ruled by the slower ion-driven CFI. We
ascribe this behavior to the similar inertia of the beam ions and electrons.
Solving the CFI dispersion relation in the presence of ions yields a maximum
growth rate $\Gamma_{\rm w}\simeq 0.025\,\omega_{\rm pp}$ for a wavenumber
$k_{\perp}\simeq 0.5\,\omega_{\rm pp}/c$ (as before, $\omega_{\rm pp}$ denotes
the _electron_ plasma frequency of the background plasma). These values are
very close to the simulation values, namely, $\Gamma_{\rm w}^{\rm PIC}\simeq
0.024\,\omega_{\rm pp}$ (as obtained by exponentially fitting $B_{z}(t)$ over
$300<\omega_{\rm pp}t<500$) and $k_{\perp}^{\rm PIC}\simeq 0.4\,\omega_{\rm
pp}/c$ (as measured from the spatial Fourier spectrum of $B_{z}$).

Figure 6: Temporal evolution of the simulated mean $B$-field strength (black
curves) compared to various saturation criteria for run (i1) defined in Table
2. Top panel: particle (green dashed-dotted curve) and trapping (blue dotted
curve) limits as applied to the beam ions. Bottom panel: spatial magnetization
(red dashed curve) and trapping (blue dotted curve) limits as applied to the
plasma ions. All curves are in units of $m_{e}c\omega_{\rm pp}/e$.

As in run (a), the CFI saturates through transverse trapping of the beam
particles (electrons and ions). This is consistent with the fact that the
theoretical trapping limit, $\overline{B}_{\rm t}^{\rm b}\simeq 6$, is much
smaller than the particle limit, $\overline{B}_{\rm p}^{\rm b}\simeq 75$, both
limits being computed for the beam ions and using the initial simulation
parameters. This estimate of $\overline{B}_{\rm t}^{\rm b}$ matches well that
evaluated at saturation time ($t\simeq 400\,\omega_{\rm pp}^{-1}$)
instantaneous simulation parameters (see top panel of Fig. 6). At later times,
again similarly to run (a), the background plasma particles turn fully
magnetized, with their typical Larmor radius adjusting to the mean filament
size (bottom panel of Fig. 6). A notable difference with run (a), however, is
that the mean $B$-field strength here remains quasi-constant following
saturation (up to the final simulation time, $t=2400\,\omega_{\rm pp}^{-1}$),
rather than slowly increasing as in Fig. 1.

In short, in this asymmetric, relativistic electron-ion simulation, in which
both species share a similar inertia, we recover the general picture of the
previous section. Accordingly, the CFI saturation is determined by the
trapping limit as applied to the species with the largest inertia.

###  IV.2 Mildly relativistic regime

Figure 7: Temporal evolution of the simulated mean $B$-field strength
($B_{z}$, black curves) compared to various saturation criteria for run (i2)
defined in Table 2. Top panel: comparison of the spatial magnetization
($B_{\rm m}^{be}$, red dashed curve), particle ($B_{\rm p}^{be}$, green
dashed-dotted curve), trapping ($B_{\rm t}^{\rm be}$, blue dotted curve)
limits as applied to the beam electrons, plus the trapping limit applied to
beam ions ($B_{\rm t}^{\rm bi}$, light-blue dotted curve). Also plotted is the
saturated $B$-field from Eq. (29) ($B_{\rm c}$, magenta dotted line). Bottom
panel: spatial magnetization ($B_{\rm m}^{\rm pe}$, red dashed curve) and
trapping ($B_{\rm t}^{\rm pe}$, blue dotted curve) limits as applied to the
plasma electrons. All curves are in units of $m_{e}c\omega_{\rm pp}/e$.

We now address the case of two electron-ion plasmas counterstreaming with a
moderate Lorentz factor ($\sim 3$) in a reference frame. These two plasma
flows mainly differ in their temperatures: the beam’s electron and ion
populations are much hotter than their plasma counterparts, and for each (beam
or plasma) component, the electrons are also much hotter than the ions. In
particular, the difference in temperature between the beam ions and electrons
is justified by the fact that, according to kinetic simulations, the shock-
reflected ions have a temperature at least three times larger than their
electronic counterpart in the downstream frame (see Caprioli and Spitkovsky
(2014); Crumley  _et al._ (2019); Ligorini  _et al._ (2021a, b) and references
therein). The initial parameters for this run (i2), as expressed in the
corresponding Weibel frame, are summarized in Table 2.

In this configuration, one has $\langle p_{\rm be}\rangle\simeq 800\,m_{e}c$,
while $\langle p_{\rm bi}\rangle\simeq 3000\,m_{e}c$. A hierarchy therefore
persists between the beam electrons and ions, leading to a somewhat different
picture for the evolution of the instability and its saturation level.

Figure 7 shows that after a transient early phase ruled by oblique modes, the
CFI sets in at $t\simeq 300\,\omega_{\rm pp}^{-1}$ and rapidly saturates at
$t\simeq 400\,\omega_{\rm pp}^{-1}$ with a measured growth rate $\Gamma_{\rm
w}^{\rm PIC}\simeq 5\times 10^{-3}\,\omega_{\rm pp}$ and a dominant wave
number $k_{\perp}^{\rm PIC}\simeq 0.35\,\omega_{\rm pp}/c$. During this short
period, the $B$-field grows only by a factor of a few, likely because the
transverse trapping limit for beam electrons is already partially fulfilled,
see top panel of Fig. 7. This figure also suggests that the trapping of beam
ions contributes to the instability saturation, as would be expected from
their larger inertia.

Figure 8: Out-of-plane magnetic field ($B_{z}$) generated by the
counterstreaming of mildly relativistic electron-ion flows (i2). The magnetic
field is plotted at two different times: in the early ($\omega_{\rm pp}t=594$,
top) and late ($\omega_{\rm pp}t=990$, bottom) phases of the CFI when cavities
have started to form.

At later times ($t\gtrsim 700\omega_{\rm pp}^{-1}$), a secondary instability
develops, leading the mean $B$-field strength to rise by almost two orders of
magnitude. As shown in Fig. 8, this instability generates isolated, large-
scale magnetic filamentary structures, which are essentially filled with beam
electrons and plasma ions, and devoid of beam ions and plasma electrons. Those
structures, or “cavities”, have been observed in previous electron-ion
simulations Ruyer  _et al._ (2015b); Naseri  _et al._ (2018) and studied
recently in greater detail in Ref. Peterson  _et al._ (2021). Although the
latter paper considered a simpler setting consisting of an electron beam-
plasma system embedded in an ion background, the picture that it sketches can
be readily extended to the present problem. Specifically, the cavities are
driven by the beam electrons, which are initially overdense relative to the
plasma ($\gamma_{\rm be}n_{\rm be}/\gamma_{\rm pe}n_{\rm pe}\simeq 1.3$ at
$t=0$). As a cavity expands due to the magnetic pressure exerted by the beam
electron current, more beam electrons join and add their contribution to the
current inside the cavity, thus feeding back positively on the magnetic field.
Meanwhile, the beam ions are expelled from the cavity by the growing field,
just as the plasma electrons. The background ions accumulate in the cavity,
mainly (initially) as a result of the confining force exerted by the $E_{y}$
electric field component. This scenario is illustrated in Fig. 9, in the case
of the cavity formed at $(x,y)\simeq(450,-60)\,c/\omega_{\rm pp}$ in the
bottom panel of Fig. 8. The top left panel depicts the time evolution of the
$B$-field profile across the cavity, while the bottom left panel plots the
density profiles of the various populations of the system.

Interestingly, this secondary instability is essentially driven by one
species, here the beam electrons, and it leads to a sharp contrast between the
beam electron density inside and outside the cavity. Thus, it is not
surprising that the particle limit, as evaluated for the beam electrons,
nicely follows the evolution of the magnetic field during this nonlinear
phase111In Fig. 7 there is an offset of about an order of magnitude between
the measured value $B_{z}$ (black solid curve) and the theoretical limit
corresponding to (19) (green dashed-dotted curve). This offset is related to
the overall geometry, in particular the fact that the structures are not
space-filling while the averages are taken over the simulation box. It is
clear, however, that inside a cavity, the magnetic field is mostly carried by
the beam electrons., yet this does not cause the instability to saturate.

We also note that a key factor for this secondary instability is a clear
hierarchy between the beam ions and the beam electrons. Were they of equal
inertia, these two species would react similarly in adjacent filaments,
leading to the growth of all filaments as in the standard CFI. A comparison of
this simulation with the previous one (i1) suggests that in order for the
instability to develop, the beam electrons and ions should differ in their
inertia by at least a factor of a few.

Figure 9: Top left panel: magnetic field profiles along the transverse
direction ($y$) and at successive times, as indicated, for simulation run
(i2). The figure reveals the growth of the magnetic field as the cavity
expands. Bottom left panel: transverse profiles of the number density of the
beam and plasma components at the onset of saturation, $\omega_{\rm pp}t\simeq
1000$. Right panel: $(y,p_{x})$ phase space of the plasma ions at the same
time.

This instability causes the magnetic field to grow rapidly until saturation is
reached at $t\simeq 1100\,\omega_{\rm pp}^{-1}$. While in Peterson  _et al._
(2021), the magnetic pressure pushes a “wall” composed of background ions
initially at rest, in the present case it evacuates the beam ions, which are
relativistic. We can thus adapt the calculation of the instability growth rate
made in that study to our conditions by taking into account the inertia of the
beam ions, as follows.

Assuming that the $B$-field inside the cavity is mainly generated by the beam
electrons, the magnetic pressure acting on this wall can be expressed as

$\frac{B_{z}^{2}}{8\pi}=\frac{(4\pi en_{\rm be}\gamma_{\rm be}\beta_{\rm
be}r_{\rm c})^{2}}{8\pi}\,.$ (25)

The momentum per unit area of the wall is mainly carried by the expelled beam
ions, and so can be estimated as $\gamma_{\rm bi}n_{\rm bi}\langle p_{\rm
bi}\rangle r_{\rm c}(t)$, where $r_{\rm c}(t)$ is the instantaneous cavity
radius. Momentum balance in the transverse ($y$) direction then leads to

$\frac{d}{dt}\left(\gamma_{\rm bi}n_{\rm bi}\langle p_{\rm bi}\rangle r_{\rm
c}\frac{dr_{\rm c}}{dt}\right)=2\pi(en_{\rm be}\gamma_{\rm be}\beta_{\rm
be}r_{\rm c})^{2}\,,$ (26)

The solution to this equation grows as $r_{\rm c}\propto e^{\Gamma_{\rm c}t}$,
where the growth rate is given by

$\Gamma_{\rm c}=\frac{\Omega_{\rm pbi}}{2}\,.$ (27)

In the present case, $\Omega_{\rm pbi}\simeq 0.020\,\omega_{\rm pp}$, which is
in fair agreement with the growth rate $\Gamma^{\rm PIC}\simeq
0.015\,\omega_{\rm pp}$ measured in the simulation over the interval
$750\lesssim\omega_{\rm pp}t\lesssim 960$.

According to Peterson  _et al._ (2021), saturation is reached once the
background plasma ions are accelerated by the inductive electric field
($E_{x}$) to a point where they become relativistic ($p_{i,x}\simeq m_{i}c$)
and neutralize the electron beam current. The right panel of Fig. 9, which
displays the $(y,p_{x})$ phase space of the plasma ions at the onset of
saturation ($\omega_{\rm pp}t=990$), confirms that they have indeed attained
relativistic momenta by that time inside the cavity. Adapting again the
calculations in Ref. Peterson  _et al._ (2021), the radius of the cavity at
saturation can be expressed as

$r_{\rm c,sat}=\frac{c}{\gamma_{\rm b}^{1/2}\omega_{\rm
pbi}}=\frac{c}{\langle\gamma_{\rm bi}\rangle^{1/2}\Omega_{\rm pbi}}\,,$ (28)

recalling that $\omega_{\rm pbi}=\sqrt{4\pi n_{\rm bi}e^{2}/m_{i}}$. This
gives $r_{\rm c,sat}\simeq 10\,c/\omega_{\rm pp}$, which agrees relatively
well with the size of the structures seen in Figs. 8 and 9.

The corresponding saturated value of the magnetic field is given by

$\overline{B}_{\rm c}\simeq\gamma_{\rm
be}^{1/2}\left(\frac{m_{i}}{m_{e}}\right)^{1/2}\frac{\omega_{\rm
pbe}}{\omega_{\rm pp}}$ (29)

in normalized units. One obtains $\overline{B}_{\rm c}\simeq 11$ in correct
agreement with the observed value $\overline{B}_{z}\simeq 6$ (see top panel of
Fig. 7). Note that in Ref. Peterson  _et al._ (2021) an extra factor of $1/2$
was added in the estimation of the saturated field, which is not included
here.

##  V Late–time evolution of the beam-plasma system

We conclude by investigating briefly the late time evolution of the beam-
plasma system after the saturation of the magnetic field. It is worth noting
that in this final stage, both the beam and plasma components are expected to
relax to isotropy in the turbulence frame. This can be seen as a transition
from the two-stream collisionless system to a long-term hydrodynamical system
in which everything has been effectively mixed. In this respect, if we assume
that the beam and plasma have relaxed to the same final velocity but with
different temperatures, the conservation of energy and momentum implies:

$\displaystyle\gamma_{\rm bi}^{2}w_{\rm bi}-p_{\rm bi}+\gamma_{\rm
pi}^{2}w_{\rm pi}-p_{\rm pi}$ $\displaystyle=\gamma_{\rm f}^{2}(w_{\rm
bf}+w_{\rm pf})$ (30) $\displaystyle-p_{\rm bf}-p_{\rm pf}\,,$
$\displaystyle\gamma_{\rm bi}^{2}\beta_{\rm bi}w_{\rm bi}+\gamma_{\rm
pi}^{2}\beta_{\rm pi}w_{\rm pi}$ $\displaystyle=\gamma_{\rm f}^{2}\beta_{\rm
f}(w_{\rm bf}+w_{\rm pf})\,,$ (31)

where the subscripts i and f here refer, respectively, to the initial and
final states of the beam (b) and plasma (p) components. As before, $w$ denotes
the enthalpy density and $p$ the pressure. Note that we have neglected the
contribution of magnetic turbulence in the final state, as it is expected to
be subdominant.

In the case where the final states of the beam and plasma are relativistically
hot, and therefore share the same adiabatic index, $\widehat{\Gamma}_{\rm
f}=w_{\rm f}/(w_{\rm f}-p_{\rm f})$ ($w_{\rm f}$ and $p_{\rm f}$ are the total
final enthalpy density and pressure), the final velocity $\beta_{\rm f}$
satisfies

$\frac{\gamma_{\rm bi}^{2}w_{\rm bi}-p_{\rm bi}+\gamma_{\rm pi}^{2}w_{\rm
pi}-p_{\rm pi}}{\gamma_{\rm bi}^{2}\beta_{\rm bi}w_{\rm bi}+\gamma_{\rm
pi}^{2}\beta_{\rm pi}w_{\rm pi}}=\frac{\kappa_{\rm f}-1+\beta_{\rm
f}^{2}}{\kappa_{\rm f}\beta_{\rm f}}\,,$ (32)

where $\kappa_{\rm f}\equiv\widehat{\Gamma}_{\rm f}/(\widehat{\Gamma}_{\rm
f}-1)$.

Consider for instance the case, exemplified by run (a) of Table 1, of an
initially sub-relativistic ($\beta_{\rm p,i}\sim 0$) and cold ($p_{\rm
p,i}\sim 0$) plasma interacting with a relativistically hot beam which carries
most of the energy (i.e. $\gamma_{\rm b,i}^{2}\,p_{\rm b,i}\gg w_{\rm p,i}$).
We then have $\kappa_{\rm f}\simeq 4$ (as in the initial state), so that

$\displaystyle\beta_{\rm f}$ $\displaystyle\,\simeq\,\beta_{\rm
b}\left(1-\frac{1}{2}\,\frac{w_{\rm pi}}{\gamma_{\rm bi}^{2}p_{\rm
bi}}\right),$ $\displaystyle\gamma_{\rm f}$
$\displaystyle\,\simeq\,\gamma_{\rm b}\left(1-\,\frac{w_{\rm pi}}{p_{\rm
bi}}\right)\,.$ (33)

The second equation further assumes $w_{\rm pi}\ll p_{\rm bi}$. The Lorentz
factors are evaluated in the simulation frame of the two-stream system.

Figure 10: Time evolution of various four-velocities as extracted from run (a)
defined in Table 1. Light-green dotted curve: four-velocity of the Weibel
frame. Yellow curve: four-velocity of the plasma. Blue curve: four-velocity of
the beam. Magenta dashed line: four-velocity of the relaxed plasma and beam as
given by Eq. (33).

The above indicates that the asymptotic velocity of the relaxed components
should be close to the initial beam velocity. This behavior is illustrated in
Fig. 10, which shows the time evolution of various four-velocities as
extracted from our reference run (a). The beam four-velocity $u_{\rm b}$
(blue) indeed approaches from above the predicted asymptotic four-velocity,
$u_{\rm f}=\beta_{\rm f}\gamma_{\rm f}$ (magenta dashed line). Conversely, the
plasma four-velocity $u_{\rm p}$ is seen to increase steadily toward $u_{\rm
f}$. Also overlaid is the instantaneous four-velocity of the Weibel frame
(green dotted line), computed from the simulation data as $u_{\rm
w}=\gamma_{\rm w}\beta_{\rm w}$ with $\beta_{\rm w}=(\langle
E_{y}^{2}\rangle/\langle B_{z}^{2}\rangle)^{1/2}$ (the average is taken over
the simulation domain). Note that this quantity is not defined at early times
because of the dominance of oblique modes characterized by $\langle
E_{y}^{2}\rangle/\langle B_{z}^{2}\rangle>1$. The four-velocity of the Weibel
frame tracks that of the background plasma quite well. We note that the
convergence to the hydrodynamical regime is not attained over the time scale
of the simulation. As a matter of fact, we expect the convergence to proceed
at a slower rate as time increases. This is because relaxation takes place in
the Weibel frame, hence time dilation effects associated with the relativistic
velocity of the Weibel frame relative to the simulation frame will slow down
the apparent relaxation rate.

##  VI Conclusions

In this paper, we have investigated the saturation of the current
filamentation instability, or Weibel instability, in an asymmetric
configuration, meaning in the case in which the counterstreaming plasmas
differ in terms of velocity, temperature and density. This configuration is
notably representative of the precursor region of electron-positron or
electron-ion shocks, although the implications of our results are not
restricted to such systems. Our study relies on large-scale periodic PIC
simulations of counterstreaming flows composed of a hot dilute population
representing the beam (_e.g._ the particles reflected at the shock front) and
a relatively cold plasma _e.g._ the background plasma that is incoming toward
the shock). The parameters of our fiducial run have been directly borrowed
from a large-scale relativistic shock simulation at a position deep in the
precursor; the parameters of subsequent runs have then been varied in an _ad
hoc_ manner to explore different possible settings. We have discussed several
theoretically motivated criteria for saturation and compared them to the
simulation results.

The asymmetric counterstreaming configuration departs from its symmetric
counterpart in two important ways: (1) there exists an ambiguity as to whether
a given criterion should be applied to the beam, or to the plasma component;
(2) there exists a preferred reference frame, dubbed here the “Weibel frame”
Pelletier  _et al._ (2019), in which the instability is purely magnetic; this
reference frame does not _a priori_ coincide with that in which the total
momentum flux vanishes, as happens for the symmetric configuration. Here, we
pay particular attention to that latter point. We have set up our simulations
such that for each set of parameters characterizing the plasma flows, the
simulation frame initially coincides with the Weibel frame.

We have then compared different mechanisms as possible sources of saturation
of the magnetic field associated with the instability: magnetic trapping,
particle limit, Alfvén limit. Our general conclusion is that, for pair
plasmas, the saturation level is determined by the criterion of magnetic
trapping as applied to the (beam or plasma) component that carries the larger
inertia of the two: the growth rate is found to diminish strongly once the
quiver frequency of that component becomes comparable with, or larger than the
instability growth rate. For all studied cases, our theoretical estimates of
the instability properties, such as the maximum growth rate and associated
wave number, are consistent with those extracted from the simulations.
Consequently, it is possible to obtain reasonable analytical approximations
for the strength of the magnetic field at saturation. Furthermore, we find
that the particle limit is never fulfilled, all the more so when the component
of larger inertia is relativistically hot, as its temperature then prevents
its charged species from being fully segregated in separate filaments. We have
observed that the component of smaller inertia becomes rapidly trapped inside
the filaments, in some cases even during the linear phase of the CFI. At late
times, the Larmor radius of those particles closely follows the characteristic
filament radius and thus grows in time through coalescence. Asymptotically,
the system tends to a final state where the two fluids are effectively mixed,
drifting at the same mean velocity. However, due to relativistic time dilation
effects, this ultimate regime could not be accessed from our simulations.

We have also investigated the case of asymmetric electron-ion systems with a
mass ratio $m_{i}/m_{e}=100$. As long as there is not a clear hierarchy in
inertia between the electron and ions species of a given (beam or plasma)
component at the beginning of the simulation, the development of the
instability and the saturation proceed much as in the case of a pair plasma.
The picture and saturation criterion discussed above thus remain applicable.
However, if the electron and ion inertia differ by a factor of a few or more,
a different instability eventually supersedes the CFI. It leads to the
formation of cavities in which the beam electrons and background plasma ions
accumulate and drive magnetic field growth, while the beam ions are pushed
outwards along with the plasma electrons. This mechanism comes to an end when
the plasma ions inside the cavities, accelerated by the inductive electric
field, become capable of neutralizing the electron beam current, as discussed
recently in Ref. Peterson  _et al._ (2021).

###### Acknowledgements.

This work was supported by the ANR (UnRIP project, Grant No.
ANR-20-CE30-0030). We acknowledge GENCI-TGCC for granting us access to the
supercomputer IRENE under Grants No. 2019-A0050407666, 2020-A0080411422 and
2021-A0080411422 to run PIC simulations. We acknowledge financial support from
Centre National d’Études Spatiales (CNES). This work has also been supported
in part by the Sorbonne Université DIWINE Emergence-2019 program.

## References

  * Weibel (1959) E. S. Weibel, Phys. Rev. Lett. 2, 83 (1959).
  * Fried (1959) B. D. Fried, Phys. Fluids 2, 337 (1959).
  * Califano  _et al._ (1997) F. Califano, F. Pegoraro, and S. V. Bulanov, Phys. Rev. E 56, 963 (1997).
  * Bret  _et al._ (2004) A. Bret, M.-C. Firpo, and C. Deutsch, Phys. Rev. E 70, 046401 (2004).
  * Bret  _et al._ (2005) A. Bret, M.-C. Firpo, and C. Deutsch, Phys. Rev. E 72, 016403 (2005).
  * Bret  _et al._ (2007) A. Bret, L. Gremillet, and J. C. Bellido, Phys. Plasmas 14, 032103 (2007).
  * Bret  _et al._ (2008) A. Bret, L. Gremillet, D. Bénisti, and E. Lefebvre, Phys. Rev. Lett. 100, 205008 (2008).
  * Bret  _et al._ (2010a) A. Bret, L. Gremillet, and M. E. Dieckman, Phys. Plasmas 17, 120501 (2010a).
  * Bret  _et al._ (2010b) A. Bret, L. Gremillet, and D. Bénisti, Phys. Rev. E 81, 036402 (2010b).
  * Wiersma and Achterberg (2004) J. Wiersma and A. Achterberg, Astron. Astrophys. 428, 365 (2004).
  * Achterberg and Wiersma (2007) A. Achterberg and J. Wiersma, Astron. Astrophys. 475, 1 (2007).
  * Achterberg  _et al._ (2007) A. Achterberg, J. Wiersma, and C. A. Norman, Astron. Astrophys. 475, 19 (2007).
  * Lee and Lampe (1973) R. Lee and M. Lampe, Phys. Rev. Lett. 31, 1390 (1973).
  * Molvig (1975) K. Molvig, Phys. Rev. Lett. 35, 1504 (1975).
  * Silva  _et al._ (2002) L. O. Silva, R. A. Fonseca, J. W. Tonge, W. B. Mori, and J. M. Dawson, Phys. Plasmas 9, 2458 (2002).
  * Adam  _et al._ (2006) J. C. Adam, A. Héron, and G. Laval, Phys. Rev. Lett. 97, 205006 (2006).
  * Debayle  _et al._ (2010) A. Debayle, J. J. Honrubia, E. d’Humières, and V. T. Tikhonchuk, Phys. Rev. E 82, 036405 (2010).
  * Fiuza  _et al._ (2012) F. Fiuza, R. A. Fonseca, J. Tonge, W. B. Mori, and L. O. Silva, Phys. Rev. Lett. 108, 235004 (2012).
  * Huntington  _et al._ (2015) C. M. Huntington, F. Fiuza, J. S. Ross, A. B. Zylstra, R. P. Drake, D. H. Froula, G. Gregori, N. L. Kugland, C. C. Kuranz, M. C. Levy, C. K. Li, J. Meinecke, T. Morita, R. Petrasso, C. Plechaty, B. A. Remington, D. D. Ryutov, Y. Sakawa, A. Spitkovsky, H. Takabe, and H.-S. Park, Nat. Phys. 11, 173 (2015).
  * Lobet  _et al._ (2015) M. Lobet, C. Ruyer, A. Debayle, E. d’Humières, M. Grech, M. Lemoine, and L. Gremillet, Phys. Rev. Lett. 115, 215003 (2015).
  * Warwick  _et al._ (2017) J. Warwick, T. Dzelzainis, M. E. Dieckmann, W. Schumaker, D. Doria, L. Romagnani, K. Poder, J. M. Cole, A. Alejo, M. Yeung, K. Krushelnick, S. P. D. Mangles, Z. Najmudin, B. Reville, G. M. Samarin, D. D. Symes, A. G. R. Thomas, M. Borghesi, and G. Sarri, Phys. Rev. Lett. 119, 185002 (2017).
  * Li  _et al._ (2019) C. K. Li, V. T. Tikhonchuk, Q. Moreno, H. Sio, E. d’Humières, X. Ribeyre, P. Korneev, S. Atzeni, R. Betti, A. Birkel, E. M. Campbell, R. K. Follett, J. A. Frenje, S. X. Hu, M. Koenig, Y. Sakawa, T. C. Sangster, F. H. Seguin, H. Takabe, S. Zhang, and R. D. Petrasso, Phys. Rev. Lett. 123, 055002 (2019).
  * Fiuza  _et al._ (2020) F. Fiuza, G. F. Swadling, A. Grassi, H. G. Rinderknecht, D. P. Higginson, D. D. Ryutov, C. Bruulsema, R. P. Drake, S. Funk, S. Glenzer, G. Gregori, C. K. Li, B. B. Pollock, B. A. Remington, J. S. Ross, W. Rozmus, Y. Sakawa, A. Spitkovsky, S. Wilks, and H. S. Park, Nat. Phys. 16, 916 (2020).
  * Gruzinov (2001) A. Gruzinov, Astrophys. J. Lett. 563, L15 (2001).
  * Schlickeiser (2005) R. Schlickeiser, Plasma Phys. Control. Fusion 47, A205 (2005).
  * Medvedev and Loeb (1999) M. V. Medvedev and A. Loeb, Astrophys. J. 526, 697 (1999).
  * Gruzinov and Waxman (1999) A. Gruzinov and E. Waxman, Astrophys. J. 511, 852 (1999).
  * Silva  _et al._ (2003) L. O. Silva, R. A. Fonseca, J. W. Tonge, J. M. Dawson, W. B. Mori, and M. V. Medvedev, Astrophys. J. Lett. 596, L121 (2003).
  * Frederiksen  _et al._ (2004) J. T. Frederiksen, C. B. Hededal, T. Haugbølle, and Å. Nordlund, Astrophys. J. Lett. 608, L13 (2004).
  * Kirk and Reville (2010) J. G. Kirk and B. Reville, Astrophys. J. Lett. 710, L16 (2010).
  * Bykov and Treumann (2011) A. M. Bykov and R. A. Treumann, Astron. Astrophys. Rev. 19, 42 (2011).
  * Medvedev  _et al._ (2011) M. V. Medvedev, J. T. Frederiksen, T. Haugbølle, and Å. Nordlund, Astrophys. J. 737, 55 (2011).
  * Plotnikov  _et al._ (2013) I. Plotnikov, G. Pelletier, and M. Lemoine, Mon. Not. Roy. Astron. Soc. 430, 1280 (2013).
  * Lemoine (2013) M. Lemoine, Mon. Not. Roy. Astron. Soc. 428, 845 (2013).
  * Kelner  _et al._ (2013) S. R. Kelner, F. A. Aharonian, and D. Khangulyan, Astrophys. J. 774, 61 (2013).
  * Spitkovsky (2008a) A. Spitkovsky, Astrophys. J. Lett. 682, L5 (2008a).
  * Spitkovsky (2008b) A. Spitkovsky, Astrophys. J. Lett. 673, L39 (2008b).
  * Moiseev and Sagdeev (1963) S. S. Moiseev and R. Z. Sagdeev, J. Nucl. Energy 5, 43 (1963).
  * Sulaiman  _et al._ (2015) A. H. Sulaiman, A. Masters, M. K. Dougherty, D. Burgess, M. Fujimoto, and G. B. Hospodarsky, Phys. Rev. Lett. 115, 125001 (2015).
  * Kato and Takabe (2008) T. N. Kato and H. Takabe, Astrophys. J. Lett. 681, L93 (2008).
  * Kato and Takabe (2010) T. N. Kato and H. Takabe, Astrophys. J. 721, 828 (2010).
  * Bret  _et al._ (2013) A. Bret, A. Stockem, D. Fiuza, C. Ruyer, L. Gremillet, R. Narayan, and L. O. Silva, Phys. Plasmas 20, 042102 (2013).
  * Matsumoto  _et al._ (2015) Y. Matsumoto, T. Amano, T. N. Kato, and M. Hoshino, Science 347, 974 (2015).
  * Bohdan  _et al._ (2020a) A. Bohdan, M. Pohl, J. Niemiec, S. Vafin, Y. Matsumoto, T. Amano, and M. Hoshino, Astrophys. J. 893, 6 (2020a).
  * Bohdan  _et al._ (2020b) A. Bohdan, M. Pohl, J. Niemiec, P. J. Morris, Y. Matsumoto, T. Amano, and M. Hoshino, Astrophys. J. 904, 12 (2020b).
  * Ruyer  _et al._ (2015a) C. Ruyer, L. Gremillet, A. Debayle, and G. Bonnaud, Phys. Plasmas 22, 032102 (2015a).
  * Ruyer  _et al._ (2015b) C. Ruyer, L. Gremillet, and G. Bonnaud, Phys. Plasmas 22, 082107 (2015b).
  * Ruyer  _et al._ (2016) C. Ruyer, L. Gremillet, G. Bonnaud, and C. Riconda, Phys. Rev. Lett. 117, 065001 (2016).
  * Ruyer  _et al._ (2017) C. Ruyer, L. Gremillet, G. Bonnaud, and C. Riconda, Phys. Plasmas 24, 041409 (2017).
  * Lyubarsky and Eichler (2006) Y. Lyubarsky and D. Eichler, Astrophy. J. 647, 1250 (2006).
  * Lemoine  _et al._ (2006) M. Lemoine, G. Pelletier, and B. Revenu, Astrophys. J. 645, L129 (2006).
  * Kato (2007) T. N. Kato, Astrophys. J. 668, 974 (2007).
  * Martins  _et al._ (2009) S. F. Martins, R. A. Fonseca, L. O. Silva, and W. B. Mori, Astrophys. J. Lett. 695, L189 (2009).
  * Nishikawa  _et al._ (2009) K.-I. Nishikawa, J. Niemiec, P. E. Hardee, M. Medvedev, H. Sol, Y. Mizuno, B. Zhang, M. Pohl, M. Oka, and D. H. Hartmann, Astrophys. J. Lett. 698, L10 (2009).
  * Keshet  _et al._ (2009) U. Keshet, B. Katz, A. Spitkovsky, and E. Waxman, Astrophys. J. 693, L127 (2009).
  * Sironi  _et al._ (2013) L. Sironi, A. Spitkovsky, and J. Arons, Astrophys. J. 771, 54 (2013).
  * Plotnikov  _et al._ (2018) I. Plotnikov, A. Grassi, and M. Grech, Mon. Not. Roy. Astron. Soc. 477, 5238 (2018).
  * Tomita  _et al._ (2019) S. Tomita, Y. Ohira, and R. Yamazaki, Astrophys. J. 886, 54 (2019).
  * Lemoine  _et al._ (2019a) M. Lemoine, L. Gremillet, G. Pelletier, and A. Vanthieghem, Phys. Rev. Lett. 123, 035101 (2019a), arXiv:1907.07595 .
  * Lemoine  _et al._ (2019b) M. Lemoine, A. Vanthieghem, G. Pelletier, and L. Gremillet, Phys. Rev. E 100, 033209 (2019b).
  * Lemoine  _et al._ (2019c) M. Lemoine, G. Pelletier, A. Vanthieghem, and L. Gremillet, Phys. Rev. E 100, 033210 (2019c).
  * Pelletier  _et al._ (2019) G. Pelletier, L. Gremillet, A. Vanthieghem, and M. Lemoine, Phys. Rev. E 100, 013205 (2019), arXiv:1907.07750 .
  * Vanthieghem  _et al._ (2020) A. Vanthieghem, M. Lemoine, I. Plotnikov, A. Grassi, M. Grech, L. Gremillet, and G. Pelletier, Galaxies 8, 33 (2020), arXiv:2002.01141 [astro-ph.HE] .
  * Davidson  _et al._ (1972) R. C. Davidson, D. A. Hammer, I. Haber, and C. E. Wagner, Phys. Fluids 15, 317 (1972).
  * Kato (2005) T. N. Kato, Phys. Plasmas 12, 080705 (2005).
  * Califano  _et al._ (1998) F. Califano, F. Pegoraro, S. V. Bulanov, and A. Mangeney, Phys. Rev. E 57, 7048 (1998).
  * Okada and Ogawa (2007) T. Okada and K. Ogawa, Phys. Plasmas 14, 072702 (2007).
  * Shvets  _et al._ (2009) G. Shvets, O. Polomarov, V. Khudik, C. Siemon, and I. Kaganovich, Phys. Plasmas 16, 056303 (2009).
  * Cagas  _et al._ (2017) P. Cagas, A. Hakim, W. Scales, and B. Srinivasan, Phys. Plasmas 24, 112116 (2017).
  * Grassi  _et al._ (2017) A. Grassi, M. Grech, F. Amiranoff, F. Pegoraro, A. Macchi, and C. Riconda, Phys. Rev. E 95, 023203 (2017).
  * Takamoto  _et al._ (2018) M. Takamoto, Y. Matsumoto, and T. N. Kato, Astrophys. J. Lett. 860, L1 (2018).
  * Takamoto  _et al._ (2019) M. Takamoto, Y. Matsumoto, and T. N. Kato, Astrophys. J. 877, 137 (2019).
  * Jaroschek  _et al._ (2004) C. H. Jaroschek, H. Lesch, and R. A. Treumann, Astrophys. J. 616, 1065 (2004).
  * Medvedev  _et al._ (2004) M. V. Medvedev, M. Fiore, R. A. Fonseca, L. O. Silva, and W. B. Mori, Astrophys. J. 618, L75 (2004).
  * Milosavljevic and Nakar (2006) M. Milosavljevic and E. Nakar, Astrophys. J. 641, 978 (2006).
  * Chang  _et al._ (2008) P. Chang, A. Spitkovsky, and J. Arons, Astrophys. J. 674, 378 (2008).
  * Lemoine (2015) M. Lemoine, J. Plasma Phys. 81, 455810101 (2015).
  * Vanthieghem  _et al._ (2018) A. Vanthieghem, M. Lemoine, and L. Gremillet, Phys. Plasmas 25, 072115 (2018).
  * Yang  _et al._ (1994) T. Y. B. Yang, J. Arons, and A. B. Langdon, Phys. Plasmas 1, 3059 (1994).
  * Istomin (2011) Y. N. Istomin, Plasma Phys. Control. Fusion 53, 015006 (2011).
  * Bret  _et al._ (2010c) A. Bret, L. Gremillet, and D. Bénisti, Phys. Rev. E 81, 036402 (2010c).
  * Kaang  _et al._ (2009) H. H. Kaang, C.-M. Ryu, and P. H. Yoon, Phys. Plasmas 16, 082103 (2009).
  * Lefebvre  _et al._ (2003) E. Lefebvre, N. Cochet, S. Fritzler, V. Malka, M.-M. Aléonard, J.-F. Chemin, S. Darbon, L. Disdier, J. Faure, A. Fedotoff, O. Landoas, G. Malka, V. Méot, P. Morel, M. Rabec LeGloahec, A. Rouyer, C. Rubbelynck, V. Tikhonchuk, R. Wrobel, P. Audebert, and C. Rousseaux, Nucl. Fusion 43, 629 (2003).
  * Sironi  _et al._ (2015) L. Sironi, U. Keshet, and M. Lemoine, Space Sci. Rev. 191, 519 (2015).
  * Caprioli and Spitkovsky (2014) D. Caprioli and A. Spitkovsky, Astrophys. J. 783, 91 (2014).
  * Crumley  _et al._ (2019) P. Crumley, D. Caprioli, S. Markoff, and A. Spitkovsky, Mon. Not. Roy. Astron. Soc. 485, 5105 (2019).
  * Ligorini  _et al._ (2021a) A. Ligorini, J. Niemiec, O. Kobzar, M. Iwamoto, A. Bohdan, M. Pohl, Y. Matsumoto, T. Amano, S. Matsukiyo, Y. Esaki, and M. Hoshino, Mon. Not. Roy. Astron. Soc. 501, 4837 (2021a).
  * Ligorini  _et al._ (2021b) A. Ligorini, J. Niemiec, O. Kobzar, M. Iwamoto, A. Bohdan, M. Pohl, Y. Matsumoto, T. Amano, S. Matsukiyo, and M. Hoshino, Mon. Not. Roy. Astron. Soc. 502, 5065 (2021b).
  * Naseri  _et al._ (2018) N. Naseri, S. G. Bochkarev, P. Ruan, V. Y. Bychenkov, V. Khudik, and G. Shvets, Phys. Plasmas 25, 012118 (2018).
  * Peterson  _et al._ (2021) J. R. Peterson, S. Glenzer, and F. Fiuza, Phys. Rev. Lett. 126, 215101 (2021).