1. Introduction
Shock waves are fundamental processes in fluids. They have been the subject of numerous studies for nearly two centuries (Salas Reference Salas2007). 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 Reference Zel'dovich and Raizer2002).
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 Reference Bale, Mozer and Horbury2003; Schwartz et al. Reference Schwartz, Henley, Mitchell and Krasnoselskikh2011). Such a shock cannot be mediated by collisions. It is mediated by collective plasma electromagnetic effects (Sagdeev Reference Sagdeev1966). 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 Reference Gurnett and Bhattacharjee2005; Goedbloed, Keppens & Poedts Reference Goedbloed, Keppens and Poedts2010; Thorne & Blandford Reference Thorne and Blandford2017). 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 Reference Bret2020). The first, because in the absence of collisions, a plasma can maintain a stable anisotropy in the presence of an ambient magnetic field (Hasegawa Reference Hasegawa1975; Gary Reference Gary1993). 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 Reference Drury1983; Blandford & Eichler Reference Blandford and Eichler1987).
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 Reference Weibel1959; Silva, Afeyan & Silva Reference Silva, Afeyan and Silva2021). 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 Reference Hudson1970; Karimabadi, Krauss-Varban & Omidi Reference Karimabadi, Krauss-Varban and Omidi1995; Erkaev, Vogl & Biernat Reference Erkaev, Vogl and Biernat2000; Gerbig & Schlickeiser Reference Gerbig and Schlickeiser2011). 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 Reference Bret and Narayan2018), 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 Reference Haggerty, Bret and Caprioli2022).
The model was then applied to the case of a perpendicular shock (Bret & Narayan Reference Bret and Narayan2019), and finally to the cases of a switch-on and of an oblique shock (Bret & Narayan Reference Bret and Narayan2022a,Reference Bret and Narayanb). 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 (Reference Bret and Narayan2022b), 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$.
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 (Reference Bret and Narayan2018) 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, energyFootnote 1 ) are therefore completed by
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.Footnote 2 Imposing marginal stability then provides again a fourth equation allowing us to fully determine the properties of the downstream.
Bret & Narayan (Reference Bret and Narayan2018), 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 Reference Guo, Sironi and Narayan2017, Reference Guo, Sironi and Narayan2018). 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 Reference Shalaby2024).
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. (Reference Haggerty, Bret and Caprioli2022). Such a result stands in contrast with MHD, where the density jump should always be 4, regardless of the field strength.Footnote 3
The perpendicular case was treated in Bret & Narayan (Reference Bret and Narayan2019). There, the direction perpendicular to the flow is eventually parallel to the field, so that the counterpart of (2.1) is
The switch-on shocks, where the field is oblique in the downstream only, was treated in Bret & Narayan (Reference Bret and Narayan2022a). The model has also been solved for a parallel or a perpendicular shock, with an anisotropic upstream (Bret Reference Bret2023a,Reference Bretb).
Finally, the general case where the field may be oblique in both the upstream and the downstream was treated in Bret & Narayan (Reference Bret and Narayan2022b), namely, Paper 1.
In Bret & Narayan (Reference Bret and Narayan2022a,Reference Bret and Narayanb), 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
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 Reference Landau and Lifshitz1981, § 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:
• 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.
• The resulting state of the plasma in the downstream is called ‘stage 1’.
• If the downstream field is strong enough to stabilize stage 1, then this is the end state of the downstream.
• If the downstream field is too weak to stabilize stage 1, then
• If stage 1 is firehose unstable, it migrates to the firehose instability threshold. This is ‘stage 2 – firehose’.
• 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:
(a) It is restricted to strong sonic shocks, namely $T_1=0$, and to non-relativistic pair plasmas. These restrictions are still considered here.
(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.
(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:
where $v_{A,i}$ is the Alfvén velocity
The parameter $\sigma$ is frequently used in simulations of collisionless shocks like Haggerty et al. (Reference Haggerty, Bret and Caprioli2022), 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:
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
is a root of the polynomial (A9) in Appendix A. The MHD density jump is then given by
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]$.
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 Reference Kennel, Blandford and Wu1990; Farris et al. Reference Farris, Russell, Fitzenreiter and Ogilvie1994; Falle & Komissarov Reference Falle and Komissarov2001; Kulsrud Reference Kulsrud2005). 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 Reference Falle and Komissarov1997; Kulsrud Reference Kulsrud2005). 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 Reference Falle and Komissarov1997) or ‘extraneous’ (Kulsrud Reference Kulsrud2005) 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
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. Reference Goedbloed, Keppens and Poedts2010, 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. Reference Sa̧dowski, Narayan, Tchekhovskoy and Zhu2013, Reference Sa̧dowski, Narayan, McKinney and Tchekhovskoy2014). 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 (Reference Kulsrud2005) for example.
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 (Reference Hudson1970) and are reported in Appendix B with the notations of Bret & Narayan (Reference Bret and Narayan2022a,Reference Bret and Narayanb).
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 (Reference Abraham-Shrauner1967) that the Alfvén velocity readsFootnote 4
where $v_A$ is still given by (2.5) and
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
whereFootnote 5
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. Reference Hellinger, Trávníček, Kasper and Lazarus2006; Bale et al. Reference Bale, Kasper, Howes, Quataert, Salem and Sundkvist2009; Maruca, Kasper & Bale Reference Maruca, Kasper and Bale2011).
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 Reference Hasegawa1975; Gary Reference Gary1993)
Yet, a different criterion is given in Abraham-Shrauner (Reference Abraham-Shrauner1967), namely
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.
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$.
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.Footnote 6 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
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 Reference Buneman1993; Spitkovsky Reference Spitkovsky2005).Footnote 7 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
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 Reference Sironi, Spitkovsky and Arons2013; Brown et al. Reference Brown, Juno, Howes, Haggerty and Constantinou2023).
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. (Reference Bret, Pe'er, Sironi, Sa̧dowski and Narayan2017).
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. Reference Farris, Russell, Fitzenreiter and Ogilvie1994; Russell & Farris Reference Russell and Farris1995; Balogh & Treumann Reference Balogh and Treumann2013, § 2.3.6.) only one ‘possible’ detection of an interplanetary switch-on shock was reported in Feng et al. (Reference Feng, Lin, Chao, Wu, Lyu and Lee2009). Also, Russell & Farris (Reference Russell and Farris1995) reported the detection of only one switch-on shock among the International Sun-Earth Explorer data.Footnote 8
Still, what about these exceptions, since our collisionless scenario should rule them out? Several explanations are possible:
• The detections were faulty. Hence the adjective ‘possible’ associated with one of them.
• 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.
• 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. Reference Bret, Pe'er, Sironi, Sa̧dowski and Narayan2017; Haggerty et al. Reference Haggerty, Bret and Caprioli2022) or electron/ion (Niemiec et al. Reference Niemiec, Pohl, Bret and Wieland2012; Zeković Reference 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 (Reference Bret and Narayan2018) and Haggerty et al. (Reference Haggerty, Bret and Caprioli2022).
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:
(a) Exclude stage 1 solutions with $r<1$. Then select the one with the $\theta _2$ closest to $\theta _1$.
(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$.
(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,Footnote 9 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 (Reference Bret and Narayan2022b). 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).
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).
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 (Reference Bret and Narayan2022b). The grey dashed line pictures the MHD result, evolutionary filtered.
Figure 9 eventually interpolates between the parallel case treated in Bret & Narayan (Reference Bret and Narayan2018), and the perpendicular one treated in Bret & Narayan (Reference Bret and Narayan2019). 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 (Reference Bret and Narayan2019), 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 Reference Bret and Narayan2018, Reference Bret and Narayan2019, Reference Bret and Narayan2022a,Reference Bret and Narayanb; Bret Reference Bret2023a,Reference Bretb). 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. (Reference Haggerty, Bret and Caprioli2022), 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 Reference Bret2021).
As can be seen in figure 9, our results differ from the isotropic MHD ones in 2 ways:
• The ranges of $\sigma$ without solutions differ.
• 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 Reference Axford, Leer and Skadron1977; Bell Reference Bell1978a,Reference Bellb; Blandford & Ostriker Reference Blandford and Ostriker1978; Caprioli, Haggerty & Blasi Reference Caprioli, Haggerty and Blasi2020; Haggerty & Caprioli Reference Haggerty and Caprioli2020).
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 Reference Shalaby2024).
Acknowledgements
Thanks are due to A. Velikovich and F. Fraschetti for valuable inputs.
Editor Luís O. Silva thanks the referees for their advice in evaluating this article.
Funding
A.B. acknowledges support by the Ministerio de Economía y Competitividad of Spain (Grant No. PID2021-125550OB-I00). R.N. acknowledges support from the NSF Grant No. AST-1816420. C.C.H. was partially supported by NSF FDSS grant AGS-1936393 as well as NASA grants 80NSSC20K1273 and 80NSSC23K0099.
Simulations were performed on computational resources provided by TACC's Stampede2 and Purdue's ANVIL through ACCESS (formally XSEDE) allocation TG-AST180008. Some of the work was supported by the Geospace Environment Modeling Focus Group “Particle Heating and Thermalization in Collisionless Shocks in the Magnetospheric MultiScale Mission (MMS) Era”.
R.N. and A.B. thank the Black Hole Initiative (BHI) at Harvard University for support. The BHI is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Isotropic MHD conservation equations for an oblique shock
The MHD conservation equations for an oblique shock and a fluid of adiabatic index $\gamma =5/3$ read (Kulsrud Reference Kulsrud2005)
where $m$ is the mass of the particles, $k_B$ the Boltzmann constant and
After some algebra, $\bar {T}_2 \equiv \tan \theta _2$ is found as the solution ofFootnote 10
with
Appendix B. Anisotropic MHD conservation equations for an oblique shock
The conservation equations for anisotropic temperatures were derived in Hudson (Reference Hudson1970) and Erkaev et al. (Reference Erkaev, Vogl and Biernat2000). They have been re-derived in Bret & Narayan (Reference Bret and Narayan2022a) with the present notations. They are formally valid even for anisotropic upstream temperatures, with $T_{\parallel 1} \neq T_{\perp 1}$. Writing them for $T_{\parallel 1}=T_{\perp 1}\equiv T_1$, they read
where
It can be checked that setting $T_{\parallel 2} = T_{\perp 2} = T_2$ gives back the MHD equations, ((A1)–(A6)).
Appendix C. Main quantities for stages 1 and 2
Implementing the algorithm described in figure 7 requires computing the density ratio $r$, the anisotropy $A_2$ and the downstream Alfvén Mach number $M_{A2x}$, for stage 1, stage 2-firehose and stage 2-mirror. The results are presented below.
C.1 Stage 1
Solving the system of ((B1)–(B6)) with prescription (2.3) allows us to derive a polynomial for the quantify $T_2=\tan \theta _2$ defined in (3.1). It has been derived in Paper 1.Footnote 11 Using now the expression (4.1) for the Alfvén velocity, the Alfvén Mach number for stage 1 reads
From (2.3), the anisotropy of stage 1 is simply
The density ratio $r$ is the same as in Paper 1, namely,
where
C.2 Stage 2 firehose
If stage 1 has $A_2 < 1-2/\beta _{\parallel 2}$, then stage 2-firehose is the end state. Since the stability criterion differs from that used in Paper 1 by the factor 2, the properties of stage 2 firehose change with respect to Paper 1. A polynomial equation can still be derived for $T_2=\tan \theta _2$ as
with
The density jump is then given by
As for the Alfvén Mach number, we found in § 4 that the Alfvén speed vanishes on the firehose instability threshold. Therefore, stage 2-firehose has
As a consequence, when the system eventually settles in stage 2-firehose with $\sigma _x < 1$ (i.e. $M_{A1x}>1$), it is evolutionary.
The anisotropy of stage 2-firehose is also modified with respect to Paper 1, due to the modified stability criterion. It now reads
C.3 Stage 2 mirror
Since the stability threshold for the mirror instability is the same as in Paper 1, the polynomial for $T_2$, the anisotropy $A_2$ and the density ratio are also the same. The Alfvén Mach number reads here