Self-organization of collisionless shocks: a ‘phase transition’ from a planar stationary profile to a rippled structure

Abstract

It is common wisdom that collisionless shocks become non-planar and non-stationary at sufficiently high Mach numbers. Whatever the shock structure, the upstream and downstream fluxes of the mass, momentum and energy should be equal. At low Mach numbers, these conservation laws are satisfied when the shock front is planar and stationary. When this becomes impossible, inhomogeneity and time dependence, presumably in the form of rippling, develop. In this study, we show that the shock structure changes as a kind of ‘phase transition’ when the Mach number is increased while other shock parameters are kept constant.

1 Introduction

Collisionless shocks, particularly fast magnetized shocks, play an important role in particle heating and acceleration in numerous systems in space plasmas. The magnetic structure of the shock front largely determines what happens to the charged particles, ions and electrons upon crossing the shock. This structure can be studied in detail only within the heliosphere, where in situ measurements are possible. Plenty of heliospheric observations (Montgomery, Asbridge & Bame 1970; Greenstadt et al. 1975, 1977, 1980; Livesey, Kennel & Russell 1982; Russell, Hoppe & Livesey 1982a; Russell et al. 1982b, 1983; Bavassano-Cattaneo et al. 1986; Mellott & Livesey 1987; Farris, Russell & Thomsen 1993; Tatrallyay et al. 1997; Bale et al. 2005; Burgess et al. 2005; Moullard et al. 2006; Lobzin et al. 2007, 2008; Lefebvre et al. 2009; Krasnoselskikh et al. 2013; Masters et al. 2013; Johlander et al. 2016; Gingell et al. 2017; Johlander et al. 2018; Dimmock et al. 2019; Liu et al. 2021) show that the shock structure depends on the Mach number, on the angle $\theta _{{\rm Bn}}$ between the upstream magnetic field and the shock normal and on the ratio of the upstream plasma pressure to the upstream magnetic pressure $\beta$. Quasi-parallel shocks, $\theta _{{\rm Bn}}<45^\circ$, are often thought to have an extended irregular shock front, although this is not universal (Burgess et al. 2005; Dimmock et al. 2023; Jebaraj et al. 2024). Oblique shocks typically have a well-defined shock transition. Very low-Mach-number shocks possess a nearly monotonic magnetic profile (Farris et al. 1993). Higher-Mach-number shocks have a substantial overshoot (Russell et al. 1982a), and even higher-Mach-number shocks are rippled (Moullard et al. 2006) or reforming (Lobzin et al. 2007; Lefebvre et al. 2009). However complicated is the shock structure, the mass, momentum and energy must be conserved, that is, the upstream and downstream fluxes of the mass, momentum and energy, properly averaged over space and time, should be equal. In other words, the developed shock structure should be able to ensure the equality of the fluxes. If a planar stationary structure is sufficient to maintain the conservation laws, the shock would be planar and stationary. If a substantial overshoot is required but otherwise planarity and stationarity are enough, an overshoot will grow (Gedalin & Sharma 2023; Gedalin et al. 2023a; Sharma & Gedalin 2023). If this is not sufficient, the shock will become rippled or even reforming, whatever the exact mechanism of the development of spatial inhomogeneity and time dependence. This is the principle of the shock self-organization. The structural changes in the shock front can be expected to occur via a series of ‘phase transitions’. To illustrate the principle, we analyse whether the conservation laws can be satisfied with a planar stationary profile, for a chosen set of the shock angle and upstream $\beta$ and gradually increasing Mach number.

2 Definitions

For simplicity, the plasma is assumed to be composed of ions $i$ (protons $p$) and electrons $e$. The upstream variables are denoted by the subscript $u$ and the downstream variables are denoted by the subscript $d$. The normal incidence frame (NIF) is the shock frame in which the upstream plasma velocity is along the shock normal. The NIF upstream plasma speed is $V_u$. The upstream Alfvén speed is $v_A=B_u/4{\rm \pi} n_um_p$, where $B_u$ is the upstream magnetic field magnitude, $n_u$ is the upstream proton number density (which is equal to the upstream electron number density) and $m_p$ is the proton mass. The Alfvén Mach number is $M_A=V_u/v_A$. The coordinates are chosen so that the upstream magnetic field is $\boldsymbol {B}_u=(B_{ux}, B_{uy}, B_{uz})=B_u(\cos \theta _{{\rm Bn}},0,\sin \theta _{{\rm Bn}})$. We also define $\beta _i=8{\rm \pi} n_uT_{iu}/B_u^2$ and $\beta _e=8{\rm \pi} n_uT_{eu}/B_u^2$, where $T_{iu}$ and $T_{eu}$ are the ion (proton) and electron upstream temperatures, respectively.

3 The method

The test particle analysis used here is described in detail by Gedalin (2016). Here, we briefly describe the principles of the method. This is not a simulation but should be considered as a numerical tool of theory. In the test particle analysis, the equations of motion of charged particles (ions here) are solved numerically in the electromagnetic fields prescribed by an appropriate model. At present, plenty of models of subcritical and supercritical shocks, including rippling and reformation, have been developed and analysed (Balikhin et al. 2008; Gedalin, Friedman & Balikhin 2015; Gedalin 2016, 2019a,b; Gedalin, Dröge & Kartavykh 2016; Dimmock et al. 2023; Gedalin, Pogorelov & Roytershteyn 2023b). The ion velocities are acquired at a dense grid, using the staying time method (see technical details in Gedalin (2016) and Gedalin et al. (2023b)), and all relevant moments of the ion distribution are calculated as functions of the position in a large space surrounding the shock transition. The model profiles used in the analysis typically depend on a small number of control parameters. This test particle analysis has several applications. It has been used to separately determine the effects of various parameters (e.g. magnetic compression and cross-shock potential) on the features and relaxation of the downstream ion distributions (see e.g. Gedalin 2015), and to separate the effect of macroscopic fields from the effect of fluctuations (Gedalin et al. 2023a). It has been used to model observed shocks, and the predictions of the analysis were found to be in surprisingly good agreement with observations (Gedalin et al. 2018; Pope, Gedalin & Balikhin 2019). The method has been used for analysing whether a supercritical shock can be stationary (Gedalin 2019a). All these tasks require the usage of the non-self-consistent approach, where various parameters can be varied independently. Here, we use the test particle analysis to study whether the conservation laws can be fulfilled if the shock profile is planar and stationary, namely the ion number density conservation:

(3.1)\begin{equation} nV_x=n_uV_u,\end{equation}

and the conservation of the momentum along the shock normal, also known as pressure balance:

(3.2)\begin{equation} p_{i,xx}+n_uT_{eu}\left(\frac{n}{n_u}\right)^{5/3}+ \frac{B^2}{8{\rm \pi}}=n_um_pV_u^2+n_uT_{eu}+n_uT_{iu}+\frac{B_u^2}{8{\rm \pi}}.\end{equation}

The procedure is as follows. First, we prescribe the model macroscopic fields $\boldsymbol {B}^{{\rm (mod)}}$ and $\boldsymbol {E}^{{\rm (mod)}}$ of the shock, depending only on the coordinate $x$ along the shock normal. Next, we trace ions across this shock by numerically solving their equations of motion in the prescribed fields. This solution provides us with the ion velocities at all positions across the shock and, thus, allows us to numerically calculate the number density $n^{{\rm (der)}}(x)$, the hydrodynamic velocity $\boldsymbol {V}^{{\rm (der)}}(x)$ and the total (dynamic and kinetic) pressure $p_{i,xx}^{{\rm (der)}}(x)$, as functions of the position $x$. Electrons are treated as a massless neutralizing fluid with the equation of state $p_e\propto n^{5/3}$. The approach maintains (3.1) automatically. Using (3.2) we derive

(3.3)\begin{gather} B^{{\rm (der)}}=\left[8{\rm \pi}\left(W-p_{i,xx}^{{\rm (der)}}-n_uT_{eu} \left(\frac{n^{{\rm (der)}}}{n_u}\right)^{5/3}\right)\right]^{1/2}, \end{gather}

(3.4)\begin{gather}W=n_um_pV_u^2+n_u(T_{eu}+T_{iu})+\frac{B_u^2}{8{\rm \pi}}, \end{gather}

for the chosen upstream parameters $n_u$, $V_u$, $T_{eu}$, $T_{iu}$ and $B_u$. The derived magnetic field magnitude $B^{{\rm (der)}}$ is ultimately compared with the magnitude $B^{{\rm (mod)}}$ of the model magnetic field used for the ion tracing. If there is substantial disagreement between the two, we conclude that the model is not satisfactory.

The shock parameter space is multi-dimensional, so it is reasonable to keep most of the parameters constant and vary only one or two. Here, the shock angle is kept constant, $\theta _{{\rm Bn}}=60^\circ$, for all Mach numbers tested during the study. We also keep constant $\beta _i=\beta _e=0.5$. The basic shock profile is chosen as follows (Gedalin 2016):

(3.5)\begin{gather} B_z^{{\rm (mod)}}=B_u\sin\theta_{{\rm Bn}}\left(\frac{R+1}{2} +\frac{R-1}{2}\tanh\frac{3x}{D}\right), \end{gather}

(3.6)\begin{gather}B_x^{{\rm (mod)}}=B_u\cos\theta_{{\rm Bn}}, \end{gather}

(3.7)\begin{gather}B_y^{{\rm (mod)}}=\frac{c\cos\theta_{{\rm Bn}}}{M_A\omega_{{\rm pi}}} \frac{{\rm d}B_z^{{\rm (mod)}}}{{\rm d} x}, \end{gather}

where the coordinate $x$ is along the shock normal. The parameter $R$ is related to the magnetic compression via

(3.8)\begin{equation} \frac{B_d}{B_u}=\sqrt{R^2\sin^2\theta_{{\rm Bn}} +\cos^2\theta_{{\rm Bn}}},\end{equation}

and the NIF electric field is $E_x^{{\rm (mod)}}=-K_E({\rm d}B_z^{{\rm (mod)}}/{{\rm d} x})$, $E_y^{{\rm (mod)}}=V_uB_u\sin \theta _{{\rm Bn}}/c$, $E_z^{{\rm (mod)}}\,{=}\,0$. The coefficient $K_E$ is determined by the value of the cross-shock potential:

(3.9)\begin{equation} e\phi^{{\rm (mod)}}={-}\int E_x^{{\rm (mod)}}{{\rm d} x}=s(m_pV_u^2/2), \end{equation}

where $s$ is one of the varied parameters of modelling. The ramp width is chosen as $D\,{=}\,c/\omega _{{\rm pi}}$, where the ion plasma frequency is $\omega _{{\rm pi}}=(4{\rm \pi} n_ue^2/m_p)^{1/2}$.

The chosen model profile reproduces the observed nearly monotonic increase of the main magnetic field $B_z$ in the ramp, the non-coplanar magnetic field $B_y$ confined within the ramp and the cross-shock electric field $E_x$ inside the ramp (see e.g. Greenstadt et al. 1975, 1980; Scudder et al. 1986; Farris et al. 1993; Dimmock et al. 2012). The approximate relations between the field components were analytically derived (Jones & Ellison 1987; Gosling, Winske & Thomsen 1988; Schwartz et al. 1988; Gedalin 1996) and successfully used for comparison with observations (Gedalin et al. 2022b). Note that ions are not sensitive to fine details of the shock front or small-scale fields (Gedalin 2020).

The magnetic compression $B_d/B_u$ is not a free parameter but is taken from the solution of the Rankine–Hugoniot relations taking into account that the downstream pressure is not necessarily isotropic (Abraham-Shrauner 1967; Hudson 1970; Chao & Goldstein 1972; Sanderson 1976; Lyu & Kan 1986; Kennel 1988; Vogl et al. 2001; Livadiotis 2019; Gedalin et al. 2022a; Haggerty, Bret & Caprioli 2022). Table 1 contains the pairs $M_A, B_d/B_u$ for the chosen $\theta _{{\rm Bn}}$, $\beta =\beta _i+\beta _e$ and anisotropy ratio $A=p_{d,\parallel }/p_{d,\perp }$, where $\parallel$ and $\perp$ refer to the direction of the downstream magnetic field. If needed, an overshoot will be added to the profile (see below).

Table 1. Alfvénic Mach number $M_A$ versus magnetic compression $B_d/B_u$ (rounded up to first digit after the decimal point), according to Rankine–Hugoniot relations with various anisotropies $A=T_\parallel /T_\perp$.

Initially, there are 80 000 Maxwellian distributed ions starting their motion towards the shock from a position far upstream, where the flow velocity is $\boldsymbol {V}=(V_u,0,0)$. Catching ions crossing thin layers, the distribution function $f(\boldsymbol {v},x)$ is derived, from which the number density $n^{{\rm (der)}}=\int f(\boldsymbol {v},x)\, {\rm d}^3\boldsymbol {v}$ and the ion pressure $p^{{\rm (der)}}_{i,xx}=m_p\int v_x^2f(\boldsymbol {v},x)\, {\rm d}^3\boldsymbol {v}$ are calculated (Gedalin 2016).

We vary the parameters $s$ and $B_d/B_u$ (the latter within the fork given in the table) to achieve convergence of the far downstream $B^{{\rm (der)}}$ to the initially chosen $B_d$. If such convergence can be achieved, we further check whether the profiles $B^{{\rm (mod)}}(x)$ and $B^{{\rm (der)}}(x)$ agree reasonably.

4 Dependence on the Mach number

4.1 $M_A=2$

For a Mach number of $M_A=2.0$, the best agreement between the derived and model downstream fields is achieved for $s=0.3$ and $B_d/B_u=1.55$. The derived and model profiles are shown in figure 1. The magnetic compression is weaker than would be expected in the isotropic case. Indeed, in such ion tracings, the downstream ion distributions are significantly anisotropic (Gedalin et al. 2022a). Such anisotropies are observed at the Earth bow shock as well. The far downstream values agree very well. There exists a transitional region where weak oscillations of the magnetic field appear in the derived profile. These oscillations are related to the gradual kinematic collisionless relaxation of the gyrating directly transmitted ions behind the shock front (Balikhin et al. 2008; Ofman et al. 2009; Pope et al. 2019), which is seen in the reduced distribution function $f(x,v_x)$ in figure 2. The relative amplitude of these oscillations $|B-B_d|/B_d$ does not exceed 15 % and they are not expected to affect the ion motion, so there is no reason to try to improve the agreement by adding an overshoot.

Figure 1. The model $B^{{\rm (mod)}}(x)$ (black line) and derived $B^{{\rm (der)}}(x)$ (blue line) magnetic field magnitudes for $M_A=2$, $s=0.3$ and $B_d/B_u=1.55$. The far downstream values agree very well. There are weak downstream magnetic field oscillations in the derived profile.

Figure 2. The reduced distribution function $f(x,v_x)$ (log scale). The model magnetic field magnitude is shown by the black line.

4.2 $M_A=2.5$

For a Mach number of $M_A=2.5$, the amplitude of the oscillations becomes substantial, $|B-B_d|/B_d>30\,\%$, and it is necessary to take into account the overshoot and the undershoot in the model profile. Each oscillation is added using the localized function (Gedalin & Ganushkina 2022; Gedalin et al. 2023a,b)

(4.1)\begin{gather} \Delta B_z=aB_u\sin\theta_{{\rm Bn}} g(x,x_l,w_l,x_r,w_r), \end{gather}

(4.2)\begin{gather}g(x,x_l,w_l,x_r,w_r)=\left(1+\tanh\frac{3(x-x_l)}{w_l}\right) \left(1-\tanh\frac{3(x-x_r)}{w_r}\right), \end{gather}

where the parameters are adjusted for both the overshoot and the undershoot. The best agreement was achieved using $a=1$, $x_l =0.8$, $x_r= 0.8$, $w_l= 1.25$, $w_r = 0.95$ for the overshoot and $a=-0.95$, $x_l =1.55$, $x_r= 1.55$, $w_l= 0.95$, $w_r = 0.95$ for the undershoot, while $s=0.25$ and $B_d/B_u=1.86$. The addition to $B_y$ and $E_x$ is obtained following the same prescription as for the non-structured shock. The model and derived profiles are shown in figure 3. The far downstream values of the magnetic field magnitude agree well. The overshoot and undershoot also agree rather well. The corresponding reduced distribution function is shown in figure 4. Given the success of adjusting the structured profile, there is no reason to seek further improvements by modelling more oscillations. It is clear that a structured planar stationary profile is sufficient for maintaining pressure balance to good precision.

Figure 3. The model $B^{{\rm (mod)}}(x)$ (black line) and derived $B^{{\rm (der)}}(x)$ (blue line) magnetic field magnitudes for $M_A=2.5$, $s=0.25$ and $B_d/B_u=1.86$. The far downstream values agree well. The overshoot and undershoot also agree rather well.

Figure 4. The reduced distribution function $f(x,v_x)$ (log scale). The model magnetic field magnitude is shown by the black line.

4.3 $M_A=3$ and $M_A=3.5$

Until now, no noticeable ion reflection has occurred. That is, even if there was a very small fraction of reflected ions, they had no effect on the shock profile. The overshoot, undershoot and subsequent oscillations were produced due to the deceleration of the directly transmitted ions at the shock crossing and their postshock gyration and gradual gyrophase mixing. At $M_A=3$, ion reflection is already substantial. Yet, using the same procedure as above, it appears possible to achieve a good agreement between the derived and model profiles, as can be seen in figure 5. In this case, the best agreement was found to be at $s=0.65$ and $B_d/B_u=2.05$. The reflected ions are seen in figure 6. Note that the overshoot and other oscillations are still due to the deceleration and gyration of the directly transmitted ions, while the reflected ions play an important role in the regulation of the overshoot strength (Gedalin & Sharma 2023).

Figure 5. The model $B^{{\rm (mod)}}(x)$ (black line) and derived $B^{{\rm (der)}}(x)$ (blue line) magnetic field magnitudes for $M_A=3$, $s=0.65$ and $B_d/B_u=2.05$. The far downstream values agree well. The overshoot and undershoot also agree rather well.

Figure 6. The reduced distribution function $f(x,v_x)$ (log scale). The model magnetic field magnitude is shown by the black line.

For $M_A=3.5$ we can similarly adjust the profile with $B_d/B_u=2.25$ and $s=0.52$. The results are shown in figures 7 and 8. The differences between the two cases, $M_A=3$ and $M_A=3.5$, are quantitative and not qualitative. In both cases, a planar stationary structured shock profile is capable of ensuring the pressure balance to a good approximation. Note that in both cases, the presence of reflected ions is seen mainly in the increasing spread of the distribution function in the velocity space, that is, heating enhancement.

Figure 7. The model $B^{{\rm (mod)}}(x)$ (black line) and derived $B^{{\rm (der)}}(x)$ (blue line) magnetic field magnitudes for $M_A=3.5$, $s=0.52$ and $B_d/B_u=2.35$. The far downstream values agree well. The overshoot and undershoot also agree rather well.

Figure 8. The reduced distribution function $f(x,v_x)$ (log scale). The model magnetic field magnitude is shown by the black line.

4.4 $M_A=4$

It was not possible to achieve agreement on the far downstream values of the derived and modelled magnetic field, whatever structure modifications or parameter variations have been tried. Although it is always possible that we were not lucky enough in this try-and-check quest, there are good chances that this failure indicates the requirement of significant non-planarity and/or non-stationarity of the shock front to ensure that the conservation laws are satisfied. Since the shock is quasi-perpendicular, we expect the development of rippling (Gingell et al. 2017).

5 Discussion

For simplicity of the discussion, we assume perpendicular geometry, cold upstream plasma and isotropic downstream plasma and ignore electron heating. Then, the conservation laws give

(5.1)\begin{gather} \frac{n_d}{n_u}=\frac{B_d}{B_u}=\frac{V_u}{V_d}=R, \end{gather}

(5.2)\begin{gather}n_um_pV_u^2+\frac{B_u^2}{8{\rm \pi}}=n_dm_pV_d^2+\frac{B_d^2}{8{\rm \pi}}+n_dT_d, \end{gather}

(5.3)\begin{gather}\frac{1}{2}n_um_pV_u^3+\frac{B_u^2V_u}{4{\rm \pi}}=\frac{1}{2}n_dm_pV_d^3+ \frac{B_d^2V_d}{4{\rm \pi}}+\frac{5}{2}n_dT_dV_d. \end{gather}

Let $Y=T_d/m_pV_u^2$, then

(5.4)\begin{gather} 1+\frac{1}{2M_A^2}=\frac{1}{R}+\frac{R^2}{2M_A^2}+RY, \end{gather}

(5.5)\begin{gather}1+\frac{2}{M_A^2}=\frac{1}{R^2}+\frac{2R}{M_A^2}+5Y. \end{gather}

The expressions (5.4) and (5.5) determine both the magnetic compression $R$ and the ion heating $Y$ as functions of $M_A$, shown in figure 9.

Figure 9. The magnetic compression $R=B_d/B_u$ (left) and the dimensionless heating $Y=T_d/m_pV_u^2$ (right) as functions of the Mach number $M_A$, as derived from Rankine–Hugoniot relations in the case $\theta _{{\rm Bn}}=0$.

We should be able to obtain the same parameters from a model using additional parameters and/or structural elements. In a non-structured low-Mach-number shock, there are only directly transmitted ions. If the cross-shock potential is $\phi =s(m_pV_u^2/2e)$, the ion velocity upon crossing the shock would be $V_u(\sqrt {1-s},0,0)$ while the drift speed is $V_u/R$. The gyration speed is $v_g=V_u|\sqrt {1-s}-1/R|$ and the dimensionless downstream temperature $Y_{{\rm model}}=\frac {1}{2}(\sqrt {1-s}-1/R)^2$. Figure 10 shows that only for $s=0.4$ can the model dimensionless downstream temperature $Y_{{\rm model}}$ be made equal to the Rankine–Hugoniot required $Y$, and only for $M_A\lesssim 2.5$. If there were an overshoot, the maximum overshoot magnetic field could be estimated using $B_{\max }/B_u=\sqrt {2M_A^2(1-\sqrt {1-s}) +1}$ (Gedalin, Russell & Dimmock 2021). For $s=0.4$ and $M_A=2.5$, this estimate indicates that there is no overshoot or it is negligible.

Figure 10. The model dimensionless downstream temperature $Y_{{\rm model}}$ (blue) for $s=0.2,0.3,0.4, 0.5$ versus the Rankine–Hugoniot required $Y$ (black).

For higher Mach numbers, the presence of reflected ions is necessary to ensure the required heating to satisfy the Rankine–Hugoniot relations. The two main model parameters affecting the ion distributions and heating are the cross-shock potential and the maximum overshoot magnetic field. Ion reflection cannot be discussed for the cold upstream ion distribution since, in this case, all ions are either directly transmitted or reflected. For an actual thermal upstream distribution, reflected are those ions that are in the tail of the distribution. Therefore, the number of reflected ions is not large and sensitive to the cross-shock potential and the maximum magnetic field in the overshoot (Sharma & Gedalin 2023). The increase of each of these parameters increases the number of reflected ions and thus increases the contribution of the reflected ions in the downstream pressure. The cross-shock potential decreases the downstream gyration energy of the directly transmitted ions, thus decreasing their contribution to the total pressure. It also decreases the directly transmitted ion pressure upon the ramp crossing, where gyration has only started and no gyrophase averaging occurs. As a result, the magnetic pressure should substantially increase, resulting in the overshoot. The effect of the overshoot on the directly transmitted ions is in effectively increasing the local magnetic field in which they start to gyrate and decreasing the local drift speed so that in the overshoot–undershoot region, the ion pressure fluctuates. The minimum of the ion pressure corresponds to the maximum of the overshoot, while the maximum of the ion pressure leads to the magnetic dip in the undershoot. In this region, the directly transmitted ions are still gyrating as a beam so that the minimum and maximum pressures are closely related, and the growth of the overshoot is accompanied by the depression in the undershoot. The magnetic pressure cannot become negative, which limits the overshoot strength in the case of a planar stationary shock. When the Rankine–Hugoniot relations require too strong an overshoot and, therefore, too deep a magnetic depression in the undershoot, that means that a planar stationary profile is no longer capable of maintaining the pressure balance, and the shock has to become non-planar and/or time-dependent.

6 Simulations

For illustrative purposes, we have performed several low-cost two-dimensional hybrid simulations using the dHybrid code (Gargaté et al. 2007). A grid size of $2000\times 20$ was chosen, with cell spacing of $0.2(c/\omega _{{\rm pi}})$ in each direction. The longest dimension is along the shock normal, $x$, with the open boundary at $x=0$ and reflection boundary at $x=400(c/\omega _{{\rm pi}})$. Ions are injected at the open boundary, and a shock is produced due to the reflection of these ions off the reflection boundary. The shortest dimension, $z$, is along the main magnetic field, with periodic conditions at both boundaries. The shock angle is $\theta _{{\rm Bn}}=60^\circ$, that is, the initial magnetic field is $B_u(\sqrt {3}/2,0,1/2)$. Initially, $\beta _i=\beta _e=0.1$, lower than in the test particle analysis, to avoid unnecessary attempts to make detailed comparisons. The objective of the simulations is to observe the transition to rippling with an increase of the Mach number while keeping other parameters constant. The task does not require high precision, so we use 36 particles per cell and a time step of $10^{-3} \varOmega _u^{-1}$.

Figure 11 shows the reduced distribution function $f(x,v_x)$, on a log scale, at time $t=70 \varOmega _u^{-1}$ from the beginning of the run with the upstream plasma speed $v_{{\rm in}}=2v_A$ in the frame of the simulation. The shock speed in this frame is $v_{{\rm sh}}\approx 0.8v_A$, so the Mach number is $M_A\approx 2.8$. The downstream gyration of the directly transmitted ions is very clear, together with the gradual relaxation. There are no reflected ions. Figure 12 shows the two-dimensional patterns of $B_z(x,z)$ (figure 12a) and $B_x(x,z)$ (figure 12b), at the same $t=70 \varOmega _u^{-1}$. There are no signs of spatial inhomogeneity along the shock front in $B_z$. Fluctuations of the normal component of the magnetic field, $B_x(x,z)$, are the clearest signatures of the developing inhomogeneity (Gedalin & Ganushkina 2022). In this case, the level of these fluctuations is negligible.

Figure 11. The reduced distribution function $f(x,v_x)$ at time $t=70 \varOmega _u^{-1}$ from the beginning of the run with $v_{{\rm in}}=2v_A$. The coordinates are measured in the ion inertial lengths $c/\omega _{{\rm pi}}$, the velocity is measured in Alfvén velocities $v_A$ and the magnetic field is normalized on the upstream magnetic field $B_u$.

Figure 12. The main component of the magnetic field $B_z(x,z)$ (a) and the normal component of the magnetic field $B_x(x,z)$ (b) at time $t=70 \varOmega _u^{-1}$ from the beginning of the run.

Figure 13 shows the reduced distribution function $f(x,v_x)$ at time $t=70 \varOmega _u^{-1}$ from the beginning of the run with the upstream plasma speed $v_{{\rm in}}=3v_A$ in the frame of the simulation. The shock speed in this frame is $v_{{\rm sh}}\approx 0.8v_A$, so the Mach number is $M_A\approx 3.8$. The downstream gyration of the directly transmitted ions is very clear, together with the gradual relaxation. The overshoot is stronger, and there are reflected ions contributing to the downstream pressure. Figure 14 shows two-dimensional patterns of $B_z(x,z)$ (figure 14a) and $B_x(x,z)$ (figure 14b) at the same $t=70 \varOmega _u^{-1}$. The main component $B_z$ does not show any signs of inhomogeneity along the shock front. There are weak but noticeable fluctuations of the normal component of the magnetic field. Such fluctuations have been observed in low-Mach-number shocks, which otherwise were planar and stationary, according to all other signatures (Gedalin et al. 2022b). Figure 14 means that rippling starts to develop but does not affect the ion motion yet.

Figure 13. The reduced distribution function $f(x,v_x)$ at time $t=70\ \varOmega _u^{-1}$ from the beginning of the run with $v_{{\rm in}}=3v_A$.

Figure 14. The main component of the magnetic field $B_z(x,z)$ (a) and the normal component of the magnetic field $B_x(x,z)$ (b) at time $t=70 \varOmega _u^{-1}$ from the beginning of the run.

Figure 15 shows the reduced distribution function $f(x,v_x)$ at time $t=70 \varOmega _u^{-1}$ from the beginning of the run with the upstream plasma speed $v_{{\rm in}}=3.5v_A$ in the frame of the simulation. The shock speed in this frame is $v_{{\rm sh}}\approx 0.9v_A$, so the Mach number is $M_A\approx 4.4$. There is a strong overshoot in the profile, ensuring strong ion reflection. The reflected ions contribute substantially to the downstream ion pressure. Figure 16 shows the two-dimensional patterns of $B_z(x,z)$ (figure 16a) and $B_x(x,z)$ (figure 16b) at the same $t=70 \varOmega _u^{-1}$. Rippling is now clearly seen both in $B_z$ and in the large fluctuations of $B_x$. Thus, for the chosen parameters of the upstream state, the phase transition from a planar stationary profile to a rippled profile occurs in the vicinity of Mach number $M_A\approx 4$. The performed simulations confirm the theoretical predictions done with the use of the test particle analysis. More sophisticated simulations with a larger simulation box, more particles per cell and a smaller grid spacing would probably refine the details of the shock profile, ion distributions and the Mach number at which the transition occurs but would not change the physical picture drastically. We leave such simulations for future work.

Figure 15. The reduced distribution function $f(x,v_x)$ at time $t=70 \varOmega _u^{-1}$ from the beginning of the run with $v_{{\rm in}}=3.5v_A$.

Figure 16. The main component of the magnetic field $B_z(x,z)$ (a) and the normal component of the magnetic field $B_x(x,z)$ (b) at time $t=70 \varOmega _u^{-1}$ from the beginning of the run.

7 Conclusions

We used test particle analysis in a model shock profile, keeping the values of $\theta _{{\rm Bn}}$ and $\beta$ constant and Mach number increasing from $M_A=2$ to $M_A=4$. We have shown that momentum conservation (pressure balance) can be ensured if the shock is planar, stationary and non-structured, for $M_A=2$ and $M_A=2.5$. At higher Mach numbers, $M_A=3$ and $M_A=3.5$, the shock can still be planar and stationary but switches to a structured ‘phase’, in which the overshoot and undershoot play an important role in shaping the ion distributions in such a way so as to ensure pressure balance. For $M_A=4$, no planar stationary structure was found that could be consistent with the momentum conservation. We argue that between $M_A=3.5$ and $M_A=4$, the shock should undergo a ‘phase transition’, and the shock front should become rippled. The conclusion is supported by simple two-dimensional hybrid simulations.

Our usage of the term ‘phase transitions’ should not be confused with true thermodynamic phase transitions. The system under study is not in thermal equilibrium, and there is no analogy of the coexistence of phases, latent heat or order parameter. On the other hand, observations show that there are various types of shock structure, also known as ‘phases’. Our results indicate that the shock parameter space is divided into ranges, in each of which shocks can exist in only one of these ‘phases’, that is, have a certain type of structure. We have seen that, given the shock angle and the upstream $\beta$, the ‘phase’ in which the shock can exist changes when the Mach number increases. Such a ‘phase transition’ is not expected to be abrupt when the Mach number exceeds some threshold but may occur rapidly within a narrow range of Mach numbers.