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 (Reference Gedalin2016). 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. Reference Balikhin, Zhang, Gedalin, Ganushkina and Pope2008; Gedalin, Friedman & Balikhin Reference Gedalin, Friedman and Balikhin2015; Gedalin Reference Gedalin2016, Reference Gedalin2019a,Reference Gedalinb; Gedalin, Dröge & Kartavykh Reference Gedalin, Dröge and Kartavykh2016; Dimmock et al. Reference Dimmock, Gedalin, Lalti, Trotta, Khotyaintsev, Graham, Johlander, Vainio, Blanco-Cano, Kajdič, Owen and Wimmer-Schweingruber2023; Gedalin, Pogorelov & Roytershteyn Reference Gedalin, Pogorelov and Roytershteyn2023b). The ion velocities are acquired at a dense grid, using the staying time method (see technical details in Gedalin (Reference Gedalin2016) and Gedalin et al. (Reference Gedalin, Pogorelov and Roytershteyn2023b)), 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 Reference Gedalin2015), and to separate the effect of macroscopic fields from the effect of fluctuations (Gedalin et al. Reference Gedalin, Dimmock, Russell, Pogorelov and Roytershteyn2023a). 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. Reference Gedalin, Zhou, Russell, Drozdov and Liu2018; Pope, Gedalin & Balikhin Reference Pope, Gedalin and Balikhin2019). The method has been used for analysing whether a supercritical shock can be stationary (Gedalin Reference Gedalin2019a). 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 Reference Gedalin2016):

(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. Reference Greenstadt, Scarf, Russell, Formisano and Neugebauer1975, Reference Greenstadt, Russell, Gosling, Bame, Paschmann, Parks, Anderson, Scarf, Anderson, Gurnett, Lin, Lin and Rème1980; Scudder et al. Reference Scudder, Aggson, Aggson, Mangeney, Lacombe and Harvey1986; Farris et al. Reference Farris, Russell and Thomsen1993; Dimmock et al. Reference Dimmock, Balikhin, Krasnoselskikh, Walker, Bale and Hobara2012). The approximate relations between the field components were analytically derived (Jones & Ellison Reference Jones and Ellison1987; Gosling, Winske & Thomsen Reference Gosling, Winske and Thomsen1988; Schwartz et al. Reference Schwartz, Thomsen, Bame and Stansberry1988; Gedalin Reference Gedalin1996) and successfully used for comparison with observations (Gedalin et al. Reference Gedalin, Golbraikh, Russell and Dimmock2022b). Note that ions are not sensitive to fine details of the shock front or small-scale fields (Gedalin Reference Gedalin2020).

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 Reference Abraham-Shrauner1967; Hudson Reference Hudson1970; Chao & Goldstein Reference Chao and Goldstein1972; Sanderson Reference Sanderson1976; Lyu & Kan Reference Lyu and Kan1986; Kennel Reference Kennel1988; Vogl et al. Reference Vogl, Biernat, Erkaev, Farrugia and Mühlbachler2001; Livadiotis Reference Livadiotis2019; Gedalin et al. Reference Gedalin, Golan, Pogorelov and Roytershteyn2022a; Haggerty, Bret & Caprioli Reference Haggerty, Bret and Caprioli2022). 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 Reference Gedalin2016).

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. Reference Gedalin, Golan, Pogorelov and Roytershteyn2022a). 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. Reference Balikhin, Zhang, Gedalin, Ganushkina and Pope2008; Ofman et al. Reference Ofman, Balikhin, Russell and Gedalin2009; Pope et al. Reference Pope, Gedalin and Balikhin2019), 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 Reference Gedalin and Ganushkina2022; Gedalin et al. Reference Gedalin, Dimmock, Russell, Pogorelov and Roytershteyn2023a,Reference Gedalin, Pogorelov and Roytershteynb)

(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 Reference Gedalin and Sharma2023).

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.

6 Simulations

For illustrative purposes, we have performed several low-cost two-dimensional hybrid simulations using the dHybrid code (Gargaté et al. Reference Gargaté, Bingham, Fonseca and Silva2007). 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 Reference Gedalin and Ganushkina2022). 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. Reference Gedalin, Golbraikh, Russell and Dimmock2022b). 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.