1 Introduction
When a shock propagates into a neutral fluid, upstream particles slow down at the shock front as a result of collisions with particles in the slower-moving downstream gas. In fact, binary collisions are the only possible microscopic mechanism for an upstream particle to slow down. As a consequence, the shock front is a few mean-free-paths thick (Zel’dovich & Raizer Reference Zel’dovich and Raizer2002).
In situ measurements of the Earth’s bow shock within the solar wind show that its front is far smaller than the mean free path of the ions at the same location, which is comparable to an astronomical unit (Bale, Mozer & Horbury Reference Bale, Mozer and Horbury2003; Schwartz et al. Reference Schwartz, Henley, Mitchell and Krasnoselskikh2011). Such shocks, where the mean free path is much larger than the front, have been dubbed ‘collisionless shocks’. Instead of being sustained by binary collisions, these shocks are mediated by collective plasma effects acting on much shorter time and length scales than binary Coulomb collisions (Petschek Reference Petschek1958; Sagdeev Reference Sagdeev1966).
Collisionless shocks are believed to occur in a wide variety of astrophysical settings: active galactic nuclei, pulsar wind nebulae, planetary environments, supernova remnants, etc. The absence of collisions allows particles to gain energy without sharing it immediately with other particles. As a result, such shocks have been found to be excellent particle accelerators and now count among the main candidates for the production of high energy cosmic rays (Blandford & Eichler Reference Blandford and Eichler1987; Sironi, Keshet & Lemoine Reference Sironi, Keshet and Lemoine2015; Marcowith et al. Reference Marcowith, Bret, Bykov, Dieckman, Drury, Lembège, Lemoine, Morlino, Murphy and Pelletier2016). They are also believed to play a role in the generation of gamma-ray bursts (Mészáros & Rees Reference Mészáros and Rees2014; Pe’er Reference Pe’er2015) and fast radio bursts (Falcke & Rezzolla Reference Falcke and Rezzolla2014; Lyubarsky Reference Lyubarsky2014).
Starting with the pioneering work of Sagdeev in the 1960s (Sagdeev Reference Sagdeev1966), our knowledge of collisionless shocks has grown tremendously, particularly in the past decade thanks to the advent of large scale particle-in-cell (PIC) simulations (Spitkovsky Reference Spitkovsky, Bulik, Rudak and Madejski2005; Martins et al. Reference Martins, Fonseca, Silva and Mori2009). However, as recently as the 1990s, there were still doubts about the very existence of collisionless shocks (Sagdeev & Kennel Reference Sagdeev and Kennel1991). While the Earth bow shock measurements have definitely eliminated these doubts, the micro-physics of collisionless shock formation, and the mechanism of particle acceleration, are still under investigation.
Given the omnipresence of collisionless shocks and their important role in many phenomena, especially in astrophysics, the conditions for such shocks to form are worthy of investigation. A detailed understanding is all the more important that electrostatic collisionless shocksFootnote 1 have been observed in the laboratory (Ahmed et al. Reference Ahmed, Dieckmann, Romagnani, Doria, Sarri, Cerchez, Ianni, Kourakis, Giesecke and Notley2013), while the production of Weibel mediated shocks such as the ones discussed here, is expected within the next few years (Huntington et al. Reference Huntington, Fiúza, Ross, Zylstra, Drake, Froula, Gregori, Kugland, Kuranz and Levy2015; Lobet et al. Reference Lobet, Ruyer, Debayle, d’Humières, Grech, Lemoine and Gremillet2015; Park et al. Reference Park, Ross, Huntington, Fiuza, Ryutov, Casey, Drake, Fiksel, Froula and Gregori2016). Note that the ‘Weibel instability’ we refer to is sometimes labelled the ‘filamentation instability’ or the ‘beam Weibel’ instability (Silva et al. Reference Silva, Fonseca, Tonge, Mori and Dawson2002; Deutsch et al. Reference Deutsch, Bret, Firpo and Fromy2005; Hill et al. Reference Hill, Key, Hatchett and Freeman2005). It is the instability of two counter-streaming flows with respect to perturbations with wave vectors normal to the flow.
When a collisionless shock forms from the encounter of two plasma shells, the downstream plasma may be thermalized by collisionless processes (see Bret (Reference Bret2015a ) and references therein). As a consequence, the equations of magnetohydrodynamics (MHD) can be applied, so that both collisionless shocks and MHD shocks can in principle be analysed using the same fluid approachFootnote 2 . For the case of a flow-aligned field, MHD prescribes that the fluid and the field are decoupled (Majorana & Anile Reference Majorana and Anile1987), so that the very same shock should form, regardless of the field intensity.
Here, we present a specific example of departure from this expected MHD behaviour. We consider the encounter of two collisionless cold pair plasmas. A flow-aligned magnetic field is present, and the system is relativistic. In § 3, we explain the predictions of MHD for this system. Then, in § 4, we describe a series of simulations using the PIC technique. These simulations work at the microscopic level, and show a departure from the MHD predictions beyond a critical magnetization. In § 5, we present a micro-physics analysis of the shock formation process explaining the departure from MHD.
Figure 1. Set-up of the system considered. Two collisionless cold pair plasmas of density
$n_{0}$
and initial Lorentz factor
$\unicode[STIX]{x1D6FE}_{0}$
collide over a flow-aligned magnetic field
$\boldsymbol{B}_{0}$
.
2 System considered
The system considered is shown schematically in figure 1. Two identical pair plasma shells of density
$n_{0}$
head toward each other with initial velocity
$\pm \boldsymbol{v}_{0}$
and Lorentz factor
$\unicode[STIX]{x1D6FE}_{0}=(1-v_{0}^{2}/c^{2})^{-1/2}$
. The whole system is embedded in an external field
$\boldsymbol{B}_{0}\Vert \boldsymbol{v}_{0}$
and aligned with the
$x$
axis. We denote by ‘upstream frame’ the frame of reference of the right shell, and by ‘downstream frame’ the frame where the total momentum is zero. When a shock forms, these frames become the upstream and downstream frames of the shock, respectively. The strength of the magnetic field is measured by the magnetization parameter,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=\frac{B_{0}^{2}/4\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{0}n_{0}mc^{2}},\end{eqnarray}$$
where all quantities are measured in the downstream frame.
3 MHD predictions
An MHD plasma sustains 3 kinds of modes: the slow mode, Alfvén mode and fast mode (Kulsrud Reference Kulsrud2005). The phase velocities of these modes satisfy the hierarchy
$v_{slow}<v_{Alfven}<v_{fast}$
. Because of this hierarchy, the Alfvén mode is sometimes dubbed the ‘intermediate mode’ (Kulsrud Reference Kulsrud2005). In the cold limit considered here,
$v_{slow}\rightarrow 0$
and
$v_{fast}\rightarrow v_{Alfven}$
.
A ‘fast shock’ has its front moving faster than the upstream fast mode, while a ‘slow shock’ only moves faster than the upstream slow mode. For fast shocks, the shock front also propagates faster than the downstream Alfvén speed; in slow shocks, it propagates slower. An intermediate regime exist, where the flow is super-Alfvénic upstream and sub-Alfvénic downstream (crossing of the Alfvénic point, see e.g. Kirk & Duffy (Reference Kirk and Duffy1999)). However, such solutions of the MHD jump equations do not survive when produced and are called ‘extraneous’ (Kulsrud Reference Kulsrud2005); they typically split into a pair of ‘fast’ and ‘slow’ shocks.
For a flow-aligned field, the fluid motion decouples from the field (Majorana & Anile Reference Majorana and Anile1987). The shock formed is therefore the same, regardless of the magnetization parameter
$\unicode[STIX]{x1D70E}$
. Nevertheless, its front velocity can still be compared to the phase speeds of the three modes. In the present cold limit, and for
$\unicode[STIX]{x1D6FE}_{0}\rightarrow \infty$
, the shock is expected to be ‘fast’ for
$\unicode[STIX]{x1D70E}<2/3$
, ‘slow’ for
$\unicode[STIX]{x1D70E}>2$
, and ‘extraneous’ in between (see appendix A). Figure 2 shows these limits for a range of
$\unicode[STIX]{x1D6FE}_{0}$
in the
$(\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6FE}_{0})$
plane.
The MHD predictions for the present system are therefore very clear: the same shock should form regardless of the
$\unicode[STIX]{x1D70E}$
parameter, simply because the fluid and the field are perfectly decoupled here. The MHD simulations run in appendix A confirm this conclusion.
Figure 2. The MHD thresholds for slow and fast shocks are represented by the thick black lines, with the thin vertical lines showing the large
$\unicode[STIX]{x1D6FE}_{0}$
limits,
$\unicode[STIX]{x1D70E}=2/3$
, 2. Extraneous shocks occur in between the two thick black lines. The Weibel instability governs systems located above and to the left of the orange curve. The Weibel filaments at saturation are able to stop the incoming flow, and initiate shock formation, only for systems to the left of the blue curve ((5.5) with
$\unicode[STIX]{x1D705}=2/3$
).
4 PIC simulations
We now turn to PIC simulations to conduct a micro-physical, i.e. kinetic, analysis of the system under scrutiny. We use the three-dimensional (3-D) electromagnetic PIC code TRISTAN-MP (Spitkovsky Reference Spitkovsky, Bulik, Rudak and Madejski2005), which is a parallel version of the publicly available code TRISTAN (Buneman Reference Buneman, Matsumoto and Omura1993) that has been optimized for studying relativistic collisionless shocks (Spitkovsky Reference Spitkovsky2008a ,Reference Spitkovsky b ; Sironi & Spitkovsky Reference Sironi and Spitkovsky2009, Reference Sironi and Spitkovsky2011; Sironi et al. Reference Sironi, Spitkovsky and Arons2013). We employ simulations in 2-D computational domains, but all three components of particle velocities and electromagnetic fields are tracked (see more details in appendix B).
We probe the regime
$\unicode[STIX]{x1D6FE}_{0}=10$
and
$0<\unicode[STIX]{x1D70E}<3$
. Note that the parameter space in Bret (Reference Bret2016a
) is parameterized is terms of
$\unicode[STIX]{x1D6FE}_{0}$
and
$\unicode[STIX]{x1D6FA}_{B}$
, the latter being related to the present
$\unicode[STIX]{x1D70E}$
through
$\unicode[STIX]{x1D6FA}_{B}=\sqrt{2\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D70E}}$
Footnote
3
. For clarity, table 1 gives the values of the values of
$\unicode[STIX]{x1D6FA}_{B}$
of Bret (Reference Bret2016a
) corresponding to the values of
$\unicode[STIX]{x1D70E}$
sampled here.
Figure 3. Shock structure from a series of 2-D PIC simulations with
$\unicode[STIX]{x1D6FE}_{0}=10$
,
$0.2\leqslant \unicode[STIX]{x1D70E}\leqslant 3$
, at
$\unicode[STIX]{x1D714}_{p}t=450$
. We plot the
$y$
-averaged density profile (a) and a measure of the plasma anisotropy (b), as defined in (4.1). The vertical dashed line indicates the position of the front, assuming that it propagates at
$c/3$
. The angle between the field and the flow is
$\unicode[STIX]{x1D703}=0$
.
Table 1. Values of the parameter
$\unicode[STIX]{x1D6FA}_{B}=\sqrt{2\unicode[STIX]{x1D6FE}_{0}\unicode[STIX]{x1D70E}}$
used in Bret (Reference Bret2016a
) corresponding to the values of
$\unicode[STIX]{x1D70E}$
sampled here and for
$\unicode[STIX]{x1D6FE}_{0}=10$
.
Figure 3 shows the
$y$
-integrated density profile of the system for
$\unicode[STIX]{x1D6FE}_{0}=10$
(figure 3
a), at a relatively early time,
$\unicode[STIX]{x1D714}_{p}t=450$
, where
$\unicode[STIX]{x1D714}_{p}^{2}=4\unicode[STIX]{x03C0}n_{0}q^{2}/\unicode[STIX]{x1D6FE}_{0}m$
. The magnetization parameter varies from 0.2 to 3, as indicated in the legend. In figure 3(b), we quantify the isotropization of the particle distribution function by plotting the ratio
$\unicode[STIX]{x1D711}$
between the momentum dispersion along the transverse directions (
$y$
and
$z$
) as compared to the longitudinal direction
$x$
, namely,
$$\begin{eqnarray}\unicode[STIX]{x1D711}=\frac{\text{Var}(p_{y})+\text{Var}(p_{z})}{2\text{Var}(p_{x})}.\end{eqnarray}$$
We notice that, for
$\unicode[STIX]{x1D70E}\lesssim 0.4$
, the shock structure is independent of the magnetization, in line with the MHD predictions. In the downstream (left of the vertical dashed line), the density approaches the value predicted by the MHD jump conditions (
${\sim}4.2$
for
$\unicode[STIX]{x1D6FE}_{0}=10$
, and
${\sim}4$
in the limit
$\unicode[STIX]{x1D6FE}_{0}\gg 1$
). Correspondingly, the shock speed approaches the value
${\sim}c/3$
predicted by MHD (indicated by the dashed line). Figure 3(b) shows that for low magnetizations the downstream plasma is nearly isotropic. However, for higher magnetizations (
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
), the downstream density is lower than the value predicted by MHD. Consequently, the shock speed is faster than the MHD prediction
${\sim}c/3$
.
It is noteworthy that the width of the density jump increases with
$\unicode[STIX]{x1D70E}$
. For small values, the shock front is
${\sim}70c/\unicode[STIX]{x1D714}_{p}$
thick. But for
$\unicode[STIX]{x1D70E}=3$
, the transition region between the ‘upstream’ and the ‘downstream’ is
${\sim}300c/\unicode[STIX]{x1D714}_{p}$
.
Why do the results for
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
deviate from MHD? One might think then, that because the PIC simulations are limited to early times, the shock has not formed yet. How much time should the formation of a shock take? For the present system, the growth rate
$\unicode[STIX]{x1D6FF}_{W}$
of the Weibel instability is given by (Stockem, Lerche & Schlickeiser Reference Stockem, Lerche and Schlickeiser2006; Bret Reference Bret2016a
),
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{W}=\unicode[STIX]{x1D714}_{p}\sqrt{2\unicode[STIX]{x1D6FD}_{0}^{2}-\unicode[STIX]{x1D70E}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FD}_{0}=v_{0}/c$
. The shock formation time typically amounts to a few tens of
$e$
-folding times (Bret et al.
Reference Bret, Stockem, Fiúza, Ruyer, Gremillet, Narayan and Silva2013, Reference Bret, Stockem, Narayan and Silva2014). With the parameters used here,
$20\unicode[STIX]{x1D6FF}_{W}^{-1}$
is at most
$28\unicode[STIX]{x1D714}_{p}^{-1}$
for
$\unicode[STIX]{x1D70E}=1.5$
(
$\unicode[STIX]{x1D6FF}_{W}$
vanishes for
$\unicode[STIX]{x1D70E}>2\unicode[STIX]{x1D6FD}_{0}^{2}$
). Therefore, the time
$t=450\unicode[STIX]{x1D714}_{p}^{-1}$
to which the simulations in figure 3 have been run, exceeds by a factor of 15 the slowest expected shock formation time.
Figure 4. Same as in figure 3, but at a later time:
$\unicode[STIX]{x1D714}_{p}t=3600$
.
To verify the above argument, we have evolved the simulations to much longer times:
$\unicode[STIX]{x1D714}_{p}t=3600$
(figure 4). We again find that the density profile strongly varies with
$\unicode[STIX]{x1D70E}$
, contrary to the MHD prescriptions. For magnetizations
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
, the system settles into a quasi-stationary state which does not satisfy the usual MHD jump conditions. Ultimately, the fact that the density jump and the shock speed do not agree with the MHD jump conditions is related to the lack of isotropy in the downstream plasma. As shown in figures 3(b), 4(b), for
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
, the downstream particle distribution is hotter along the longitudinal direction than in the transverse directionsFootnote
4
. For large
$\unicode[STIX]{x1D70E}$
, we find downstream
$\unicode[STIX]{x1D711}<1$
(left of the vertical dashed lines), which means that the flow is not isotropized even at late times.
The width of the density jump is again worth emphasizing. For small values of
$\unicode[STIX]{x1D70E}$
, the shock front on figure 4 is still
${\sim}70c/\unicode[STIX]{x1D714}_{p}$
thick. But for
$\unicode[STIX]{x1D70E}=3$
, the transition region is now
${\sim}2000c/\unicode[STIX]{x1D714}_{p}$
.
The micro-physical analysis discussed next in § 5 predicts that the departure from the MHD behaviour we just observed, is
$\unicode[STIX]{x1D6FE}_{0}$
independent at large
$\unicode[STIX]{x1D6FE}_{0}$
. This prediction has been successfully tested in appendix B by running a series of PIC simulations with
$\unicode[STIX]{x1D6FE}_{0}=30$
. We also confirm in appendix B that these results are not restricted to a perfectly flow-aligned field (
$\unicode[STIX]{x1D703}=0$
) but survive even for a misaligned field.
5 Micro-physics of the shock formation
From the discussion so far, it appears that the observed departure of the system under consideration from the predictions of MHD, boils down to the non-isotropization of the downstream particle distribution function, even at late times. The following kinetic analysis of the shock formation process allows us to understand why isotropization fails.
Weibel shocks are mediated by purely collective phenomena. When the two plasma shells start interpenetrating, the overlapping region turns unstable to counter-streaming instabilities. Many linear instabilities compete (Bret et al.
Reference Bret, Gremillet and Dieckmann2010), but the Weibel (filamentation) instability, with a
$\boldsymbol{k}$
normal to the flow, grows faster than all others, provided (Bret Reference Bret2016a
),
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{0}>\sqrt{\frac{2}{4/3-\unicode[STIX]{x1D70E}}}.\end{eqnarray}$$
The line corresponding to this limit is shown in figure 2 by the orange curve. The Weibel instability dominates the unstable spectrum of the system for all points of the
$(\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6FE}_{0})$
plane above and to the left this line.
The micro-physics of shock formation depends on the ability of the Weibel instability to form magnetic filaments capable of blocking the plasma that keeps entering the overlapping region. In the case of un-magnetized pair plasmas, for example, this condition is already met at saturation of the Weibel instability (Bret et al. Reference Bret, Stockem, Fiúza, Ruyer, Gremillet, Narayan and Silva2013, Reference Bret, Stockem, Narayan and Silva2014). As a result, the density quickly builds up in the overlapping region, and a shock forms. Distribution functions are quickly isotropized in the overlapping region and MHD considerations apply.
The magnetic filaments generated by the Weibel instability are of the form,
$$\begin{eqnarray}\boldsymbol{B}_{f}=B_{f}\sin (ky)\boldsymbol{e}_{z},\end{eqnarray}$$
where
$k$
is the fastest growing wavenumber. When there is no external magnetic field, an analysis of the motion of a particle of mass
$m$
and charge
$q$
in such filaments (Bret Reference Bret2015b
) shows that it is stopped inside if
$$\begin{eqnarray}k^{-1}>\frac{v_{0}}{\unicode[STIX]{x1D714}_{B_{f}}},\quad \text{with }\unicode[STIX]{x1D714}_{B_{f}}=\frac{qB_{f}}{\unicode[STIX]{x1D6FE}_{0}mc},\end{eqnarray}$$
where
$\boldsymbol{v}_{0}$
is the initial, flow-aligned velocity of the particle and
$\unicode[STIX]{x1D6FE}_{0}$
is its Lorentz factor. Although the model from which this conclusion is derived is highly simplified, the condition is consistent with the results of PIC simulations (Bret et al.
Reference Bret, Stockem, Narayan and Silva2014).
How is (5.3) modified in the presence of a flow-aligned magnetic field? One would expect a guiding field to suppress the transverse scattering of particles and to thereby help particles go through the filaments without stopping. Indeed, analysis shows that regardless of their initial velocity or initial position along the
$y$
axis, all particles stream through the filaments whenever (Bret Reference Bret2016b
)
$$\begin{eqnarray}B_{0}>{\textstyle \frac{1}{2}}B_{f}.\end{eqnarray}$$
Since
$B_{f}$
arises from the growth of the Weibel instability, its magnitude can be quantified (Stockem et al.
Reference Stockem, Lerche and Schlickeiser2006; Bret Reference Bret2016a
). Therefore, the above criterion can eventually be expressed in terms of
$\unicode[STIX]{x1D70E}$
and
$\unicode[STIX]{x1D6FD}_{0}^{2}=1-1/\unicode[STIX]{x1D6FE}_{0}^{2}$
, giving (see details in appendix C),
$$\begin{eqnarray}\unicode[STIX]{x1D70E}>\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FD}_{0}^{2},\end{eqnarray}$$
where
$\unicode[STIX]{x1D705}=2/3$
if equipartition is assumed at saturation of the Weibel instability. The boundary corresponding to the criterion (5.5) with
$\unicode[STIX]{x1D705}=2/3$
is shown in figure 2 by the blue curve.
The region between the bounds corresponding to (5.1) and (5.5), i.e. the region between the orange and blue lines in figure 2, corresponds to a range of parameters where the Weibel instability governs the linear phase of the overlapping region, but the filaments at saturation are not strong enough to stop the flow. The expected consequence, as indeed observed in our PIC simulations, is that the flows are not trapped in the overlapping region but keep streaming through. Isotropization is not achieved and MHD does not apply.
The reader may have noticed that for
$\unicode[STIX]{x1D70E}=2$
and 3, the Weibel instability does not govern the linear phase of the initial interaction between to two shells. How is it then that the system still fails to follow MHD? We conjecture that the analysis described above, where we were able to quantify all the steps because the Weibel instability is well understood, must be a particular case of the following more general argument. Instead of the Weibel filaments described by (5.2), consider a turbulent electromagnetic perturbation
$\sum _{\boldsymbol{k}}\boldsymbol{E}_{\boldsymbol{k}}+\boldsymbol{B}_{\boldsymbol{k}}$
(with
$\langle \boldsymbol{E}_{\boldsymbol{k}}\rangle =\langle \boldsymbol{B}_{\boldsymbol{k}}\rangle =\mathbf{0}$
) that is present in the overlapping region and that can potentially isotropize the incoming flow. Consider also a superimposed, flow-aligned field
$\boldsymbol{B}_{0}$
. In the limit
$B_{0}=0$
, the incoming flow is isotropized, and usual MHD applies. In the opposite limit
$B_{0}\rightarrow \infty$
, the incoming flow is strongly guided by the mean field, and will ignore the weaker underlying turbulence. Hence, MHD prescriptions are violated. When does the switch from one regime to the other happen? We conjecture that particles will tend to follow the mean field instead of being randomized whenever the energy density
$B_{0}^{2}/8\unicode[STIX]{x03C0}$
of the mean field exceeds a fraction of order unity of the turbulent energy
${\mathcal{E}}_{T}$
. Now, if the turbulence is caused by an instability of the counter-streaming flows, its energy will be a fraction of the flow energy density, i.e.
${\mathcal{E}}_{T}\lesssim \unicode[STIX]{x1D6FE}_{0}n_{0}mc^{2}$
. As a consequence, the system will depart from MHD beyond a critical value of
$(B_{0}^{2}/8\unicode[STIX]{x03C0})/\unicode[STIX]{x1D6FE}_{0}n_{0}mc^{2}=\unicode[STIX]{x1D70E}/2$
. We thus conclude that, regardless of which instability is initially triggered in the overlapping region, the MHD behaviour is inhibited for values of
$\unicode[STIX]{x1D70E}$
greater than about unity. This is indeed what is observed in our PIC simulations.
6 Conclusions
In summary, we have found a departure from MHD behaviour when two collisionless pair plasma shells with a flow-aligned magnetic field collide. While MHD stipulates that the very same shock should form regardless of the
$\unicode[STIX]{x1D70E}$
parameter, the micro-physics analysis of the shock formation allows to understand why the standard shock formation scenario can be jeopardized beyond a critical magnetization.
PIC simulations have confirmed the theoretical analysis. The results are similar when considering an angle
$\unicode[STIX]{x1D703}=5^{\circ }$
between the field and the flow (see appendix B). This shows that the observed MHD departure is not a ‘Dirac delta’ effect, strictly restricted to
$\unicode[STIX]{x1D703}=0$
.
What about an electron/proton plasma? It is difficult at this stage to draw definite conclusions about that case. When protons are accounted for instead of positrons, the asymmetric role of electrons and protons results in an upstream current which, in the presence of a flow-aligned magnetic field, is likely to trigger the Bell instability (Bell Reference Bell2004). This instability is not triggered here because of the symmetric role of electrons and positrons. But if excited, the upstream Bell turbulence, when transported downstream, could help isotropizing the flow. Yet, in spite of some differences with pair plasmas (Stockem Novo et al.
Reference Stockem Novo, Bret, Fonseca and Silva2015), shock formation in electron/proton plasmas eventually still boils down to the capability of an instability generated turbulence to stop the flow. If the conjecture enounced at the end of § 5 turns out to be valid, we could recover a
$\unicode[STIX]{x1D70E}$
threshold for the validity of MHD in electron/proton plasmas as well, since the energy of the downstream turbulence should remain a fraction of the upstream kinetic energy. Further studies will be necessary to sort out this important issue.
Would it be possible to modify MHD so that it keeps fitting the kinetic results for
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
? A tentative pathway, beyond the scope of this work, would be to include the downstream anisotropy within the MHD analysis. Indeed, bottom panels of figures 3 and 4 clearly show that the downstream is not isotropized because of the magnetic field. One could therefore try to quantify this anisotropy in terms of the field before inserting the corresponding temperature anisotropy in the Rankine–Hugoniot jump conditions analysis (Karimabadi, Krauss-Varban & Omidi Reference Karimabadi, Krauss-Varban and Omidi1995; Vogl et al.
Reference Vogl, Biernat, Erkaev, Farrugia and Mühlbachler2001; Gerbig & Schlickeiser Reference Gerbig and Schlickeiser2011).
Future work will also explore in detail the angular dependence of our results, together with the expected consequences for astrophysics.
Acknowledgements
A.B. acknowledges grants ENE2013-45661-C2-1-P, PEII-2014-008-P and ANR-14-CE33-0019 MACH. A.P. acknowledges support by the European Union Seventh Framework Program (FP7/2007-2013) under grant agreement no. 618499, and support from NASA under grant no. NNX12AO83G. RN’s research was supported in part by NASA grant TCAN NNX14AB47G. OS acknowledges support by NASA through Einstein Post-doctoral Fellowship number PF4-150126 awarded by the Chandra X-ray Center, operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. Thanks to S. Naoz and V. Malka for valuable inputs.
Appendix A. MHD predictions
A.1 MHD criterion for slow shocks
A key quantity is the speed of the shock front relative to the phase velocity of the fast mode. Consider first the upstream fast mode. Its phase velocity
$v_{fast}$
, which is also the phase velocity of the Alfvén mode, is given in the upstream frame by (Kirk & Duffy Reference Kirk and Duffy1999; Keppens & Meliani Reference Keppens and Meliani2008)Footnote
5
,
$$\begin{eqnarray}v_{fast}^{2}=c^{2}\frac{\unicode[STIX]{x1D70E}^{\prime }}{1+\unicode[STIX]{x1D70E}^{\prime }},\quad \text{with }\unicode[STIX]{x1D70E}^{\prime }=\frac{B_{0}^{2}}{4\unicode[STIX]{x03C0}n_{0}^{\prime }mc^{2}},\end{eqnarray}$$
where the prime stands for quantities measured in the upstream frame. Note that the field has no prime because it is aligned with the axis of the Lorentz transformations.
Regarding the velocity of the shock front, it is in two dimensions and for large
$\unicode[STIX]{x1D6FE}_{0}$
,
$\unicode[STIX]{x1D6FD}_{s}\equiv v_{s}/c=1/3$
in the downstream frame (Kirk & Duffy Reference Kirk and Duffy1999). In this frame, the upstream propagates at
$\boldsymbol{v}_{0}$
(figure 1). Therefore, the speed
$v_{s}^{\prime }$
of the shock front in the upstream frame is
$$\begin{eqnarray}\frac{v_{s}^{\prime }}{c}\equiv \unicode[STIX]{x1D6FD}_{s}^{\prime }=\frac{\unicode[STIX]{x1D6FD}_{s}+\unicode[STIX]{x1D6FD}_{0}}{1+\unicode[STIX]{x1D6FD}_{s}\unicode[STIX]{x1D6FD}_{0}}\sim 1-\frac{1}{4\unicode[STIX]{x1D6FE}_{0}^{2}}\quad \text{if }\unicode[STIX]{x1D6FE}_{0}\gg 1.\end{eqnarray}$$
The condition
$v_{fast}>v_{s}^{\prime }$
reads therefore (
$\unicode[STIX]{x1D6FE}_{0}\gg 1$
),
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x1D70E}^{\prime }}{1+\unicode[STIX]{x1D70E}^{\prime }}>\left(1-\frac{1}{4\unicode[STIX]{x1D6FE}_{0}^{2}}\right)^{2}\\ \displaystyle \Rightarrow \unicode[STIX]{x1D70E}^{\prime }>2\unicode[STIX]{x1D6FE}_{0}^{2}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
We need now to express the left-hand side of this equation in terms of the
$\unicode[STIX]{x1D70E}$
parameter (2.1). Since the field
$\boldsymbol{B}_{0}$
is aligned with the axis of the Lorenz transformation, it does not change. Regarding the density, there is a relativistic bunching between
$n_{0}^{\prime }$
, the upstream density in the upstream frame, and
$n_{0}$
, the upstream density in the downstream frame, with
$$\begin{eqnarray}n_{0}=\unicode[STIX]{x1D6FE}_{0}n_{0}^{\prime }.\end{eqnarray}$$
Combining with (2.1), (A 3), we finally obtain,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}>2.\end{eqnarray}$$
Therefore, MHD shocks with
$\unicode[STIX]{x1D70E}>2$
have to be of the slow type. This limit is pictured on figure 2 in the
$(\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6FE}_{0})$
phase space, together with a refined calculation accounting for the
$\unicode[STIX]{x1D6FE}_{0}$
dependence of the threshold.
Our calculation has been conducted in two dimensions in order to compare with the PIC simulations. In three dimensions, one has
$\unicode[STIX]{x1D6FD}_{s}=1/4$
and (A 5) reads
$\unicode[STIX]{x1D70E}>5/3$
instead. We now turn to the same analysis, but for the downstream.
A.2 MHD criterion for extraneous shocks
Measured in the downstream, the downstream Alfvén and fast mode velocities, read
$$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{fast}^{2}=\frac{\unicode[STIX]{x1D70E}_{w}}{1+\unicode[STIX]{x1D70E}_{w}},\end{eqnarray}$$
with,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{w}=\frac{B_{0}^{2}}{4\unicode[STIX]{x03C0}w},\quad \text{and}\quad w=nmc^{2}+\frac{\hat{\unicode[STIX]{x1D6FE}}}{\hat{\unicode[STIX]{x1D6FE}}-1}nk_{B}T,\end{eqnarray}$$
where
$\hat{\unicode[STIX]{x1D6FE}}=4/3$
is the adiabatic index,
$w$
the enthalpy density and
$n$
the downstream density in its own frame. Writing (Service Reference Service1986),
$$\begin{eqnarray}nk_{B}T={\textstyle \frac{4}{3}}\unicode[STIX]{x1D6FE}_{0}^{2}n_{0}mc^{2},\end{eqnarray}$$
and neglecting
$nmc^{2}$
in the enthalpy, we obtain
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{w}=\frac{3}{16}\frac{B_{0}^{2}/4\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{0}^{2}n_{0}mc^{2}}.\end{eqnarray}$$
Since in the downstream frame, the shock front propagates at
$c/3$
, we need to compare
$\unicode[STIX]{x1D6FD}_{fast}^{2}$
with
$1/9$
. The equation
$\unicode[STIX]{x1D6FD}_{fast}^{2}=1/9$
gives,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{w}=\frac{1}{8}\Rightarrow \frac{B_{0}^{2}/4\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{0}^{2}n_{0}mc^{2}}=\unicode[STIX]{x1D70E}=\frac{2}{3}.\end{eqnarray}$$
The corresponding
$\unicode[STIX]{x1D70E}$
limit is pictured on the same figure 2. Table 2 summarizes the MHD conditions for super-Alfvénic flow upstream and downstream, (A 5), (A 10). As it appears, the interval
$2/3<\unicode[STIX]{x1D70E}<2$
defines a range of extraneous solutions. For
$\unicode[STIX]{x1D70E}<2/3$
, MHD gives a fast shock and for
$\unicode[STIX]{x1D70E}>2$
, it gives a slow shock.
Table 2. Summary of the MHD conditions for super-Alfvénic gflow upstream and downstream for
$\unicode[STIX]{x1D6FE}_{0}\gg 1$
, (A 5), (A 10). The interval
$2/3<\unicode[STIX]{x1D70E}<2$
pertains to extraneous solutions.
Figure 5. MHD simulations of the collision of two cold pair plasmas with
$\unicode[STIX]{x1D6FE}_{0}=10$
. The extraneous regime requires a tilted field in order to appear.
A.3 MHD simulations
We performed 1-D MHD simulations of colliding magnetized flows in order to confirm the previous analysis. We used the general relativistic magnetohydrodynamical code KORAL (Sa̧dowski et al.
Reference Sa̧dowski, Narayan, McKinney and Tchekhovskoy2014). The domain was filled initially with gas of uniform comoving frame density,
$\unicode[STIX]{x1D70C}=1$
, and velocities set-up so that the gas collides at
$x=0$
and on both sides of the contact surface the velocities correspond to the given Lorentz factor
$\unicode[STIX]{x1D6FE}_{0}$
. The gas magnetization is described by the
$\unicode[STIX]{x1D70E}$
parameter and tilt angle
$\unicode[STIX]{x1D703}$
(zero tilt angle corresponds to field perfectly parallel to the flow). We ran a few simulations for the flow-aligned case (
$\unicode[STIX]{x1D703}=0$
) as well as for a
$5^{\circ }$
tilt. The results are shown in figure 5. For the flow-aligned field, the exact decoupling of the fluid motion from the field (Majorana & Anile Reference Majorana and Anile1987) results in a shock independent from
$\unicode[STIX]{x1D70E}$
. For the
$5^{\circ }$
tilt, the field couples to the fluid and the extraneous regime is unravelled as the shock splits into 2 sub-shocks (Kulsrud Reference Kulsrud2005) for intermediate
$\unicode[STIX]{x1D70E}$
parameters.
Appendix B. PIC simulations
B.1 Simulation set-up
The shock is set-up by reflecting a cold ‘upstream’ flow from a conducting wall located at
$x=0$
. The interaction between the incoming beam (that propagates along
$-\boldsymbol{e}_{x}$
) and the reflected beam triggers the formation of a shock, which moves away from the wall along
$+\boldsymbol{e}_{x}$
. This set-up is equivalent to the head-on collision of two identical plasma shells (figure 1), which would form a forward and reverse shock and a contact discontinuity. Here, we follow only one of these shocks, and replace the contact discontinuity with the conducting wall. The simulation is performed in the ‘wall’ frame, where the ‘downstream’ plasma behind the shock is at rest.
We use a rectangular simulation box in the
$x,y$
plane, with periodic boundary conditions in the
$y$
direction. The incoming plasma is injected through a ‘moving injector’, which recedes from the wall along
$+\boldsymbol{e}_{x}$
at the speed of light. The simulation box is expanded in the
$x$
direction as the injector approaches the right boundary of the computational domain. This permits us to save memory and computing time, while following the evolution of all the upstream regions that are causally connected with the shock.
In two dimensions, each computational cell is initialized with 16 electrons and 16 positrons. The relativistic electron skin depth for the incoming plasma (
$c/\unicode[STIX]{x1D714}_{p}$
, with
$\unicode[STIX]{x1D714}_{p}^{2}=4\unicode[STIX]{x03C0}n_{0}q^{2}/\unicode[STIX]{x1D6FE}_{0}m$
) is resolved with 10 computational cells and the simulation time step is
$\unicode[STIX]{x0394}t=0.045\unicode[STIX]{x1D714}_{p}^{-1}$
. The computational domain is typically
${\sim}102c/\unicode[STIX]{x1D714}_{p}$
wide (corresponding to 1024 cells). The simulations extend typically up to
$\unicode[STIX]{x1D714}_{p}t\sim 3600$
, corresponding to a box length of
${\sim}3600c/\unicode[STIX]{x1D714}_{p}$
, or
${\sim}36\,000$
cells.
The incoming stream is injected along
$-\boldsymbol{e}_{x}$
with bulk Lorentz factor
$\unicode[STIX]{x1D6FE}_{0}$
. The incoming plasma is cold, with thermal spread
$k_{B}T/mc^{2}\ll 1$
. We vary the upstream Lorentz factor between
$\unicode[STIX]{x1D6FE}_{0}=10$
and 30, but we find that our results are nearly insensitive to the choice of the Lorentz factor, as long as the flow is ultra-relativistic. The upstream flow is seeded with a uniform magnetic field
$\boldsymbol{B}_{0}$
.
Figure 6. Same as figure 3, but for
$\unicode[STIX]{x1D6FE}_{0}=30$
.
Figure 6 shows the result of a series of simulations similar to those displayed on figure 3, but for
$\unicode[STIX]{x1D6FE}_{0}=30$
. The weak
$\unicode[STIX]{x1D6FE}_{0}$
-dependence, expected from the micro-physics analysis of § 5, is confirmed.
B.2 Oblique shocks
In most of our studies, the upstream field is aligned with the flow (i.e. we study a ‘parallel shock’), but we also consider quasi-parallel shocks in which the upstream field makes an angle
$\unicode[STIX]{x1D703}=5^{\circ }$
with the flow direction of propagationFootnote
6
.
In figure 7, we present the results of a suite of 2-D simulations of nearly parallel shocks, with obliquity
$\unicode[STIX]{x1D703}=5^{\circ }$
. In addition to the
$y$
-averaged profile of the particle number density (figure 7
a) and the measure of particle anisotropy (figure 7
b), we present the
$y$
-averaged profiles of the transverse magnetic field
$B_{z}$
, in units of the upstream field
$B_{0}$
. For such small field obliquity, most of the conclusions presented above for the case of parallel shocks still hold. In particular, for
$\unicode[STIX]{x1D70E}\gtrsim 0.6$
the downstream flow is no longer isotropic. As a result, the density jump is lower than the MHD predictions.
Figure 7. Shock structure from a series of 2-D PIC simulations with
$\unicode[STIX]{x1D6FE}_{0}=10$
,
$0.2\leqslant \unicode[STIX]{x1D70E}\leqslant 3$
, at
$\unicode[STIX]{x1D714}_{p}t=450$
. We plot the
$y$
-averaged density profile (a), the
$y$
-averaged profile of the transverse field
$B_{z}$
(b), and a measure of the plasma anisotropy (c), as defined in the text. The vertical dashed line indicates the position of the shock front, assuming that it propagates at
$c/3$
. The angle between the flow and the field is
$\unicode[STIX]{x1D703}=5^{\circ }$
.
Appendix C. Value of the parameter
$\unicode[STIX]{x1D705}$
in (5.5)
The criterion (5.5) has been obtained in Bret (Reference Bret2016b
) assuming equipartition at saturation of the Weibel instability. We here elaborate on the parameter
$\unicode[STIX]{x1D705}$
beyond what is done in Bret (Reference Bret2016b
).
For the present system, the growth rate
$\unicode[STIX]{x1D6FF}_{W}$
of the Weibel instability is given by (4.2). The instability grows the field given by (5.2) until (Davidson et al.
Reference Davidson, Hammer, Haber and Wagner1972),
$$\begin{eqnarray}\frac{q\sqrt{\langle \boldsymbol{B}_{f}^{2}\rangle }}{\unicode[STIX]{x1D6FE}_{0}mc}=\frac{qB_{f}/\sqrt{2}}{\unicode[STIX]{x1D6FE}_{0}mc}=\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FF}_{W},\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}=1$
means equipartition, as evidenced by (C 4), (C 5) belowFootnote
7
. The magnetic filaments become unable to trap the particles, regardless of their initial
$y$
position or velocities, for
$$\begin{eqnarray}B_{0}>{\textstyle \frac{1}{2}}B_{f}.\end{eqnarray}$$
We now replace
$B_{f}$
in the equation above by its expression from (C 2). We then express
$\unicode[STIX]{x1D6FF}_{W}$
from (4.2) and find that
$B_{0}>B_{f}/2$
is equivalent to,
$$\begin{eqnarray}\unicode[STIX]{x1D70E}>\frac{\unicode[STIX]{x1D702}^{2}}{1+\unicode[STIX]{x1D702}^{2}/2}\unicode[STIX]{x1D6FD}_{0}^{2}\equiv \unicode[STIX]{x1D705}\unicode[STIX]{x1D6FD}_{0}^{2}.\end{eqnarray}$$
The vertical asymptote defined by the blue curve on figure 2 relies therefore on the parameter
$\unicode[STIX]{x1D702}$
. The constraints on its value stem from the degree of equipartition reached at saturation of the Weibel instability. In this respect, let us compute the total amount of magnetic energy contained in the magnetic field at saturation, and compare it to the upstream kinetic energy. Using (C 1) to express
$B_{f}^{2}$
we find,
$$\begin{eqnarray}\frac{(B_{f}^{2}/2+B_{0}^{2})/8\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{0}n_{0}mc^{2}}=\unicode[STIX]{x1D702}^{2}\unicode[STIX]{x1D6FD}_{0}^{2}+\frac{\unicode[STIX]{x1D70E}}{2}(1-\unicode[STIX]{x1D702}^{2}).\end{eqnarray}$$
Accounting only for the energy contained in the Weibel field, one finds,
$$\begin{eqnarray}\frac{(B_{f}^{2}/2)/8\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FE}_{0}n_{0}mc^{2}}=\unicode[STIX]{x1D702}^{2}(\unicode[STIX]{x1D6FD}_{0}^{2}-\unicode[STIX]{x1D70E}/2).\end{eqnarray}$$
In the relativistic regime where
$\unicode[STIX]{x1D6FD}_{0}\sim 1$
, (C 4) gives unity for
$\unicode[STIX]{x1D702}=1$
and is monotonically increasing for
$\unicode[STIX]{x1D702}\in [0,1]$
and
$\unicode[STIX]{x1D70E}<2$
. Note that the regime
$\unicode[STIX]{x1D70E}>2\unicode[STIX]{x1D6FD}_{0}^{2}$
is irrelevant in the present context since (4.2) implies the instability vanishes there. (C 5) implies that the relative amount of energy contained in the Weibel field varies like
$\unicode[STIX]{x1D702}^{2}$
. For
$\unicode[STIX]{x1D702}=1$
, (C 3) give
$\unicode[STIX]{x1D705}=2/3$
, which is the values considered on figure 1.
