Can We Trust MHD Jump Conditions for Collisionless Shocks?

Antoine Bret

doi:10.3847/1538-4357/aba68e

1. Introduction

Since their discovery during the 19th century (Johnson & Cheret 1998; Salas 2007), shockwaves have been the object of innumerable investigations. The fluid equations first used to describe them operate under the assumption that the mean free path of the particles is much smaller than any other dimension of the system under scrutiny. With such a prominent role given to binary collision to randomize the flow at the microscopic level, it is reasonable to assume (1) that the pressure is isotropic in both the upstream and the downstream and (2) that all the matter upstream goes downstream, together with the energy and the momentum it carries.³ The second assumption allows us to apply the conservation laws between the upstream and the downstream, while the first assumption allows us to write these laws using fluid mechanics or magnetohydrodynamics (MHD) equations. From there, one derives the jump conditions for the density, pressure, magnetic field, etc. (Fitzpatrick 2014; Goedbloed et al. 2019).

Contrary to fluid shockwaves, where dissipation at the shock front is provided by binary collisions, collisionless shockwaves are mediated by collective plasma effects on length scales much shorter than the mean free path (Sagdeev 1966; Tidman & Krall 1971; Balogh & Treumann 2013). A good example is the Earth's bow-shock in the solar wind, where the shock front is about 100 km thick, while the mean free path at the same location is of the order of the Sun–Earth distance (Bale et al. 2003; Schwartz et al. 2011).

In the absence of binary collisions to isotropize the flow, to what extent can we assume isotropic pressures? Also, given the mean free path is much larger than the shock front, to what extent can we assume all the matter upstream goes downstream, together with the momentum and energy it carries? Indeed, it turns out that these two assumptions are far from obvious in a collisionless environment. As a consequence, it also is not obvious whether the fluid or MHD jump conditions derived for a collisional fluid are still valid.

Note that we hereafter refer to MHD jump conditions derived considering isotropic pressures. Several authors adapted them to the case of anisotropic pressures, considering the anisotropy degree as a free parameter (Karimabadi et al. 1995; Erkaev et al. 2000; Vogl et al. 2001; Gerbig & Schlickeiser 2011). Yet the goal of the present paper is to compare jump conditions (mainly the density jump) of collisionless shocks with the simple and well known "isotropic MHD" jump conditions, also frequently referred to as "Rankine–Hugoniot" (RH) jump conditions, even though William Rankine and Pierre Hugoniot derived these relations for a neutral fluid. In the sequel, we shall use "isotropic MHD," "MHD," or "RH" interchangeably.

Two processes have been identified that can trigger a non-RH density jump,

1. An external magnetic field B₀ can sustain stable anisotropies, breaking the isotropy assumption of MHD. Its strength is characterized by the σ parameter,

$\begin{eqnarray}&&\sigma =\displaystyle \frac{{B}_{0}^{2}/4\pi }{({\gamma }_{1}-1){n}_{1}({\displaystyle \sum }_{i}{m}_{i}){c}^{2}},\end{eqnarray} \tag{ 1 }$

where n₁ and ${\gamma }_{1}$ are, respectively, the upstream density and Lorentz factor (measured in the downstream frame). The m_i's are the masses of the species composing the plasma.

2. As they accelerate particles, collisionless shocks generate a population that goes back and forth around the front, breaking the "everything upstream goes downstream" assumption. As we shall see in Section 3.2, the process can be characterized by the parameter,

$\begin{eqnarray}&&\alpha =\displaystyle \frac{{F}_{E}}{\tfrac{1}{2}{n}_{1}{v}_{1}^{3}},\end{eqnarray} \tag{ 2 }$

where F_E is the energy fluxes escaping the RH budget and v₁ is the upstream velocity.

Particle-in-cell (PIC) simulations are undoubtedly the tool par excellence to study nonlinear collisionless phenomena like collisionless shocks. Because they operate from first principles at the microscopic level, they are inherently kinetic. We therefore present in Section 2 a literature review of PIC simulations of collisionless shocks, magnetized or not, in pair or electron/ion plasmas. The observations gathered will then feed Section 3 where departures from MHD density jump are modeled.

Most of the simulations found in the literature use the "reflecting wall" technique to produce a shock. There, a semi-infinite plasma is sent against a reflecting wall where it bounces back to interact with itself. The present work focuses on this technique. In this reflecting scheme, the simulations are therefore performed in the downstream frame of the formed shock. By design, such a scheme can only simulate shocks formed by the encounter of two identical plasmas.

Noteworthily, the less represented "injection method" allows us to study shocks produced by the collision of any two kinds of plasmas (different compositions and/or different densities). Shocks arising from the interaction of a jet with a standing plasma can be studied with this scheme. For example Nishikawa et al. (2009) could study the interaction of a diluted relativistic pair jet with a unmagnetized pair plasma. While many "reflecting wall papers" studied shocks in pair plasmas (see Sections 2.1 and 2.2.1), a density ratio different from unity (Nishikawa et al. 2009 has 0.676) is only achievable with the injection method. Still with the injection method, Ardaneh et al. (2016) studied the interaction of an electron jet with an unmagnetized electron/ion plasma, and commented on the differences between the reflected and injected schemes. Magnetized systems have also been explored, with Dieckmann et al. (2019), for example, considering a pair jet colliding with an electron–proton plasma over a guiding magnetic field.

It appears that the injection method truly allows for an extensive exploration of the possible shocks. The reflected wall scheme restricts the dimension of the parameters' phase space, and to date accounts for more studies, which is why we focus on it here. Yet it would be interesting to extend the current analysis to the injection scheme.

Defining now the density ratio between the shock upstream (subscript "1") and downstream (subscript "2") like,

$\begin{eqnarray}&&r=\displaystyle \frac{{n}_{2}}{{n}_{1}},\end{eqnarray} \tag{ 3 }$

we shall model Δ_RH, the relative departure from the RH jump r_RH, defined by,

$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{RH}}(\sigma ,\alpha )\equiv \displaystyle \frac{r-{r}_{\mathrm{RH}}}{{r}_{\mathrm{RH}}}.\end{eqnarray} \tag{ 4 }$

Beyond the elaboration of a full theory of the density jump accounting for the effects listed above, our present goal is mainly to determine when the RH density jump does apply to collisionless shocks.

2. Literature Review

We conducted a literature review of PIC simulations of collisionless shocks. We selected 68 articles, where (1) a shock structure was clearly obtained, with a downstream significantly longer than the overshoot region, if any, right behind the front, and (2) the density jump can be related to its MHD counterpart with reasonable accuracy, whether explicitly or implicitly. The medium where the shock propagates is homogeneous (see Tomita et al. 2019 for an inhomogeneous case). Save a few exceptions, like Sironi & Spitkovsky (2009a), Stockem et al. (2012), Guo et al. (2018), and Plotnikov et al. (2018), the density jump was not explicitly compared to its MHD counterpart, as this was not the main goal of the article. It is then possible that a few percent discrepancy between the two went unnoticed for some articles.

2.1. Unmagnetized Shocks

All the articles examined but Keshet et al. (2009) and Stockem et al. (2012) pertaining to the unmagnetized regime, relativistic or not, display a shock in agreement with the MHD requirements.

For shocks in pair plasmas (Spitkovsky 2005, 2008a; Kato 2007; Chang et al. 2008; Keshet et al. 2009; Sironi & Spitkovsky 2009b; Bret et al. 2013, 2014; Dieckmann et al. 2016; Dieckmann & Bret 2017, 2018; Li et al. 2017; Lemoine et al. 2019a, 2019b, 2019c; Pelletier et al. 2019; Vanthieghem et al. 2020), the longest simulation (Keshet et al. 2009) was ran up to ${\rm{11,925}}{\omega }_{p}^{-1}$ , where ${\omega }_{p}$ is the electronic plasma frequency.

In Stockem et al. (2012) the authors carefully measured the departure from MHD, as the very goal of the paper was to "assess the impact of non-thermally shock-accelerated particles on the MHD jump conditions of relativistic shocks." Pushing the simulation up to $2395{\omega }_{p}^{-1}$ , a +7% departure for the density jump was found.

Although Keshet et al. (2009) did not precisely measure the density jump, Figure 1 of their article shows an MHD jump at $2250{\omega }_{p}^{-1}$ and a slight departure (+3.5%) at ${\rm{11,925}}{\omega }_{p}^{-1}$ due to energy leakage in accelerated particles (see Section 3.2).

${{\rm{\Delta }}}_{\mathrm{RH}}\lt -4 \% $ — **Figure 1.** Values of Δ_RH defined by Equation (4) for the six articles dealing with magnetized shocks in pair plasmas, in terms of θ_n and σ. Red means negative departure ${{\rm{\Delta }}}_{\mathrm{RH}}\lt -4 \%$ . Blue means positive departure. Green means the shock fits the MHD density jump to within 4%. The labels correspond to the following references: [1]—Spitkovsky (2005), [2]—Sironi & Spitkovsky (2009a), [3]—Bret et al. (2017), [4]—Gallant et al. (1992), [5]—Iwamoto et al. (2017), [6]—Plotnikov et al. (2018). Parametric studies are indicated by vertical lines. [4, 5, 6] are all for *θ_n* = 90^°.
Download figure:
Standard image High-resolution image

For shocks in electron/ion plasmas (Kato & Takabe 2008; Spitkovsky 2008b; Martins et al. 2009; Dieckmann et al. 2010b; Fiuza et al. 2012; Niemiec et al. 2012; Stockem et al. 2014a, 2014b; Ruyer et al. 2015, 2017; Stockem Novo et al. 2015; Naseri et al. 2018; Moreno et al. 2020), the longest simulation time was $4111{\omega }_{\mathrm{pi}}^{-1}$ (Niemiec et al. 2012). No significant departure from the MHD jump was detected in any article.

The case of pair plasmas suggests that departure from MHD due to accelerated particles requires running the simulation for several thousands of electronic plasma frequencies to be perceptible. In electron/ion plasmas, this translates into running the simulation for several thousands of ionic plasma frequencies. The longest run examined in this respect was $4111{\omega }_{\mathrm{pi}}^{-1}$ (Niemiec et al. 2012), where ${\omega }_{\mathrm{pi}}$ is the ionic plasma frequency. Yet the density jump is not measured accurately enough.⁴ Pushing simulations beyond this timescale for non-MHD effects to become clear requires Hybrid simulations as discussed in Section 2.3.

An important feature observed is related to accelerated particles. Their effect on the shock is not steady. As specified in Keshet et al. (2009), "simulations do not reach a steady state; rather, an increasing fraction of shock energy is transferred to energetic particles and magnetic fields throughout the simulation time domain." We shall comment further on this point in Section 3.2.

2.2. Magnetized Shocks

In magnetized media, we have one more source of departure from MHD. There, collisionless shocks can still accelerate particles, which will break the "everything upstream goes downstream" MHD assumption. But now in addition, the field can sustain stable pressure anisotropies and prompt departures from isotropic MHD. Considering the shock behavior strongly depends on the external field strength and on its orientation, the physics of magnetized collisionless shocks is extremely rich.

The upstream field strength is measured by the σ parameter defined by Equation (1). The field orientation is measured by the angle θ_n it makes with the shock front normal. In the nonrelativistic regime, this angle is Lorentz invariant in the direction of the shock propagation up to order ${(v/c)}^{2}$ , where v is the speed of the frame to which the field is transformed. In the relativistic regime, a perpendicular or a parallel field remained so in any frame. The only articles mentioned here where an oblique field is considered in a relativistic setting are Sironi & Spitkovsky (2009a, 2011). There, the angle of the upstream field with the shock normal is given in the simulation frame, that is, the downstream frame.

2.2.1. Magnetized Shocks in Pair Plasmas

Figure 1 summarizes the results for the six articles falling into the present category in terms of θ_n and σ (see references in the caption of Figure 1).

In [2] (Sironi & Spitkovsky 2009a) the departure is about −3% for ${\theta }_{n}\leqslant 30^\circ$ and goes down to −13% for 45°.

In [3] (Bret et al. 2017), the field is parallel and introduces a downstream anisotropy responsible for the departure from the MHD jump. For σ = 3, the jump was reduced by −35%, an effect all the more interesting than for a parallel shock the MHD jump does not depend on the field.⁵

The departure in [1] (Spitkovsky 2005) is directly related to the perpendicular field, and the density jump is said to saturate at two instead of four for larger σ's.

The perpendicular shocks of References [4, 5, 6] are intriguing. The σ-ranges of MHD departure of [4, 5] and [6] do not overlap. [4] (Gallant et al. 1992) finds an increase in the density jump reaching a maximum of +27% for σ = 0.1. [5] (Iwamoto et al. 2017) finds a decrease of the density jump reaching −9% for σ = 0.3. Finally, [6] (Plotnikov et al. 2018) also finds a decrease reaching −4% for $\sigma ={2.10}^{-3}$ .

Table 1 summarizes the main features of these works. We shall comment further on these results in Section 3.1.

Table 1. Simulation Results for References [1, 4, 5, 6] of Figure 1

References	[1]	[4]	[5]	[6]
Δ_RH,max	−50%	+27%	−9%	−4%
σ_m	>0.1	0.1	0.3	2.10⁻³
γ₁	15	40 & 10⁶	40	10

Note. Δ_RH,max is the maximum relative deviation from the RH density jump (see Equation (4)). This deviation is reached for σ = σ_m. γ₁ is the upstream Lorentz factor measured in the downstream frame.

Download table as: ASCII Typeset image

2.2.2. Magnetized Shocks in Electron/Ion Plasmas

Here 22 articles were analyzed, from parallel to normal orientations and σ's ranging from $6\times {10}^{-5}$ (Kato & Takabe 2010) to 0.25 (Dieckmann et al. 2010a).

Figure 2 illustrates the results in the $({\theta }_{n},\sigma )$ phase space. Besides some PIC simulations performed in Guo et al. (2018), all the simulations fulfilled the RH density jump. As specified earlier, a comparison with RH was not the point of some works so that a discrepancy of a few percent may have escaped the analysis.

Guo et al. (2018) performed a detailed comparison with the RH jump for 16 simulations.⁶ Discrepancies in RH range from −0.6% (run "Ms5beta8") to −7% (run "Ms3beta8"). Only discrepancies <−5% have been highlighted in red in Figure 2. The dispersion observed for some identical values of σ stems from different values of the upstream parameters ${\beta }_{p0}=16\pi {n}_{0}{k}_{B}{T}_{0}/{B}_{0}^{2}$ (see Section 3.1).

In summary all of the RH-departures in the examined articles come from the field and decrease the density jump.⁷ An increase stemming from accelerated particles seems to require a few ${10}^{3}{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ (see Section 2.3) to be observed, while the longest simulation in the present section was run up to $559{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ (Fang et al. 2019), where ${{\rm{\Omega }}}_{\mathrm{ci}}$ is the ionic cyclotron frequency.

2.3. Hybrid Results

Hybrid codes treat part of the medium as a fluid, and the rest through the PIC method. In some, the fluid part is the plasma while the PIC method is devoted to accelerated particles (Bai et al. 2015; Casse et al. 2018; van Marle et al. 2018). In others, the electrons are the fluid part, while the ions are dealt with using PIC (Sugiyama 2011; Gargaté & Spitkovsky 2012; Guo & Giacalone 2013; Caprioli & Spitkovsky 2014a; Caprioli & Haggerty 2019; Haggerty & Caprioli 2019).

The advantage of the method is clearly that it allows us to run the simulations longer at a similar computational cost. Such a feature is necessary to render the back-reaction of accelerated particles on the shock itself. Indeed, among the articles examined, the only ones that ran the simulations longer than 10³ ion cyclotron periods ${{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ were Hybrid, with the two longest simulations, Bai et al. (2015), and Haggerty & Caprioli (2019), pushing the computation up to ${\rm{10,800}}{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ and $6000{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ respectively.

Haggerty & Caprioli (2019) and Caprioli & Spitkovsky (2014a) indicate a downstream Maxwellian reaching only 80% of the expected MHD temperature after 6000 and $2500{{\rm{\Omega }}}_{\mathrm{ci}}^{-1}$ respectively, due to "energy leakage" into accelerated particles. Regarding the density jump, Figure 3 shows its increase in terms of the simulated time lengths for various works. A significant difference (from +0% to +75%) is noticeable between Haggerty & Caprioli (2019) ⁸ and Caprioli & Haggerty (2019), due to the way the fluid electrons are modeled. Such is a challenge of Hybrid simulations: giving up a first principle (PIC) description of the electrons is the price to pay for longer simulation times. Much depends then of the fluid closure implemented, as evidenced by these two references.

**Figure 3.** Values of Δ_RH defined by Equation (4) in terms of the simulated time lengths for various Hybrid references. The labels correspond to the following references: [1]—Guo & Giacalone (2013), [2]—Caprioli & Spitkovsky (2014a), [3]—Gargaté & Spitkovsky (2012), [4]—Casse et al. (2018), [5]—Caprioli & Haggerty (2019), [6]—Sugiyama (2011), [7]—Bai et al. (2015).
Download figure:
Standard image High-resolution image

3. Modeling of the Density Jump

Deviations from the RH density jump observed so far are small (except maybe Caprioli & Haggerty 2019, see Section 2.3). We can therefore devise a first-order modeling of ${{\rm{\Delta }}}_{\mathrm{RH}}(\sigma ,\alpha )$ , the departure from RH given by Equation (4), writing

$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{RH}}\left(\sigma ,\alpha \right)\sim \alpha \displaystyle \frac{\partial {{\rm{\Delta }}}_{\mathrm{RH}}}{\partial \alpha }+\sigma \displaystyle \frac{\partial {{\rm{\Delta }}}_{\mathrm{RH}}}{\partial \sigma },\end{eqnarray} \tag{ 5 }$

where all derivatives are considered in (σ,α) = (0, 0).

The α parameter determines the departure from RH due to accelerated particles. From Equation (2), we see it is proportional to F_E, the escaping energy flux. In turn, it is known that the ability of shocks to accelerate particles depends on their magnetization, the angle of the field with the shock normal or the sonic Mach number ${{ \mathcal M }}_{1}$ . In addition, F_E is also an increasing function of time. Therefore, strictly speaking, we should write $\alpha =\alpha (\sigma ,{\theta }_{n},{{ \mathcal M }}_{1},t)$ . A theory of the density jump accounting for all these parameters is out of the scope of this work. Yet, among them the time variable is prominent since,

$\begin{eqnarray}&&\alpha (\sigma ,{\theta }_{n},{{ \mathcal M }}_{1},t=0)=0,\,\,\forall (\sigma ,{\theta }_{n},{{ \mathcal M }}_{1}).\end{eqnarray} \tag{ 6 }$

We shall therefore focus on the time dependence of α, thus deriving its order of magnitude instead of its precise value. We will elaborate further on the effects of accelerated particles in Section 3.2.

3.1. Field Effect on the Density Jump

To what extent can we assume isotropic distribution functions in a collisionless plasma? Here it seems relevant to single out the magnetized and unmagnetized cases.

In a unmagnetized plasma, an anisotropic distribution function is Weibel unstable (Weibel 1959; Kalman et al. 1968). Although Weibel's result was only obtained for Maxwellian distribution functions with anisotropic temperatures, it seems reasonable to conjecture that any anisotropic distribution function is unstable (see Silva et al. 2019 and T. Silva 2020, in preparation, for efforts to obtain mathematical proof). We could also refer to the ample literature on collisionless "Maxwellianization" (see Bret 2015 and references therein), starting with the "Langmuir paradox" (Langmuir 1925).

Indeed, observations of the solar wind show that in the small field limit (high ${\beta }_{\parallel }$ ), the protons' temperature becomes isotropic (Bale et al. 2009; Maruca et al. 2011; Schlickeiser et al. 2011). It seems therefore that past the front turbulence, the downstream anisotropy of a collisionless shock should relax to isotropy on a timescale related to the instability growth rate.

The magnetized case is different because a magnetized Vlasov plasma can sustain stable anisotropies (Gary 1993). A strong enough field B₀ can therefore maintain an anisotropic upstream and/or an anisotropic downstream. An isotropic upstream can turn anisotropic as it goes downstream depending on the magnetization parameter σ.

How should a residual downstream anisotropy modify the density jump? A hint can be given by the fact that the field tends to reduce the degrees of freedom D of the plasma. It therefore increases its macroscopic adiabatic index ${\rm{\Gamma }}=1+2/D$ . And since the RH density jump is a decreasing function of Γ, the presence of the field should lower it. This has been seen in the literature review.

A quantitative assessment of this reduction implies the need to determine the downstream anisotropy in terms of the upstream properties. As noted earlier, the effect is especially interesting for the parallel case because the MHD jump of a parallel shock is σ-independent. Making an ansatz on the kinetic evolution of the plasma through the front, Bret & Narayan (2018) derived for a parallel shock in a nonrelativistic pair plasma (strong shock limit, upstream Γ = 5/3),⁹

$\begin{eqnarray}\begin{array}{rcl}r & = & \displaystyle \frac{1}{2}\left(-\sigma +\sqrt{(\sigma -9)(\sigma -1)}+5\right)\\ & = & 4-\displaystyle \frac{4}{3}\sigma +{ \mathcal O }({\sigma }^{2}).\end{array}\end{eqnarray} \tag{ 7 }$

After some algebra applied to Equation (5), we get

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Delta }}}_{\mathrm{RH}}}{\partial \sigma }=-\displaystyle \frac{1}{3}.\end{eqnarray} \tag{ 8 }$

The perpendicular case is quite different as the MHD jump is already reduced by the field like,¹⁰

$\begin{eqnarray}&&r=4-18\sigma +{ \mathcal O }({\sigma }^{2}).\end{eqnarray} \tag{ 9 }$

Applying the same method as for the parallel case, Bret & Narayan (2019) derived for the perpendicular one (see the Appendix),

$\begin{eqnarray}&&r=4-\displaystyle \frac{86}{3}\sigma +{ \mathcal O }({\sigma }^{2}),\end{eqnarray} \tag{ 10 }$

where $86/3\sim 28.6$ represents a steeper decline than that in the MHD (9). Here, we can write in Equation (5) after some algebra,

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Delta }}}_{\mathrm{RH}}}{\partial \sigma }=-\displaystyle \frac{8}{3}.\end{eqnarray} \tag{ 11 }$

At the present stage it is premature to accurately contrast the model with the simulations, since a full blown theory should account for the composition of the plasmas (pair or e/i) and relativistic effects. Yet the orders of magnitude can be checked at least with Guo et al. (2018) and Bret et al. (2017). The latter is relativistic while the former is not.

The values¹¹ of Δ_RH obtained in Guo et al. (2018) for a perpendicular shock can be fitted by ${{\rm{\Delta }}}_{\mathrm{RH}}\sim -1.14\sigma$ .

The values¹² of Δ_RH obtained in Bret et al. (2017) for a parallel shock can be fitted by ${{\rm{\Delta }}}_{\mathrm{RH}}\sim -0.13\sigma$ . The orders of magnitude fit Equation (8) well for the parallel case and Equation (11) for the perpendicular one.

These results are displayed in Figure 4. The dispersion around the fit for Guo et al. (2018; red) is important due to the various values of β_p0 explored. The dispersion around the fit for Bret et al. (2017) is far less important because all the runs dealt with plasmas that were initially cold.

${{\rm{\Delta }}}_{\mathrm{RH}}=a\sigma $ — **Figure 4.** Values of −Δ_RH defined by Equation (4) in terms of σ, for Guo et al. (2018; red) and Bret et al. (2017; blue). Dashed lines stem from Equations (8) and (11). Plain lines are best fits of the form ${{\rm{\Delta }}}_{\mathrm{RH}}=a\sigma$ , with $a\in {\mathbb{R}}$ .
Download figure:
Standard image High-resolution image

We may finally comment on the "departure windows" observed in Figure 1 for pair perpendicular shocks, and on the +27% increase in the jump observed in Gallant et al. (1992; reference [4] in Figure 1). This increase was attributed to the emission of electromagnetic waves at the shock front. As commented in Plotnikov et al. (2018) (Reference [6] of the present Figure 1), this should be a dimension effect as Gallant et al. (1992) includes the only 1D simulation out of the four references.

The "departure windows"¹³ in [5, 6] (Iwamoto et al. 2017; Plotnikov et al. 2018) were left unexplained. Besides their weak amplitude (−9% and −4%, respectively,¹⁴ for the maximum departure) they may simply arise from the following process: at small σ the jump is in agreement with RH, as evidenced in all the simulations and discussed in Section 3.1. Then the field generates an anisotropy that triggers a negative departure from the isotropic MHD jump. Yet, for even higher σ's, the growing anisotropy drives the jump to 2, and MHD does the same (see Figure 6 in Gallant et al. 1992 and Figure 2 in Plotnikov et al. 2018). We may therefore expect departure windows for perpendicular shocks since MHD and PIC simulations have the same limits in both the weak and the strong field limits.

This is confirmed in Figure 5 where we again plotted Figure 4 for Guo et al. (2018), simply joining the points that share the same β_p0. An MHD-departure window is evidenced for each curve, centered around values of σ similar to those observed in Figure 1 for pair plasmas.

Such a feature is absent in parallel shocks since MHD is insensitive to the field. In parallel geometry, we can expect a departure threshold beyond which kinetic effects progressively drive the collisionless density jump away from its MHD counterpart.

3.2. Accelerated Particle Effect on the Density Jump

In a fluid shock, particles constantly share their energy with each other through binary collisions. In collisionless shocks, it has been known for a long time that particles can be accelerated (Axford et al. 1977; Krymskii 1977; Bell 1978a, 1978b; Blandford & Ostriker 1978). By going back and forth around the front and/or escaping upstream, these particles escape the RH budget. Bret & Pe'er (2018a, 2018b) derived some requirements for particle accelerations that are all fulfilled in the present cases.

A simple calculation derived from Berezhko & Ellison (1999) allows us to conclude that accelerated particles should increase the density jump. We outline it here for completeness. Considering the nonrelativistic regime for simplicity, we start writing the conservation equations between the upstream and the downstream with subscripts "1" and "2," respectively. Labeling ${n}_{i},{v}_{i},{P}_{i},{\rm{\Gamma }}$ the density, velocity pressure, and adiabatic index of the fluid we have,

$\begin{eqnarray}&&{n}_{2}{v}_{2}={n}_{1}{v}_{1}-{F}_{m},\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&{n}_{2}{v}_{2}^{2}+{P}_{2}={n}_{1}{v}_{1}^{2}+{P}_{1}-{F}_{p},\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&\displaystyle \frac{1}{2}{n}_{2}{v}_{2}^{3}+\displaystyle \frac{{\rm{\Gamma }}}{{\rm{\Gamma }}-1}{P}_{2}{v}_{2}=\displaystyle \frac{1}{2}{n}_{1}{v}_{1}^{3}+\displaystyle \frac{{\rm{\Gamma }}}{{\rm{\Gamma }}-1}{P}_{1}{v}_{1}-{F}_{E},\end{eqnarray} \tag{ 14 }$

where ${F}_{m},{F}_{p},{F}_{E}$ are the mass, momentum, and energy fluxes escaping the RH budget because of accelerated particles. It turns out that ${n}_{1}{v}_{1}\gg {F}_{m}$ and ${n}_{1}{v}_{1}^{2}\gg {F}_{p}$ , while F_E in Equation (14) is not negligible with respect to ${n}_{1}{v}_{1}^{3}$ (see Berezhko & Ellison 1999 and references therein). We can therefore neglect ${F}_{m},{F}_{p}$ in Equations (12) and (13) and solve the system for n₂. From Equation (12) we derive ${v}_{2}\,=({n}_{1}/{n}_{2}){v}_{1}$ . Using this expression to eliminate v₂ from Equations (13) and (14) allows us to derive two different expressions for P₂. Equaling them yields a second degree polynomial in n₂ that can be solved exactly. The shock solution reads,

$\begin{eqnarray}&&r=\displaystyle \frac{1+{\rm{\Gamma }}{{ \mathcal M }}_{1}^{2}+\sqrt{{{ \mathcal M }}_{1}^{4}\left(\alpha \left({{\rm{\Gamma }}}^{2}-1\right)+1\right)-2{{ \mathcal M }}_{1}^{2}+1}}{(1-\alpha )({\rm{\Gamma }}-1){{ \mathcal M }}_{1}^{2}+2},\end{eqnarray} \tag{ 15 }$

where α is defined by Equation (2) and,

$\begin{eqnarray}&&{{ \mathcal M }}_{1}^{2}=\displaystyle \frac{{n}_{1}{v}_{1}^{2}}{{\rm{\Gamma }}{P}_{1}},\end{eqnarray} \tag{ 16 }$

is the upstream sonic Mach number. Clearly the density jump (15) can be arbitrarily high as $\alpha \to 1$ . Such a feature can be elaborated further by modeling α (Berezhko & Ellison 1999; Vink et al. 2010).

In the strong shock limit ${{ \mathcal M }}_{1}\to \infty$ it reduces to,

$\begin{eqnarray}\begin{array}{rcl}{r}_{\infty } & = & \displaystyle \frac{\sqrt{\alpha ({{\rm{\Gamma }}}^{2}-1)+1}+{\rm{\Gamma }}}{(1-\alpha )({\rm{\Gamma }}-1)},\\ & = & \displaystyle \frac{{\rm{\Gamma }}+1}{{\rm{\Gamma }}-1}+\alpha \displaystyle \frac{{{\rm{\Gamma }}}^{2}+2{\rm{\Gamma }}+1}{2({\rm{\Gamma }}-1)}+{ \mathcal O }({\alpha }^{2})\\ & = & 4+\displaystyle \frac{16}{3}\alpha +{ \mathcal O }({\alpha }^{2})\,\,\mathrm{for}\,\,{\rm{\Gamma }}=5/3.\end{array}\end{eqnarray} \tag{ 17 }$

After some algebra applied to Equation (5), we find

$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Delta }}}_{\mathrm{RH}}}{\partial \alpha }=+\displaystyle \frac{4}{3},\end{eqnarray} \tag{ 18 }$

so that the relative deviation stemming from accelerated particles should read $+\tfrac{4}{3}\alpha$ . Now, α is not constant in time because the energy F_E poured into cosmic rays is not. For example, the maximum energy of accelerated particles grows like ${t}^{1/2}$ for relativistic shocks (Sironi et al. 2013; Plotnikov et al. 2018) and t for nonrelativistic ones (Caprioli & Spitkovsky 2014b).

These results allow us to phenomenologically assess the α coefficient in (17). Density jump departures stemming from accelerated particles were noted in Keshet et al. (2009) and Stockem et al. (2012; see Section 2.1).¹⁵ For magnetized shocks, such departures were noted in the references featured in Figure 3, among which Caprioli & Haggerty (2019) stands out as the only study in which the density jump is clearly evaluated at various times. These results suggest altogether that the departure reaches a few percent for run times of the order of $5\times {10}^{3}$ time units. For magnetized shocks, this "time unit" is the gyroperiod. For the unmagnetized shocks in pair plasmas, it is the inverse electronic plasma frequency.

The orders of magnitude derived from unmagnetized shocks in pairs (Keshet et al. 2009; Stockem et al. 2012) and magnetized parallel shocks in electron/ion (Caprioli & Haggerty 2019), are similar. To what extent can they be generalized to any shock? Indeed, these three studies are but a sample of the possible combinations achievable varying the parameters $(\sigma ,{\theta }_{n},{{ \mathcal M }}_{1})$ in pairs and electron/ion. It turns out that the acceleration efficiency has been studied extensively in the magnetized pair (Sironi & Spitkovsky 2009a) or electron/ion plasmas (Sironi & Spitkovsky 2011; Caprioli & Spitkovsky 2014a; Guo et al. 2014a, 2018), and it was found that the aforementioned order of magnitude is representative of the full range of possibilities (L. Sironi 2020, private communication).

Assuming that the departure grows like $\beta {t}^{\kappa }$ with $\kappa \gt 0$ , we then set,

$\begin{eqnarray}&&\displaystyle \frac{4}{3}\beta {t}^{\kappa }={ \mathcal O }({10}^{-1})\,\,\mathrm{for}\,\,t={ \mathcal O }(5\times {10}^{3}),\end{eqnarray} \tag{ 19 }$

so that,

$\begin{eqnarray}\begin{array}{rcl}\alpha & = & { \mathcal O }({5}^{-\kappa }{10}^{-1-3\kappa }){t}^{\kappa },\\ & = & { \mathcal O }({10}^{-1-3.7\kappa }){t}^{\kappa },\end{array}\end{eqnarray} \tag{ 20 }$

where t is measured in the dominant unit of the simulation. Note that the scaling of the maximum energy of the accelerated particles does not have to translate to the energy flux F_E (see Equation (2)) leaking into accelerated particles, since the most energetic ones are but a few. In all the papers examined, the only one from which it was possible to extract a time dependent variation of F_E is Caprioli & Haggerty (2019). Their Figure 2 suggests $\kappa ={ \mathcal O }(1)$ . Further works would be welcome to narrow down the value of κ and the timescale of $5\times {10}^{3}$ set by Equation (19).

3.3. Summary

Gathering the results obtained in Section 3.1 for the field effects, and in Section 3.2 for accelerated particles, we can complete Equation (5) and write,

$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{RH}}(\sigma ,\alpha )\sim +{ \mathcal O }({10}^{-1-3.7\kappa })t-\sigma \left\{\begin{array}{l}1/3\,\,{\theta }_{n}=0,\\ 8/3\,\,{\theta }_{n}=\pi /2.\end{array}\right.\end{eqnarray} \tag{ 21 }$

Since we are interested in the order of magnitude of the field correction, we can aggregate the results for ${\theta }_{n}=0$ and π/2 and propose,

$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{RH}}(\sigma ,\alpha )\sim +{ \mathcal O }({10}^{-1-3.7\kappa }){t}^{\kappa }-\sigma { \mathcal O }(1),\end{eqnarray} \tag{ 22 }$

where κ is of order unity. As commented above, the first term, $+{ \mathcal O }({10}^{-1-3.7\kappa }){t}^{\kappa }$ , should be representative, in order of magnitude, of the full spectrum of shocks populating the $(\sigma ,{\theta }_{n},{{ \mathcal M }}_{1})$ phase space.

To which extent can we make the same claim for the correction term ∝σ, due to the field driven anisotropy?

The field correction obviously vanishes for σ = 0, $\forall ({\theta }_{n},{{ \mathcal M }}_{1})$ . Regarding the θ_n variation, Figures 1 and 2 present the results of nine simulations at various angles in pair plasmas, and 31 in electron/ion plasmas. No difference of order of magnitude has been detected with respect to Equation (22) for the coefficient of σ.

As for the incidence of ${{ \mathcal M }}_{1}$ , all the pair shocks featuring Figure 1 are strong, that is, ${{ \mathcal M }}_{1}\gg 1$ . As for those in electron/ion presented in Figure 2, they span sonic Mach numbers ranging from 2 (Guo et al. 2018; Ha et al. 2018; Kang et al. 2019) to ${{ \mathcal M }}_{1}\gg 1$ . Again no deviation, order of magnitude wise, has been detected from the coefficient of the σ correction reported in Equation (22).

The only part of the phase space parameter that has not been tested in the literature is the deviation from RH in weak shocks in pair plasmas. If the model developed in Bret & Narayan (2018, 2019, 2020) is further confirmed by PIC simulation, then this gap of weak shock in pairs will be filled.

4. Summary and Conclusion

The present work represents an attempt to determine when the RH density jump can be applied to collisionless shocks. A tentative answer is given by Equation (22), which is valid as long as Δ_RH, the relative departure from the RH density jump, is small.

The departure is the sum of two terms. One, positive, arises from the accelerated particles escaping the RH budget, and grows with time. The second is negative and stems from the pressure anisotropy sustained by the field.

The literature review clearly evidences positive departures arising from particle acceleration, and negative ones from field driven anisotropies. We did not find studies considering both effects together, contemplating, for example, the possibility that they compensate for each other.

What about the long times and/or strong field regime, beyond the validity of Equation (22)?

Can accelerated particles drive a density jump that is arbitrarily high on the long run? Some theoretical models suggest so (Berezhko & Ellison 1999; Vink et al. 2010) together with some simulations (Caprioli & Haggerty 2019). Yet, since the maximum energy of these particles saturates with time (Sironi et al. 2013), the energy flux F_E escaping the RH budget in Equation (14) could also saturate.

As for the large σ regime, we have to distinguish parallel shocks from perpendicular ones.

For parallel shocks, the MHD jump is insensitive to the field. The jump departure is therefore going to increase until it saturates since a large σ will eventually trigger an anisotropy yielding an asymptotic density jump. For example, for a strong shock in 3D with Γ = 5/3, the MHD jump remains four while the collisionless one tends toward two, resulting in the largest possible departure ${{\rm{\Delta }}}_{\mathrm{RH}}=-50 \%$ .

For perpendicular shocks, the large σ jumps are identical for MHD and collisionless shocks. The negative jump departure should therefore reach a maximum for intermediate values of σ and vanish before and after. Such a "departure window" has been retrieved for relativistic pair shocks (see Iwamoto et al. 2017; Plotnikov et al. 2018 and Figure 1 of the present work) and nonrelativistic electron/ion shocks (Guo et al. 2018). In the later case, Figure 5 shows that the location and magnitude of the window depend on the upstream proton temperature parameter β_p0. At any rate, Δ_RH never goes below −10% in any of the aforementioned studies.

The dimensionality of the works involved in the present study may be its main limitation. Even though the formalism of Bret & Narayan (2018, 2019, 2020) is 3D, simulations under scrutiny feature at best three velocity dimensions but only two spatial dimensions (2D3V). To what extent can the conclusions be generalized to 3D space? The reduced number of spatial dimensions can have important effects on particle acceleration and/or the external field effect. Comparisons between 2D and 3D results found, for example, that some particle trapping occurs in 2D and not in 3D (Cruz et al. 2017; Trotta & Burgess 2018). Also, considering only two spatial dimensions necessary excludes waves and instabilities with a wavevector k oriented along the excluded spatial dimension. Indeed, theoretical explorations of the full unstable k-spectrum of some beam-plasma (or Weibel-like) instabilities, found it is truly 3D (Kalman et al. 1968; Dieckmann et al. 2008; Bret 2014; Novo et al. 2016).

Finally, Matsumoto et al. (2017; injection scheme, Mach number =22.8) explicitly compared electron acceleration in 2D and 3D simulations for quasi-perpendicular electron/ion shocks. They found the acceleration to be more efficient in 3D than in 2D, be it with an out-of-plane or an in-plane field.¹⁶

A very similar system was studied in Guo et al. (2014a; reflecting wall, Mach number =3). There, 2D simulations were also compared to 3D ones. Yet, in contrast with Matsumoto et al. (2017), it was found that the 2D in-plane field configuration "is a good choice to capture the acceleration physics of the full 3D problem." Perhaps the discrepancy with Matsumoto et al. (2017) is due to the different Mach numbers or the different methods. Therefore, it seems that, as concluded in Bohdan et al. (2017), "true 3D simulations are urgently needed to resolve this issue."

This work has been achieved under projects ENE2016-75703-R from the Spanish Ministerio de Economía y Competitividad and SBPLY/17/180501/000264 from the Junta de Comunidades de Castilla-La Mancha. Thanks are due to Colby Haggerty, Thales Silva, Arno Vanthieghem, Laurent Gremillet, Damiano Caprioli, and Lorenzo Sironi for valuable input.

Appendix: Proof of Equation (10)

We start from Equation (25) of Bret & Narayan (2019) where the parameter χ₁ is proportional to the sonic Mach number. This equation already has the adiabatic index Γ =5/3. The equation for the density jump r in the strong shock limit ${\chi }_{1}\to \infty$ is given by the coefficient of ${\chi }_{1}^{2}$ of Equation (25) in Bret & Narayan (2019). The result is the third-degree polynomial,

$\begin{eqnarray}&&P(r)\equiv (2r(r+2)-5)r\sigma +2(r-5)r+8=0,\end{eqnarray} \tag{ A1 }$

for which it can be checked that for σ = 0, the two roots are r = 1 and 4. We then set r = 4+kσ with $k\in {\mathbb{R}}$ , and perform a Taylor expansion of $P(4-k\sigma )$ up to first order in σ. The zeroth order vanishes, and the first order does as well if k = −86/3, proving Equation (10).