Density jump for oblique collisionless shocks in pair plasmas: physical solutions

Abstract

Collisionless shocks are frequently analysed using the magnetohydrodynamics (MHD) formalism, even though MHD assumes a small mean free path. Yet, isotropy of pressure, the fruit of binary collisions and assumed in MHD, may not apply in collisionless shocks. This is especially true within a magnetized plasma, where the field can stabilize an anisotropy. In a previous article (Bret & Narayan, J. Plasma Phys., vol. 88, no. 6, 2022b, p. 905880615), a model was presented capable of dealing with the anisotropies that may arise at the front crossing. It was solved for any orientation of the field with respect to the shock front. Yet, for some values of the upstream parameters, several downstream solutions were found. Here, we complete the work started in Bret & Narayan (J. Plasma Phys., vol. 88, no. 6, 2022b, p. 905880615) by showing how to pick the physical solution out of the ones offered by the algebra. This is achieved by 2 means: (i) selecting the solution that has the downstream field obliquity closest to the upstream one. This criterion is exemplified on the parallel case and backed up by particle-in-cell simulations. (ii) Filtering out solutions which do not satisfy a criteria already invoked to trim multiple solutions in MHD: the evolutionarity criterion, that we assume valid in the collisionless case. The end result is a model in which a given upstream configuration results in a unique, or no downstream configuration (as in MHD). The largest departure from MHD is found for the case of a parallel shock.

1. Introduction

Shock waves are fundamental processes in fluids. They have been the subject of numerous studies for nearly two centuries (Salas 2007). When the frequency of collisions between particles is high, the thickness of the shock front is of the order of a few mean free paths, since binary collisions are ultimately the only microscopic mechanism capable of transferring some of the kinetic energy of the upstream medium into heat in the downstream (Zel'dovich & Raizer 2002).

However, in situ observations of the bow shock of the Earth's magnetosphere in the solar wind have shown that its front is approximately a hundred kilometres thick, while the mean free path at the same location is of the order of the Sun–Earth distance (Bale, Mozer & Horbury 2003; Schwartz et al. 2011). Such a shock cannot be mediated by collisions. It is mediated by collective plasma electromagnetic effects (Sagdeev 1966). This type of shock is known as a ‘collisionless shock’.

Strictly speaking, collisionless shocks should be studied at the kinetic level, using the Vlasov equation, since the absence of collisions can even make it difficult to define a local velocity, as is the case, for example, in counter streaming systems. Due to the complexity involved in solving the Vlasov equation, collisionless shocks, and in particular the density, temperature or velocity jumps they present, are often interpreted via magnetohydrodynamics (MHD).

Yet MHD relies on hydrodynamics and as such entails the same hypothesis of small mean free path (Gurnett & Bhattacharjee 2005; Goedbloed, Keppens & Poedts 2010; Thorne & Blandford 2017). This hypothesis of small mean free path has 2 consequences that are important for the calculation of the density jump around a shock. The first consequence is that pressure is isotropic, both before and after the shock. In fact, even if a fluid is subjected to pressure anisotropy during its transport from the upstream to the downstream, binary collisions will restore isotropy on a time scale of the order of the collision frequency, well below the macroscopic times involved. The second consequence is that all of the upstream fluid passes into the downstream, along with the matter and momentum it carries.

It turns out that, in a collisionless shock, these 2 consequences can be invalidated (Bret 2020). The first, because in the absence of collisions, a plasma can maintain a stable anisotropy in the presence of an ambient magnetic field (Hasegawa 1975; Gary 1993). The second consequence, because due to the large mean free path, plasma particles can bounce off the shock front or even travel upstream from the downstream, as is the case with accelerated particles (Drury 1983; Blandford & Eichler 1987).

This article proposes a remedy to the first consequence: How to correct the MHD jump equations so that they can account for an anisotropy in the plasma?

Notably, in the absence of an ambient magnetic field, the Weibel instability ensures isotropy of a collisionless plasma (Weibel 1959; Silva, Afeyan & Silva 2021). Therefore, the MHD jump equations only need to account for anisotropies if a magnetic field is present.

Several authors have derived the MHD jump equations for the non-isotropic case (Hudson 1970; Karimabadi, Krauss-Varban & Omidi 1995; Erkaev, Vogl & Biernat 2000; Gerbig & Schlickeiser 2011). But in all of these works, while the anisotropy of the upstream is considered a free parameter, so is that of the downstream. These equations are therefore unable, on their own, to derive the density jump of a shock whose downstream is not isotropic, because they lack precisely one parameter: the anisotropy of the downstream.

In a recent series of papers, we developed a model that precisely fills this gap. It was first applied to the case of a parallel shock (Bret & Narayan 2018), i.e. a shock moving parallel to the ambient magnetic field. The assumptions made and the results obtained were confirmed by numerical simulations (Haggerty, Bret & Caprioli 2022).

The model was then applied to the case of a perpendicular shock (Bret & Narayan 2019), and finally to the cases of a switch-on and of an oblique shock (Bret & Narayan 2022a,b). The latter case, that of an oblique shock, is analytically much more complicated than the parallel and perpendicular cases, due to the complexity of the MHD jump equations for an oblique field and anisotropic temperatures.

In Bret & Narayan (2022b), hereafter referred to as Paper 1, the algebra of these equations was solved, but the solutions were left unfiltered. As a result, several coexisted for some combinations of the upstream parameters.

Here, Paper 1 is completed by filtering the algebraic solutions so that, for a given combination of upstream variables, there is no more than 1 solution for the downstream.

2. Method

The system under scrutiny is represented in figure 1. It is identical to that of Paper 1. Although the conservation equations we shall use are valid for any upstream temperatures, we here, as in Paper I, only treat the strong sonic case $T_{\perp 1}=T_{\parallel 1}=0$.

Figure 1. System under scrutiny. It is identical to that of Paper 1. Although the conservation equations used are valid for any upstream temperatures, we here, as in Paper 1, only treat the strong sonic case $T_{\perp 1}=T_{\parallel 1}=0$. We work in the frame of reference where $\boldsymbol {v}_1$ is normal to the front.

2.1. Summary of previous works

As previously said, the basic caveat of MHD is that, if a collisionlessly stable anisotropy develops at the front crossing, MHD itself cannot derive it. The jump of quantities like the density is therefore under-determined.

For completeness, we now briefly recall the results obtained in previous works.

In Bret & Narayan (2018) we presented a model capable of solving this issue for a parallel shock. We reasoned that, as it crosses the front, the plasma is compressed in the direction parallel to the motion. As a consequence, its parallel temperature increases while its perpendicular temperature remains constant. The 3 MHD conservation equations (matter, momentum, energy¹ ) are therefore completed by

(2.1)\begin{equation} T_{{\perp} 2} = T_{{\perp} 1}, \end{equation}

allowing us to derive the 4 downstream unknowns ($n_2,v_2,T_{\perp 2},T_{\parallel 2}$), in terms of the upstream variables. Note that, here, the perpendicular direction is common to the flow and the field.

Now, the state of the downstream resulting from the conservation of $T_\perp$ may be stable, or not. If it is stable, then this is the end state of the downstream. If it is unstable, the plasma migrates to the instability threshold.² Imposing marginal stability then provides again a fourth equation allowing us to fully determine the properties of the downstream.

Bret & Narayan (2018), as well as every subsequent works of ours on the same model, is limited to pair plasmas. The reason for this is that the equality of the mass of the species allows us to consider only one parallel and one perpendicular temperature. In an electron/ion plasma where electrons and ions are heated differently in the shock, a 4 temperature model would be required (Guo, Sironi & Narayan 2017, 2018). Yet, since the model eventually relies on macroscopic physics, it should also apply to electron/ion plasmas, as preliminary particle-in-cell (PIC) simulations seem to indicate (Shalaby 2024).

The model predicted, for a strong sonic parallel shock, a density jump going from 4 to 2 in the high field regime, a prediction successfully confirmed by PIC simulations in Haggerty et al. (2022). Such a result stands in contrast with MHD, where the density jump should always be 4, regardless of the field strength.³

The perpendicular case was treated in Bret & Narayan (2019). There, the direction perpendicular to the flow is eventually parallel to the field, so that the counterpart of (2.1) is

(2.2)\begin{equation} T_{{\parallel} 2} = T_{{\parallel} 1}. \end{equation}

The switch-on shocks, where the field is oblique in the downstream only, was treated in Bret & Narayan (2022a). The model has also been solved for a parallel or a perpendicular shock, with an anisotropic upstream (Bret 2023a,b).

Finally, the general case where the field may be oblique in both the upstream and the downstream was treated in Bret & Narayan (2022b), namely, Paper 1.

In Bret & Narayan (2022a,b), the closure of the MHD jump equations was achieved through an ansatz interpolating between (2.1) and (2.2). In the limit of a cold upstream with $T_1=0$, which is the regime treated in Paper 1 and hereafter, the ansatz reads

(2.3)\begin{equation} \left.\begin{array}{c@{}} T_{{\parallel} 2} = T_e\cos^2\theta_2, \\ T_{{\perp} 2} = \dfrac{1}{2}T_e\sin^2\theta_2, \end{array}\right\} \end{equation}

where $T_e$ is a parameter determined when solving the equations and $\theta _2$ is the angle of the downstream field with the shock normal (see figure 1). Equations (2.3) correctly reduce to (2.1) and (2.2) in their respective limits since $\theta _2=0$ for the parallel case, while $\theta _2={\rm \pi} /2$ for the perpendicular one. Such a scheme guarantees that both perpendicular temperatures are equal, as required by the Vlasov equation (Landau & Lifshitz 1981, § 27). Also, the total thermal energy in the 3 directions sums up to $k_BT_e$, where $k_B$ is the Boltzmann constant.

In summary, and since the 2 instabilities involved are the firehose and the mirror instabilities (see § 4), our model can be stated as follow:

1. • As the plasma goes through the shock front, its temperature normal to the flow is conserved for the parallel and perpendicular cases. This translates directly to (2.1) and (2.2), respectively. For the oblique case with cold upstream $T_1=0$, (2.3) interpolate between these 2 extremes.

2. • The resulting state of the plasma in the downstream is called ‘stage 1’.

3. • If the downstream field is strong enough to stabilize stage 1, then this is the end state of the downstream.

4. • If the downstream field is too weak to stabilize stage 1, then

   1. • If stage 1 is firehose unstable, it migrates to the firehose instability threshold. This is ‘stage 2 – firehose’.

   2. • If stage 1 is mirror unstable, it migrates to the mirror instability threshold. This is ‘stage 2 – mirror’.

2.2. Present work

Paper 1 has 3 kinds of limitations:

1. (a) It is restricted to strong sonic shocks, namely $T_1=0$, and to non-relativistic pair plasmas. These restrictions are still considered here.

2. (b) It considers the simplest expressions for the Alfvén velocity and the stability criterion of the instabilities involved in our model. Yet, more accurate expressions are required in an anisotropic plasma. The present work accounts for one of them.

3. (c) It only presents the allowed solutions to the conservation equations, plus (2.3). It does not filter these solutions according to their physical relevance. Such a filtering is the main goal of the present work.

Our purpose here is to deal with points (b) and (c) above.

Besides the variables explained on figure 1, we shall use the following dimensionless parameters:

(2.4)\begin{equation} \left.\begin{array}{c@{}} \displaystyle r = \dfrac{n_2}{n_1}, \\ \displaystyle A_i = \dfrac{T_{{\perp} i}}{T_{{\parallel} i}}, \\ \displaystyle M_{A,i} = \dfrac{v_i}{v_{A,i}}, \\ \displaystyle \sigma = \dfrac{B_1^2/8{\rm \pi}}{\dfrac{1}{2}n_1 m v_1^2} \equiv \dfrac{1}{M_{A,1}^2}, \end{array}\right\} \end{equation}

where $v_{A,i}$ is the Alfvén velocity

(2.5)\begin{equation} v_{A,i} = \frac{B_i}{\sqrt{4 {\rm \pi}n_i m}}. \end{equation}

The parameter $\sigma$ is frequently used in simulations of collisionless shocks like Haggerty et al. (2022), while the Alfvén Mach number is common in the MHD shock literature.

In addition to the Alfvén Mach number defined above, we shall often use in the sequel its following variant:

(2.6)\begin{equation} M_{Aix} \equiv \frac{v_i \cos\xi_i}{v_{A,i}\cos\theta_i}. \end{equation}

It compares the projection of the flow velocity along the shock normal ($x$ axis) with the projection of the Alfvén velocity, still along the shock normal.

Since the road map for solving our model is eventually the one already used in MHD, we start by reminding the reader how it applies there.

3. Isotropic MHD results and evolutionarity

One of the criteria used here to filter out some solutions is the so-called ‘evolutionarity criterion’, already present in isotropic MHD. It is therefore convenient to present how it operates there.

The MHD conservation equations for isotropic temperatures are reported in Appendix A. They can be used to derive an expression of the downstream field angle $\theta _2$ in terms of the upstream quantities only. Namely, the quantity

(3.1)\begin{equation} T_2 \equiv \tan\theta_2, \end{equation}

is a root of the polynomial (A9) in Appendix A. The MHD density jump is then given by

(3.2)\begin{equation} r = \frac{M_{A1}^2 T_2}{M_{A1}^2 \tan \theta_1+T_2 \cos ^2\theta_1-\sin \theta_1 \cos \theta_1}, \end{equation}

in terms of the upstream Alfvén Mach number (2.4).

Although such a derivation of the MHD density jump is uncommon in the literature, it mimics the derivation of the jump in our model. It is therefore pedagogically valuable for the present work.

Figure 2 presents the solutions for 3 different angles $\theta _1$. The upper row shows all the possible solutions. For $\theta _1=0$ and $M_{A1x}\in [1,2]$, there are 2 MHD branches, the lower one pertaining to the switch-on solutions. For $\theta _1=0.175 {\rm \pi}/2$ there are even 3 solutions for $M_{A1x}\in [1, 1.34]$.

Figure 2. Upper row: density jump $r$ for different angles $\theta _1$ in terms of $M_{A1x}=M_{A1}/\cos \theta _1$, defined by (2.6). Sometimes there is more than 1 solution. Second row: angle $\theta _2$ of the field $\textbf {B}_2$ with the $x$ axis (see figure 1). Third row: evolutionarity criterion. Some branches, namely the ones crossing the shaded areas, are to be excluded. Lower row: same as upper row, but without the branches excluded by the evolutionarity criterion. There is now but 1 solution for a given $M_{A1x}$, or none. The black circles for $\theta _1=0.5{\rm \pi} /2$ are the results of our MHD simulations (see § 3.1).

Which one will the shock choose? This question is at the heart of this work. Let us now see how it is solved in MHD.

The MHD answer relies on the notion of shock ‘evolutionarity’, which has been discussed several times in the literature (e.g. Kennel, Blandford & Wu 1990; Farris et al. 1994; Falle & Komissarov 2001; Kulsrud 2005). For given upstream and downstream boundary conditions, the MHD Rankine–Hugoniot jump conditions give a unique solution for fast shocks, where the flow speed on both sides of the shock are super-Alfvénic, and also for slow shocks, where the fluid is throughout sub-Alfvénic. These two types of shocks are stable and are described as evolutionary. Four other potential shock types, each of which has super-Alfvénic upstream fluid and sub-Alfvénic downstream fluid, have no unique solutions to the Rankine–Hugoniot relations. In these shocks, the transverse magnetic field switches sign across the shock, and the fluid equations do not provide the correct number of Alfvén waves to handle the field flip (Falle & Komissarov 1997; Kulsrud 2005). Such shocks will not arise naturally from generic initial conditions. Even if they are artificially set up, they will tend to deviate quickly from their initial configuration, typically splitting into two shocks, one fast and the other slow. These ‘forbidden’ shocks do not satisfy the evolutionarity condition. They are called ‘intermediate’ (Falle & Komissarov 1997) or ‘extraneous’ (Kulsrud 2005) shocks.

The MHD solutions presented on the upper row of figure 2 need therefore to be filtered according to the aforementioned evolutionarity criterion. The upstream Alfvén Mach number to consider for the evolutionarity analysis is not the one given by (2.4), but its variant (2.6) instead. With this definition, the downstream Alfvén Mach number reads

(3.3)\begin{equation} M_{A2x}^2 = 1+\frac{\tan \theta_1 (M_{A1}^2 \sec ^2\theta_1-1)}{T_2}. \end{equation}

It is with definition (2.6) that the Alfvén Mach number of the MHD switch-on shocks, see left plot of third row in figure 2, is found equal to 1 (Goedbloed et al. 2010, p. 853). These modes, by definition, do not have $\theta _2=0$ (second row in figure 2) nor $\xi _2=0$. Corrections (2.6) to the Alfvén Mach number (2.4) are therefore important here.

The value of $M_{A2x}$ is plotted in terms of $M_{A1x}$ in the third row of figure 2. The forbidden, non-evolutionary, zones have been coloured in the plots. They feature the non-evolutionary transition just described, namely super-Alfvénic $\rightarrow$ sub-Alfvénic, together with the reverse transition sub-Alfvénic $\rightarrow$ super-Alfvénic. As a result, the MHD branches going through these regions are non-evolutionary, hence not physical solutions of the MHD problem.

The lowest row of figure 2 is the result of this evolutionary filter applied to the upper row. There is now but 1 solution for a given $M_{A1x}$, or none, as there are $M_{A1x}$-gaps where no solutions appear. Regardless of the initial conditions, the shock formed is never found with a $M_{A1x}$ inside such gaps like, for example, $M_{A1x}\in [1 , 1.4]$ for $\theta _1=0.5{\rm \pi} /2$ (see figure 2, bottom right plot).

We shall assume in the sequel that the evolutionarity criteria also applies in the collisionless case. This will have to be checked in future works.

3.1. Magnetohydrodynamics simulations

In order to illustrate these theoretical results, we ran some MHD simulations with the code KORAL (Sa̧dowski et al. 2013, 2014). KORAL is designed for multi-dimensional radiative MHD simulations in general relativistic space–times. However, if we turn off the radiation module as well as special and general relativity, the code reduces to a multi-dimensional non-relativistic MHD code. This version of the code was used here. Some of the results for $\theta _2$, $M_{A2x}$ and $r$ in terms of $M_{A1x}$ for upstream $\theta _1=0.5{\rm \pi} /2$, are pictured in figure 2 by the black circles. They perfectly line up with the theory and avoid the predicted gap in the solutions for $M_{A1x}\in [1 , 1.4]$.

Figure 3 shows the result of a simulation of a shock in a non-evolutionary case. At $t=0$, two cold fluids of identical density and opposite velocities $\pm 1.5 \boldsymbol {x}$ collide (blue lines). The magnetic field has modulus unity and $\theta _1=0.5{\rm \pi} /2$. At $t=1.5$, two shocks formed instead of one. The non-evolutionary case gives rise to 2 sub-shocks, as predicted in Kulsrud (2005) for example.

Figure 3. Magnetohydrodynamics simulation in which the evolutionarity criterion results in the formation of 2 sub-shocks instead of one. The dashed blue lines show the initial fluid configuration at $t=0$, when two fluid columns collide at $x=0$. The solid red lines show the configuration at $t=1.5$. Each fluid column develops two shocks, a fast shock at $|x|=1.58$ and a slow shock at $|x|=0.34$. The little dip in density $n$ at $x=0$ for $t=1.5$ is a numerical artefact.

4. Instabilities involved and modified Alfvén velocity

As previously stated, the starting point for our model are the MHD conservation equations with anisotropic temperatures. They have been derived in Hudson (1970) and are reported in Appendix B with the notations of Bret & Narayan (2022a,b).

For stage 1, we solve them imposing relations (2.3). This defines stage 1. Then, mirror or firehose stability of stage 1 has to be assessed. These instabilities are discussed here.

Now, evolutionarity involves the downstream Alfvén velocity (projected onto the shock normal). As long as the plasma is isotropic, the Alfvén velocity in given by (2.5). As a result, the downstream Alfvén Mach number for the MHD switch-on shocks is exactly 1.

Yet, in an anisotropic plasma, it was found in Abraham-Shrauner (1967) that the Alfvén velocity reads⁴

(4.1)\begin{equation} c_A ={\pm} v_A \sqrt{S_\perp{-} S_\parallel{+}1} \cos\theta, \end{equation}

where $v_A$ is still given by (2.5) and

(4.2)\begin{equation} S_{{\perp},\parallel} = \frac{n k_B T_{{\perp},\parallel}}{B^2/4{\rm \pi}}. \end{equation}

In (4.1), the $\pm$ sign refers to waves propagating along the field, or in the opposite direction. As shall be checked in the sequel (§§ 5.1 and 5.2 and figure 5), with this correction to the Alfvén velocity, the switch-on solutions of our model have $M_{A2x}=1$, exactly like in MHD.

In addition, the quantity below the square root can become negative. Such a situation indicates that the Alfvén waves become unstable, which corresponds to the firehose instability. For this to happen, $S_\perp - S_\parallel +1 < 0$ is required, which can be cast under the form

(4.3)\begin{equation} A \equiv \frac{T_\perp}{T_\parallel} < 1 - \frac{2}{\beta_\parallel}, \end{equation}

where⁵

(4.4)\begin{equation} \beta_\parallel{=} \frac{n k_B T_\parallel}{B^2/8{\rm \pi}}. \end{equation}

The criterion (4.3) is therefore the one we shall use in the sequel, in order to preserve the inner coherence with the switch-on solutions of our model having $M_{A2x}=1$. This criterion slightly differs, by the factor 2, from the one commonly used in the literature and in Paper 1.

The Solar Wind data are a key test for the threshold of the firehose instability. They are indeed compatible with both criteria, namely with or without the factor 2 (Hellinger et al. 2006; Bale et al. 2009; Maruca, Kasper & Bale 2011).

The other instability considered in relation to the stability of stage 1, is the mirror instability. While the firehose instability occurs for too low an anisotropy $T_\perp /T_\parallel$, the mirror instability occurs for too high an anisotropy. The standard threshold given in the literature reads (Hasegawa 1975; Gary 1993)

(4.5)\begin{equation} A > 1 + \frac{1}{\beta_\parallel}. \end{equation}

Yet, a different criterion is given in Abraham-Shrauner (1967), namely

(4.6)\begin{equation} A > 1 + \frac{1}{A ~\beta_\parallel}. \end{equation}

However, while stage 2-firehose remains analytically tractable using (4.3), stage 2-mirror is not when using (4.6) instead of (4.5). Therefore, we adopt criterion (4.3) for the firehose instability and keep (4.5) for mirror. Note in this respect that (i) the Solar Wind data are compatible with both criteria and (ii), the inner coherence of the model does not impose a specific mirror criteria, as is the case for firehose.

Figure 4 pictures the various criteria commented here for the mirror and the firehose instabilities. Even if the corrected criterion for the mirror instability is functionally different than the one without the correction, the end result is qualitatively the same, and remains compatible with the Solar Wind data.

Figure 4. Stability diagram. The solid lines picture the criterion used here. The dashed lines pertain to other criteria discussed in § 4.

The two stability domains are disconnected, so that the plasma cannot be unstable to both at once. There is no possible competition between them.

5. Branch selection

Having specified the instabilities involved in the transition from stage 1 to stage 2, together with their respective thresholds, we can proceed to the filtering of the solutions our model offers. In this respect, it is instructive to single out the case $\theta _1=0$.

5.1. Case $\theta _1=0$

Figure 5(a) shows the solutions offered by our model for $\theta _1=0$, without any filtering. It displays the various solutions with stable stage 1, plus a branch, in green, corresponding to stage 2-firehose because stage 1 is firehose unstable for $\sigma \in [0 , 0.5]$. For $\sigma \in [0.5,1]$, there are up to 3 possible solutions: one has $r=2$ and $\theta _2 = 0$, and the other two pertain to the switch-on solutions in red, with $\theta _2 \neq 0$.

Figure 5. (a) All solutions offered by our model for $\theta _1=0$. (b) All the solutions fulfil the evolutionarity criterion. Note that this is different from the isotropic MHD problem, where some solutions for $\theta _1=0$ enter the shaded forbidden zone (figure 2, row 3). Note also that, with definition (4.1) of the Alfvén speed, $M_{A2x}=+\infty$ on the firehose threshold. This is why the green branch on the left is sent to $M_{A2x}=+\infty$ on the centre plot. Notably, $M_{A2x}=1$ for the switch-on shocks of our model (red segment on centre plot). For this analysis, the horizontal axis must be $\sigma _x=\sigma \cos ^2\theta _1$, which for $\theta _1=0$ makes no difference. (c) End result once the ‘$\theta _2$ closest to $\theta _1$’ filter has been applied. This is the result found in Bret & Narayan (2018) and Haggerty et al. (2022).

Could the evolutionarity criterion help filtering them (the answer is ‘no’)? Figure 5(b) answers the question. Here is plotted the downstream Mach number $M_{A2x}$ of each solution, as defined by (2.6), where the Alfvén speed has been corrected according to (4.1). Our switch-on solutions have exactly $M_{A2x} = 1$, while the others also satisfy the evolutionary criterion.⁶ In line with the adequate definition (2.6) of the Alfvén Mach number for evolutionarity analysis, we use on the horizontal axis the parameter

(5.1)\begin{equation} \sigma_x \equiv \sigma \cos^2\theta_1, \end{equation}

which, for the present case $\theta _1=0$, makes no difference.

Note that, with definition (4.1) of the Alfvén speed, $M_{A2x}=+\infty$ on the firehose threshold. This is the reason why the green branch on the left is sent to $M_{A2x}=+\infty$ on the centre plot. As a consequence, each time the system eventually settles in stage 2-firehose with $\sigma _x < 1$ (i.e. $M_{A1x}>1$), it is evolutionary.

As a conclusion, in the interval $\sigma \in [0.5,1]$, there are three solutions for stage 1, and all three satisfy the evolutionarity criterion. The evolutionarity criterion by itself is therefore not sufficient to trim the number of solutions down to 1, or even 0.

5.2. Particle-in-cell simulations for $\theta _1=0$

In order to check which solution the system chooses, we ran a series a PIC simulations for various $\sigma$ values, performed with the fully kinetic three-dimensional PIC code, TRISTAN-MP (Buneman 1993; Spitkovsky 2005).⁷ The surest way to tell whether the shock is of the switch-on type or not, is to plot the perpendicular component $B_{\perp 2}$ of the field in the downstream. According to (B2), with $\theta _1=0$ it simply reads

(5.2)\begin{equation} B_{{\perp} 2} = B_1 \tan \theta_2. \end{equation}

Figure 6(a) shows the expected value of $B_{\perp 2}/B_1$, in red for the switch-on solutions with $\theta _2\neq 0$, and in black for the solution with $\theta _2=0$. The circles show the values of $\sigma$ which have been simulated. Figure 6(b) shows the results of the PIC simulations, with the shock density profile on top, and the ratio $B_{\perp 2}/B_1$ below. Besides some perturbations around the shock front at low $\sigma$, $B_{\perp 2}/B_1$ does not depart from 0 in the downstream for any $\sigma$. These perturbations are due to particle acceleration and back-reaction in the precursor, and from the instabilities triggered by the interaction of the fast upstream flow with the front (Sironi, Spitkovsky & Arons 2013; Brown et al. 2023).

Figure 6. (a) Expected value of $B_{\perp 2}/B_1$, in red for the switch-on solutions with $\theta _2\neq 0$, and in black for the solution with $\theta _2=\theta _1=0$. The circles show the values of $\sigma$ which have been simulated. (b) Results of the PIC simulations, with the shock density profile on top, and the ratio $B_{\perp 2}/B_1$ below. The $\sigma$ value for a curve is given by the circle of the same colour on the left plot. The same results are obtained for $\theta _1=2^\circ$. The flow is along the $x$ axis and $d_e=c/\omega _p$ is the electron inertial length.

Simply put, PIC simulations consistently discard the switch-on solutions. The same results are obtained for $\theta _1=2^\circ$, so that we are not witnessing a singular behaviour fruit of a perfect, hence unrealistic, symmetry. A similar pattern was observed in the relativistic regime in Bret et al. (2017).

The situation for $\theta _1 \sim 0$ is here markedly different from the MHD case. In MHD, where isotropy is imposed, the evolutionarity criterion imposes switch-on shocks within some $\sigma$-range, as explained in § 3. In our model, where an anisotropy is driven by the field, the evolutionarity criterion does not impose switch-on solutions, while PIC simulations consistently choose the non-switch-on ones.

Arguably, this explains why so few detections of switch-on shocks have been made in the Solar System. Indeed, among the thousands of shocks observed (Farris et al. 1994; Russell & Farris 1995; Balogh & Treumann 2013, § 2.3.6.) only one ‘possible’ detection of an interplanetary switch-on shock was reported in Feng et al. (2009). Also, Russell & Farris (1995) reported the detection of only one switch-on shock among the International Sun-Earth Explorer data.⁸

Still, what about these exceptions, since our collisionless scenario should rule them out? Several explanations are possible:

1. • The detections were faulty. Hence the adjective ‘possible’ associated with one of them.

2. • We only solve our model for a cold upstream, namely $T_1=0$. Maybe a finite $T_1$ would have our model allowing for some switch-on shocks.

3. • The ansatz (2.3) is not accurate enough, and a better version would allow for some switch-on shocks.

At any rate, PIC simulations and observations tell switch-on shocks in collisionless plasma are rare. We therefore propose the following criterion allowing us to choose between various solutions: the system chooses the solution with $\theta _2$ closest to $\theta _1$.

Note that this ‘$\theta _2$ closest to $\theta _1$’ criterion stems from our PIC simulations, as well as others of parallel shocks in pair (Bret et al. 2017; Haggerty et al. 2022) or electron/ion (Niemiec et al. 2012; Zeković 2019) plasmas, where the same, non-switch-on branch, is consistently chosen. This is why we used the verb ‘propose’. Its robustness on longer time scales, or other shock geometries, is beyond the scope of this work and should be assessed in further works.

Figure 5(c) eventually shows the end result once the ‘$\theta _2$ closest to $\theta _1$’ filter has been applied to the $\theta _1=0$ case. There is now but one solution for any $\sigma$, which indeed is the one found in Bret & Narayan (2018) and Haggerty et al. (2022).

5.3. General algorithm for branches selection

We may now lay out the general algorithm for branches selection. The criteria used to eliminate some are, applied in this order:

1. (a) Exclude stage 1 solutions with $r<1$. Then select the one with the $\theta _2$ closest to $\theta _1$.

2. (b) In case stage 1 is unstable, is the anisotropy $A_2$ of a given stage 2 solution, negative? This was already implemented in Paper 1 and allows us to eliminate some stage 2 -firehose and -mirror solutions. This never happens with stage 1 since it has, by design from (2.3), $A_2 = \frac {1}{2}\tan ^2\theta _2$.

3. (c) Does the resulting solution fulfil the evolutionarity criterion ?

The evolutionarity criterion is applied last because it operates on time scales related to the propagation of the shock,⁹ whereas the other criteria operate on much shorter time scales, related to plasma instabilities.

The algorithm is eventually represented by the flowchart in figure 7. In case stage 1 is mirror unstable for a given pair $(\sigma,\theta _1)$, there can be a degeneracy in the choice of the stage 2-mirror states for the same pair $(\sigma,\theta _1)$. In such a case, we choose the stage 2-mirror with the $\theta _2$ closest to the $\theta _2$ of the stage 1 it comes from, as was already done in Bret & Narayan (2022b). Such a situation never happens with stage 2-firehose.

Figure 7. Flowchart of the resolution of our model. In case stage 1 is mirror unstable for a pair $(\sigma,\theta _1)$, there can be a degeneracy in the choice of stage 2-mirror for the same pair $(\sigma,\theta _1)$. In such a case, we choose the stage 2-mirror with the $\theta _2$ closest to the $\theta _2$ of the stage 1 it comes from, as was already done in Bret & Narayan (2022b). Such a situation never happens with stage 2-firehose.

According to the flowchart on figure 7, there is necessarily only 1 solution left before applying the evolutionarity criterion. There is therefore no need to apply the ‘$\theta _2$ closest to $\theta _1$’ filter in (c), since only 1 branch can make it to this stage. Yet, this does not mean the evolutionarity filter is useless, as it can simply forbid the existence of a solution in some $\sigma$ range, as is the case in MHD (see bottom row of figure 2).

6. Results

While the fruit of our algorithm has already been explained for $\theta _1=0$, it is interesting to detail how it unfolds for an oblique field, like for example $\theta _1 = 0.3{\rm \pi} /2$.

Figure 8(a) shows all possible density jumps for stage 1. Notably, $r < 1$ for $\sigma \in [1, 1.3]$. The corresponding branches are then eliminated in figure 8(b).

Figure 8. Detailed progression of the application of the filtering process described by the flow chart in figure 7, for $\theta _1=0.3{\rm \pi} /2$. In (d), the green colour indicates that these stage 1 solutions are firehose unstable. Then in (g), the green colour shows the density jump of stage 2-firehose. The red branch in ( f) with $A_2 \sim -1$ discards stage 2 firehose solutions with $r >4$, not shown in (e).

For these stage 1 solutions with $r > 1$, figure 8(c) then shows the corresponding values of $\theta _2$. The horizontal red line is at $\theta _2=\theta _1$. The result of the ‘$\theta _2$ closest to $\theta _1$’ selection rule is displayed in figure 8(d).

In figure 8(d), the green colour indicates that some stage 1 solutions at low $\sigma$ are firehose unstable. We therefore need to examine stage 2-firehose solutions. All of them are displayed in figure 8(e). Yet, figure 8( f) shows that the lower branch has $A_2 < 0$ and needs to be eliminated. As a consequence, the only physical stage 2-firehose solution available when stage 1 is firehose unstable, is the upper branch.

Replacing then the firehose unstable stage 1 solutions of figure 8(d), by the corresponding stage 2-firehose solutions, gives figure 8(g).

We finally need to apply the evolutionarity test to the solutions of figure 8(g). This is done in figure 8(h) where the downstream Alfvén Mach number $M_{A2x}$ has been computed for the solution plotted in figure 8(g), with the forbidden zones shaded. Note that, for such an analysis, the horizontal scale has to be $\sigma _x=\sigma \cos ^2\theta _1$. Then and only then has the evolutionary analysis some branches passing exactly through the point $(1,1)$, as in figure 8(h).

The stage 2-firehose branch visible in green in figure 8(g) at small $\sigma$ passes the evolutionarity test since it has $\sigma _x < 1$ and $M_{A2x} = +\infty$. The analysis shows that some stable stage 1 solutions do not pass the evolutionarity test. As a result, figure 9 for $\theta _1=0.3{\rm \pi} /2$ displays a gap without solution for $\sigma \in [1,1.17]$, that does not show in figure 8(g).

Figure 9. Result of the filtering process elaborated in this work for various angles $\theta _1 \in [0, {\rm \pi}/2]$. All of them but $\theta _1=0.3{\rm \pi} /2$ are the ones considered in figure 8 of Bret & Narayan (2022b). The grey dashed lines picture the MHD result, evolutionary filtered. The green solutions pertain to (stage 2-firehose, green), and the red ones to (stage 2-mirror, red).

The end result of such a filtering process is eventually shown in figure 9 for several angles $\theta _1$ between 0 and ${\rm \pi} /2$. For comparison, all of them but $\theta _1=0.3{\rm \pi} /2$, are the ones considered in figure 8 of Bret & Narayan (2022b). The grey dashed line pictures the MHD result, evolutionary filtered.

Figure 9 eventually interpolates between the parallel case treated in Bret & Narayan (2018), and the perpendicular one treated in Bret & Narayan (2019). While § 5 showed that some filtering is needed to get to the end result for the parallel case, no filtering at all is required for the perpendicular case. Figure 9 for $\theta _1=0.9{\rm \pi} /2$ gives the result already found in Bret & Narayan (2019), without any filtering but the $A_2 < 0$ for the mirror solutions. The main reason for this is that for the perpendicular case, there is but one stage 1 solution, which, when stable, fulfils the evolutionarity criterion.

7. Conclusion

Applying the MHD formalism to analyse collisionless shocks may be problematic in the presence of an ambient magnetic field, capable of stabilizing pressure anisotropies. In a series of recent papers, we developed a method allowing us to determine the downstream anisotropy of a collisionless shock (Bret & Narayan 2018, 2019, 2022a,b; Bret 2023a,b). The anisotropy can then be included in the MHD conservation equations for anisotropic pressures, and the modified density jump derived. The case of a parallel shock has been successfully tested by PIC simulations in Haggerty et al. (2022), confirming that for a strong enough field, the density jump can go from 4, the expected MHD value, to only 2, the anisotropy corrected one.

Once the density jump is found, all the other jumps such as pressure, temperature, entropy, etc. can be straightforwardly derived (Bret 2021).

As can be seen in figure 9, our results differ from the isotropic MHD ones in 2 ways:

1. • The ranges of $\sigma$ without solutions differ.

2. • For $\sigma$ values with a solution, our density jump is usually, although not always, lower than the MHD one.

Overall, our results and the MHD ones bear a ‘family resemblance’, the largest discrepancy being found for the parallel case $\theta _1=0$. The predicted large reduction of the density jump could have important consequences for particle acceleration, since their index scales like $(r-1)^{-1}$ (Axford, Leer & Skadron 1977; Bell 1978a,b; Blandford & Ostriker 1978; Caprioli, Haggerty & Blasi 2020; Haggerty & Caprioli 2020).

Besides the strong sonic shock assumption of the present work, namely $P_1=0$, its main limitation may lie in the composition of the plasma, here a pair plasma. As stated in the introduction, such an assumption allows us to consider only one parallel and one perpendicular temperature. Yet, our formalism being eventually macroscopic, we expect our conclusions to hold for electron/ion plasmas as well. Such a conjecture is currently being tested through PIC simulations of such plasmas (Shalaby 2024).