Kinetic inhibition of MHD-shocks in the vicinity of a parallel magnetic field

TL;DR

This paper demonstrates that a strong flow-aligned magnetic field can inhibit shock formation in collisionless pair plasmas, challenging traditional MHD predictions and revealing new plasma behavior.

Contribution

It shows that magnetic fields aligned with plasma flow can suppress shock formation, a violation of classical MHD theory, supported by micro-physics analysis and particle-in-cell simulations.

Findings

Magnetic field suppresses isotropization in plasma

Shock formation is inhibited by strong flow-aligned fields

Results extend beyond strictly parallel magnetic configurations

Abstract

According to magnetohydrodynamics (MHD), the encounter of two collisional magnetized plasmas at high velocity gives rise to shock waves. Investigations conducted so far have found that the same conclusion still holds in the case of collisionless plasmas. For the case of a flow-aligned field, MHD stipulates that the field and the fluid are disconnected, so that the shock produced is independent of the field. We present a violation of this MHD prediction when considering the encounter of two cold pair plasmas along a flow-aligned magnetic field. As the guiding magnetic field grows, isotropization is progressively suppressed, resulting in a strong influence of the field on the resulting structure. A micro-physics analysis allows to understand the mechanisms at work. Particle-in-cell simulations also support our conclusions and show that the results are not restricted to a strictly parallel field.

Tables

Table 1. Table 1: Values of the parameter Ω B = 2 ​ γ 0 ​ σ subscript Ω 𝐵 2 subscript 𝛾 0 𝜎 \Omega_{B}=\sqrt{2\gamma_{0}\sigma} used in Bret ( 2016 a ) corresponding to the σ 𝜎 \sigma ’s sampled here and for γ 0 = 10 subscript 𝛾 0 10 \gamma_{0}=10 .

$\sigma$	0.2	0.4	0.6	0.8	1	1.5	2	3
${\Omega }_B$	2.0	2.8	3.5	4.0	4.5	5.5	6.3	7.7

Table 2. Table 2: Summary of the MHD conditions for super-Alfvénic upstream and downstream for γ 0 ≫ 1 much-greater-than subscript 𝛾 0 1 \gamma_{0}\gg 1 , Eq. ( 13 , 18 ). The interval 2 / 3 < σ < 2 2 3 𝜎 2 2/3<\sigma<2 pertains to extraneous solutions.

$\sigma$			$>2$
Upstream	Super-A	Super-A	Sub-A
Downstream	Super-A	Sub-A	Sub-A
Shock type	Fast	Extraneous	Slow

Equations

$$\sigma =\frac{B_0^2/4\pi }{{\gamma }_0{n}_{0}m{c}^2},$$

$$\phi =\frac{\mathrm{Var}(p_y)+\mathrm{Var}(p_z)}{2\mathrm{Var}(p_x)}.$$

$${\delta }_W={\omega }_p\sqrt{2{\beta }_0^2-\sigma },$$

$${\gamma }_0>\sqrt{\frac{2}{4/3-\sigma }}.$$

$${\mathbf{B}}_f=B_f\sin (ky){\mathbf{e}}_z,$$

$$k^{-1}>\frac{v_0}{{\omega }_{B_f}},\mathrm{with}{\omega }_{B_f}=\frac{q{B}_f}{{\gamma }_{0}mc},$$

$$B_0>\frac{1}{2}{B}_f.$$

$$\sigma >\kappa {\beta }_0^2,$$

$$v_{\mathrm{fast}}^2=c^2\frac{{\sigma }^{\prime }}{1+{\sigma }^{\prime }},\mathrm{with}{\sigma }^{\prime }=\frac{B_0^2}{4\pi n_0^{\prime }m{c}^2},$$

$$\frac{v_s^{\prime }}{c}\equiv {\beta }_s^{\prime }=\frac{{\beta }_s+{\beta }_0}{1+{\beta }_s{\beta }_0}\sim 1-\frac{1}{4{\gamma }_0^2}\mathrm{if}{\gamma }_0\gg 1.$$

$$\frac{{\sigma }^{\prime }}{1+{\sigma }^{\prime }}$$

$$\Rightarrow {\sigma }^{\prime }$$

$$n_0={\gamma }_0{n}_0^{\prime }.$$

$$\sigma >2.$$

$${\beta }_{\mathrm{fast}}^2=\frac{{\sigma }_w}{1+{\sigma }_w},$$

$${\sigma }_w=\frac{B_0^2}{4\pi w},\mathrm{and}w=nm{c}^2+\frac{\hat{\gamma }}{\hat{\gamma }-1}n{k}_{B}T,$$

$$n{k}_{B}T=\frac{4}{3}{\gamma }_0^2{n}_{0}m{c}^2,$$

$${\sigma }_w=\frac{3}{16}\frac{B_0^2/4\pi }{{\gamma }_0^2{n}_{0}m{c}^2}.$$

$${\sigma }_w=\frac{1}{8}\Rightarrow \frac{B_0^2/4\pi }{{\gamma }_0^2{n}_{0}m{c}^2}=\sigma =\frac{2}{3}.$$

$$\frac{q\sqrt{<{\mathbf{B}}_f^2>}}{{\gamma }_{0}mc}=\frac{q{B}_f/\sqrt{2}}{{\gamma }_{0}mc}=\eta {\delta }_W,$$

$$B_0>\frac{1}{2}{B}_f.$$

$$\sigma >\frac{{\eta }^2}{1+{\eta }^2/2}{\beta }_0^2\equiv \kappa {\beta }_0^2.$$

$$\frac{(B_f^2/2+B_0^2)/8\pi }{{\gamma }_0{n}_{0}m{c}^2}={\eta }^2{\beta }_0^2+\frac{\sigma }{2}(1-{\eta }^2).$$

$$\frac{(B_f^2/2)/8\pi }{{\gamma }_0{n}_{0}m{c}^2}={\eta }^2({\beta }_0^2-\sigma /2).$$

Full text

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, 2002).

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 et al., 2003; Schwartz et al., 2011). 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, 1958; Sagdeev, 1966).

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, 1987; Sironi et al., 2015; Marcowith et al., 2016). They are also believed to play a role in the generation of gamma-ray-bursts (Mészáros & Rees, 2014; Pe’er, 2015) and fast radio bursts (Lyubarsky, 2014; Falcke & Rezzolla, 2014).

Starting with the pioneering work of Sagdeev in the 1960’s (Sagdeev, 1966), 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, 2005; Martins et al., 2009). However, as recently as the 1990’s, there were still doubts about the very existence of collisionless shocks (Sagdeev & Kennel, 1991). 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 shocks111Before they collide, two plasmas display a Debye sheath at their border, with an associated potential jump (Gurnett & Bhattacharjee, 2005). At low energy of collision, the encounter is mediated by the interaction of these sheaths, and an electrostatic shock is formed. At higher energy, the interaction is rather mediated by the counter-streaming instabilities arising from the overlapping of the plasmas (Stockem et al., 2014; Bret, 2015). If the dominant instability is the Weibel one (see conditions in Bret et al. (2010)), then a “Weibel shock” is formed. have been observed in the laboratory (Ahmed et al., 2013), while the production of Weibel mediated shocks such as the ones discussed here, is expected within the next few years (Huntington et al., 2015; Lobet et al., 2015; Park et al., 2016). Note that the “Weibel instability” we refer to is sometimes labelled “filamentation instability” of “beam Weibel” instability (Silva et al., 2002; Hill et al., 2005; Deutsch et al., 2005). 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 (2015) 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 approach222Once source of discrepancies are the accelerated particles which escape the Rankine-Hugoniot “budget” (Stockem et al., 2012; Sironi et al., 2013; Caprioli & Spitkovsky, 2014; Bret, 2015).. For the case of a flow-aligned field, MHD prescribes that the fluid and the field are decoupled (Majorana & Anile, 1987), 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 section 3, we explain the predictions of MHD for this system. Then, in section 4, we describe a series of simulations using the particle-in-cell (PIC) technique. These simulations work at the microscopic level, and show a departure from the MHD predictions beyond a critical magnetization. In section 5, we present a micro-physics analysis of the shock formation process explaining the departure from MHD.

2 System considered

The system considered is shown schematically in Fig. 1. Two identical pair plasma shells of density $n_{0}$ head toward each other with initial velocity $\pm\textbf{v}_{0}$ and Lorentz factor $\gamma_{0}=(1-v_{0}^{2}/c^{2})^{-1/2}$. The whole system is embedded in an external field $\mathbf{B}_{0}\parallel\mathbf{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 0. 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,

[TABLE]

where all quantities are measured in the downstream frame.

3 MHD predictions

An MHD plasma sustains 3 kinds of modes: slow mode, Alfvén mode, and fast mode (Kulsrud, 2005). The phase velocities of these modes satisfy the hierarchy $v_{\rm slow}<v_{\rm Alfven}<v_{\rm fast}$. Because of this hierarchy, the Alfvén mode is sometimes dubbed the “intermediate mode” (Kulsrud, 2005). In the cold limit considered here, $v_{\rm slow}\rightarrow 0$ and $v_{\rm fast}\rightarrow v_{\rm 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 eg Kirk & Duffy (1999)). However, such solutions of the MHD jump equations do not survive when produced and are called “extraneous” (Kulsrud, 2005); 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, 1987). The shock formed is therefore the same, regardless of the magnetization parameter $\sigma$. Nevertheless, its front velocity can still be compared to the phase speeds of the three modes. In the present cold limit, and for $\gamma_{0}\to\infty$, the shock is expected to be “fast” for $\sigma<2/3$, “slow” for $\sigma>2$, and “extraneous” in between (see Appendix A). Figure 2 shows these limits for a range of $\gamma_{0}$ in the $(\sigma,\gamma_{0})$ plane.

The MHD predictions for the present system are therefore very clear: the same shock should form regardless of the $\sigma$ parameter, simply because the fluid and the field are perfectly decoupled here. The MHD simulations run in Appendix A confirm this conclusion.

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 3D electromagnetic PIC code TRISTAN-MP (Spitkovsky, 2005), which is a parallel version of the publicly available code TRISTAN (Buneman, 1993) that has been optimized for studying relativistic collisionless shocks (Spitkovsky, 2008a, b; Sironi & Spitkovsky, 2009, 2011; Sironi et al., 2013). We employ simulations in 2D computational domains, but all three components of particle velocities and electromagnetic fields are tracked (see more details in Appendix B).

We probe the regime $\gamma_{0}=10$ and $0<\sigma<3$. Note that the parameter space in Bret (2016a) is parameterized is terms of $\gamma_{0}$ and $\Omega_{B}$, the later being related to the present $\sigma$ through $\Omega_{B}=\sqrt{2\gamma_{0}\sigma}$333The factor 2 comes from the fact that the plasma frequency in Bret (2016a) is the one of the electrons (or the positrons) alone, while the density $n_{0}$ in $\sigma$ is the total density of one pair beam.. For better clarity, Table 1 gives the values of the $\Omega_{B}$’s of Bret (2016a) corresponding to the $\sigma$’s sampled here.

Figure 3 shows the $y$-integrated density profile of the system for $\gamma_{0}=10$ (top panel), at a relatively early time, $\omega_{p}t=450$, where $\omega_{p}^{2}=4\pi n_{0}q^{2}/\gamma_{0}m$. The magnetization parameter varies from 0.2 to 3, as indicated in the legend. In the bottom panel, we quantify the isotropization of the particle distribution function by plotting the ratio $\varphi$ between the momentum dispersion along the transverse directions ($y$ and $z$) as compared to the longitudinal direction $x$, namely,

[TABLE]

We notice that, for $\sigma\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 $\gamma_{0}=10$, and $\sim 4$ in the limit $\gamma_{0}\gg 1$). Correspondingly, the shock speed approaches the value $\sim c/3$ predicted by MHD (indicated by the dashed line). The bottom panel in Fig. 3 shows that for low magnetizations the downstream plasma is nearly isotropic. However, for higher magnetizations ($\sigma\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$.

Noteworthily, the width of the density jump increases notably with $\sigma$. For small values, the shock front is $\sim 70c/\omega_{p}$ thick. But for $\sigma=3$, the transition region between the “upstream” and the “downstream” is $\sim 300c/\omega_{p}$.

Why do the results for $\sigma\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 $\delta_{W}$ of the Weibel instability is given by (Stockem et al., 2006; Bret, 2016a),

[TABLE]

where $\beta_{0}=v_{0}/c$. The shock formation time typically amounts to a few tens of $e$-folding times (Bret et al., 2013, 2014). With the parameters used here, $20\delta_{W}^{-1}$ is at most $28\omega_{p}^{-1}$ for $\sigma=1.5$ ($\delta_{W}$ vanishes for $\sigma>2\beta_{0}^{2}$). Therefore, the time $t=450\omega_{p}^{-1}$ to which the simulations in Fig. 3 have been run, exceeds by a factor of 15 the slowest expected shock formation time.

To verify the above argument, we have evolved the simulations to much longer times: $\omega_{p}t=3600$ (Fig. 4). We again find that the density profile strongly varies with $\sigma$, contrary to the MHD prescriptions. For magnetizations $\sigma\gtrsim 0.6$, the system settles in 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 the bottom panel of Figs. 3 & 4 , for $\sigma\gtrsim 0.6$, the downstream particle distribution is hotter along the longitudinal direction than in the transverse directions444Note that, since the particle momenta are measured in the downstream frame of the simulations, we do not expect $\varphi=1$ in the upstream medium (but rather $\varphi\propto 1/\gamma_{0}$), despite the fact that the upstream plasma is isotropic in its own rest frame.. For large $\sigma$, we find downstream $\varphi<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 $\sigma$, the shock front on Fig. 4 is still $\sim 70c/\omega_{p}$ thick. But for $\sigma=3$, the transition region is now $\sim 2000c/\omega_{p}$.

The micro-physical analysis discussed next in section 5 predicts that the departure from the MHD behavior we just observed, is $\gamma_{0}$-independent at large $\gamma_{0}$. This prediction has been successfully tested in Appendix B by running a series of PIC simulations with $\gamma_{0}=30$. We also confirm in Appendix B that these results are not restricted to a perfectly flow-aligned field ($\theta=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., 2010), but the Weibel (filamentation) instability, with a $\mathbf{k}$ normal to the flow, grows faster than all others, provided (Bret, 2016a),

[TABLE]

The line corresponding to this limit is shown in Fig. 2 by the orange curve. The Weibel instability dominates the unstable spectrum of the system for all points of the $(\sigma,\gamma_{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., 2013, 2014). 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,

[TABLE]

where $k$ is the fastest growing wave-number. 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, 2015) shows that it is stopped inside if

[TABLE]

where $\mathbf{v}_{0}$ is the initial, flow-aligned velocity of the particle and $\gamma_{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., 2014).

How is Eq. (6) 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, 2016b)

[TABLE]

Since $B_{f}$ arises from the growth of the Weibel instability, its magnitude can be quantified (Stockem et al., 2006; Bret, 2016a). Therefore, the above criterion can eventually be expressed in terms of $\sigma$ and $\beta_{0}^{2}=1-1/\gamma_{0}^{2}$, giving (see details in Appendix C),

[TABLE]

where $\kappa=2/3$ if equipartition is assumed at saturation of the Weibel instability. The boundary corresponding to the criterion (8) with $\kappa=2/3$ is shown in Figure 2 by the blue curve.

The region between the bounds corresponding to Eqs. (4) and (8), i.e., the region between the orange and blue lines in Fig. 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 my have noticed that for $\sigma=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 Eq. (5), consider a turbulent electromagnetic perturbation $\sum_{\mathbf{k}}\mathbf{E}_{\mathbf{k}}+\mathbf{B}_{\mathbf{k}}$ (with $<\mathbf{E}_{\mathbf{k}}>=<\mathbf{B}_{\mathbf{k}}>=\mathbf{0}$) that is present in the overlapping region and that can potentially isotropize the incoming flow. Consider also a superimposed, flow-aligned field $\mathbf{B}_{0}$. In the limit $B_{0}=0$, the incoming flow is isotropized, and usual MHD applies. In the opposite limit $B_{0}\to\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\pi$ 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\gamma_{0}n_{0}mc^{2}$. As a consequence, the system will depart from MHD beyond a critical value of $(B_{0}^{2}/8\pi)/\gamma_{0}n_{0}mc^{2}=\sigma/2$. We thus conclude that, regardless of which instability is initially triggered in the overlapping region, the MHD behaviour is inhibited for values of $\sigma$ 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 $\sigma$ 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 $\theta=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 $\theta=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, 2004). 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., 2015), 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 Section 5 turns out to be valid, we could recover a $\sigma$ 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 $\sigma\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-Figs. 3 & 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 et al., 1995; Vogl et al., 2001; Gerbig & Schlickeiser, 2011).

Future work will also explore in detail the angular dependence of our results, together with the expected consequences for astrophysics.

7 Acknowledgements

AB acknowledges grants ENE2013-45661-C2-1-P, PEII-2014-008-P and ANR-14-CE33-0019 MACH. AP acknowledges support by the European Union Seventh Framework Program (FP7/2007-2013) under grant agreement #618499, and support from NASA under grant #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 Smadar Naoz and Victor 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, 1999; Keppens & Meliani, 2008)555The following expression allows for the relativistic corrections needed when the field energy density $B_{0}^{2}/8\pi$ approaches the matter energy density $n_{0}mc^{2}$. In such circumstances, the field participates in the inertia of the medium and modifies the speed of the waves.,

[TABLE]

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 2D and for large $\gamma_{0}$, $\beta_{s}\equiv v_{s}/c=1/3$ in the downstream frame (Kirk & Duffy, 1999). In this frame, the upstream propagates at $\mathbf{v}_{0}$ (Fig. 1). Therefore, the speed $v_{s}^{\prime}$ of the shock front in the upstream frame is

[TABLE]

The condition $v_{\rm fast}>v_{s}^{\prime}$ reads therefore ($\gamma_{0}\gg 1$),

[TABLE]

We need now to express the left-hand-side of this equation in terms of the $\sigma$ parameter Eq. (1). Since the field $\mathbf{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

[TABLE]

Combining with Eqs. (1,A.1), we finally obtain,

[TABLE]

Therefore, MHD shocks with $\sigma>2$ have to be of the slow type. This limit is pictured on Fig. 2 in the $(\sigma,\gamma_{0})$ phase space, together with a refined calculation accounting for the $\gamma_{0}$-dependence of the threshold.

Our calculation has been conducted in 2D in order to compare with the PIC simulations. In 3D, one has $\beta_{s}=1/4$ and Eq. (13) reads $\sigma>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

[TABLE]

with,

[TABLE]

where $\hat{\gamma}=4/3$ is the adiabatic index, $w$ the enthalpy density, and $n$ the downstream density in its own frame. Writing (Service, 1986),

[TABLE]

and neglecting $nmc^{2}$ in the enthalpy, we obtain

[TABLE]

Since in the downstream frame, the shock front propagates at $c/3$, we need to compare $\beta_{\rm fast}^{2}$ with $1/9$. The equation $\beta_{\rm fast}^{2}=1/9$ gives,

[TABLE]

The corresponding $\sigma$-limit is pictured on the same Fig. 2. Table 2 summarizes the MHD conditions for super-Alfvénic upstream and downstream, Eqs. (13,18). As it appears, the interval $2/3<\sigma<2$ defines a range of extraneous solutions. For $\sigma<2/3$, MHD gives a fast shock and for $\sigma>2$, it gives a slow shock.

A.3 MHD simulations

We performed 1D 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., 2014). The domain was filled initially with gas of uniform comoving frame density, $\rho=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 $\gamma_{0}$. The gas magnetization is described by the $\sigma$ parameter and tilt angle $\theta$ (zero tilt angle corresponds to field perfectly parallel to the flow). We ran a few simulations for the flow-aligned case ($\theta=0$) as well as for a $5^{\circ}$ tilt. The results are shown in Fig. 5. For the flow-aligned field, the exact decoupling of the fluid motion from the field (Majorana & Anile, 1987) results in a shock independent from $\sigma$. 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, 2005) for intermediate $\sigma$ parameters.

Appendix B PIC simulations

B.1 Simulation Setup

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 $-\mathbf{e}_{x}$) and the reflected beam triggers the formation of a shock, which moves away from the wall along $+\mathbf{e}_{x}$. This setup is equivalent to the head-on collision of two identical plasma shells (Fig. 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 $xy$ plane, with periodic boundary conditions in the $y$ direction. The incoming plasma is injected through a “moving injector,” which recedes from the wall along $+\mathbf{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 2D, each computational cell is initialized with 16 electrons and 16 positrons. The relativistic electron skin depth for the incoming plasma ($c/\omega_{p}$, with $\omega_{p}^{2}=4\pi n_{0}q^{2}/\gamma_{0}m$) is resolved with 10 computational cells and the simulation time-step is $\Delta t=0.045~{}\omega_{p}^{-1}$. The computational domain is typically $\sim 102c/\omega_{p}$ wide (corresponding to 1024 cells). The simulations extend typically up to $\omega_{p}t\sim 3600$, corresponding to a box length of $\sim 3600c/\omega_{p}$, or $\sim 36000$ cells.

The incoming stream is injected along $-\mathbf{e}_{x}$ with bulk Lorentz factor $\gamma_{0}$. The incoming plasma is cold, with thermal spread $k_{B}T/mc^{2}\ll 1$. We vary the upstream Lorentz factor between $\gamma_{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 $\mathbf{B}_{0}$.

Figure 6 shows the result of a series of simulations similar to those displayed on Fig. 3, but for $\gamma_{0}=30$. The weak $\gamma_{0}$-dependence, expected from the micro-physics analysis of Sec. 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 $\theta=5^{\circ}$ with the flow direction of propagation666In 2D simulations, we initialize the magnetic field out of the plane of the simulations. Yet, we have verified that our results do not change for an in-plane field with the same obliquity $\theta$. In this case, we also initialize a motional electric field \mathbf{E}_{0}=-\mbox{\boldmath{\beta}}_{0}\times\mathbf{B}_{0} in the upstream medium, where \mbox{\boldmath{\beta}}_{0}=-\beta_{0}\,\mathbf{e}_{x} is the three-velocity of the injected plasma..

In Fig. 7, we present the results of a suite of 2D simulations of nearly-parallel shocks, with obliquity $\theta=5^{\circ}$. In addition to the $y$-averaged profile of the particle number density (top panel) and the measure of particle anisotropy (bottom panel), 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 $\sigma\gtrsim 0.6$ the downstream flow is no longer isotropic. As a result, the density jump is lower than the MHD predictions.

Appendix C Value of the parameter $\kappa$ in Eq. (8)

The criterion (8) has been obtained in Bret (2016b) assuming equipartition at saturation of the Weibel instability. We here elaborate on the parameter $\kappa$ beyond what is done in Bret (2016b).

For the present system, the growth-rate $\delta_{W}$ of the Weibel instability is given by Eq. (3). The instability grows the field given by Eq. (5) until (Davidson et al., 1972),

[TABLE]

where $\eta=1$ means equipartition, as evidenced by Eqs. (22,23) below777The $\sqrt{2}$ factor in Eq. (19) is not accounted for in Bret (2016b). The magnetic filaments become unable to trap the particles, regardless of their initial $y$-position or velocities, for

[TABLE]

We now replace $B_{f}$ in the equation above by its expression from Eq. (20). We then express $\delta_{W}$ from Eq. (3) and find that $B_{0}>B_{f}/2$ is equivalent to,

[TABLE]

The vertical asymptote defined by the blue curve on Fig. 2 relies therefore on the parameter $\eta$. 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 Eq. (19) to express $B_{f}^{2}$ we find,

[TABLE]

Accounting only for the energy contained in the Weibel field, one finds,

[TABLE]

In the relativistic regime where $\beta_{0}\sim 1$, Eq. (22) gives unity for $\eta=1$, and is monotonically increasing for $\eta\in[0,1]$ and $\sigma<2$. Note that the regime $\sigma>2\beta_{0}^{2}$ is irrelevant in the present context since Eq. (3) implies the instability vanishes there. Eq. (23) implies that the relative amount of energy contained in the Weibel field varies like $\eta^{2}$. For $\eta=1$, Eq. (21) give $\kappa=2/3$, which is the values considered on Fig. 1.