3.1. Non-structured profile

We start with the attempt to model the observed shock with a laminar profile with the calculated magnetic compression of $B_d/B_u=2.9$. Figure 2 presents the magnetic compression $B_d/B_u$ as a function of the Alfvénic Mach number $M_A$, obtained by solving the Rankine–Hugoniot relations (Kennel et al. Reference Kennel, Edmiston and Hada1985) for various upstream $\beta _u$. For a given $B_d/B_u$, higher $\beta _u$ requires higher $M_A$. It is seen that the measured $B_d/B_u$ and $\beta _u$ are inconsistent with the derived Mach number. Since the Rankine–Hugoniot relations are nothing but conservation laws and must be fulfilled, we will choose $B_d/B_u=2.8$, $\beta _u=12$ and $M_A=10$ for the modelling.

Figure 2. The magnetic compression $B_d/B_u$ as a function of the Alfvénic Mach number $M_A$, for various $\beta _u$, as obtained from the solution of the Rankine–Hugoniot relations (Kennel, Edmiston & Hada Reference Kennel, Edmiston and Hada1985).

We start with a monotonic magnetic field given by the expressions

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

(3.2)\begin{gather}B_x=B_u\cos\theta_{Bn} \end{gather}

(3.3)\begin{gather}B_y=\frac{c\cos\theta_{Bn}}{M_A\omega_{pi}} \frac{{\rm d} B_z}{{\rm d}\kern0.07em x} \end{gather}

(3.4)\begin{gather}\frac{B_d}{B_u}=\sqrt{R^2\sin^2\theta_{Bn}+\cos^2\theta_{Bn}}. \end{gather}

The expression (3.3) has been derived within a two-fluid study of stationary nonlinear waves (Gedalin Reference Gedalin1998). These expressions should be completed with the electric field, which takes the form: $E_z=0$, $E_y=V_uB_u\sin \theta _{Bn}/c$ and $E_x=-K_E B_y$. The coefficient $K_E$ is determined by the value of the cross-shock potential

(3.5)\begin{equation} e\phi={-}\int E_x\,{\rm d}\kern0.07em 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 _{pi}}$.

Figure 3 presents the reduced distribution function $f(x,v_x)=\int f(x,\boldsymbol {v})\,{\rm d} v_y\,{\rm d} v_z$. Heating occurs mainly due to the downstream gyration of the directly transmitted ions. Reflected ions contribute a small part. The kinematic collisionless relaxation due to the gyrophase mixing is seen very clearly. The decisive step is the calculation of the magnetic field from the conservation law

(3.6)\begin{equation} p_{xx}+ \left(\frac{n}{n_u}\right)^{5/3}(n_uT_u)+\frac{B^2}{8{\rm \pi}}=n_um_pV_u^2+2n_uT_u+ \frac{B_u^2}{8{\rm \pi}}, \end{equation}

where $n=\int f(x,v_x)\,{\rm d} v_x$. This magnetic field is shown in figure 4. It is seen that the calculated magnetic field approaches the model magnetic field well behind the shock transition. It should be mentioned that the result is sensitive to the value of $s$. At first sight, the model is successful. However, in a planar stationary shock the equation (3.6) should be valid throughout the shock, for each position $x$. However, figure 4 shows that (3.6) requires that the shock profile be structured. In particular, the shock structure should include a large overshoot. We, therefore, move on to modelling a structured planar stationary shock. Such shocks were shown to exist (Scudder et al. Reference Scudder, Aggson, Aggson, Mangeney, Lacombe and Harvey1986).

Figure 3. The reduced distribution function $f(x,v_x)=\int f(x,\boldsymbol {v})\,{\rm d} v_y\,{\rm d} v_z$ obtained by tracing 80 000 ions through the shock profile given by (3.1)–(3.3) and $s=0.32$.

Figure 4. The magnetic field calculated from (3.6).

3.2. A structured planar stationary shock model

Our objective is to model a structured shock profile, which would include at least a foot and an overshoot and would be similar to the observed shock profile, but not necessarily reproduce the latter. Such modelling requires first denoising the magnetic profile, that is, removing oscillations that are supposed not to be a part of the shock structure but, possibly, transient waves or observed due to the tangential crossing of the shock by the spacecraft. Figure 4 makes the impression that at least one magnetic maximum corresponds to the overshoot, after which the magnetic field drops to the downstream value or below it. The position of the overshoot corresponds to the position of the maximal deceleration of the directly transmitted ions in figure 3. Figure 3 also shows that the reflected ions appear as far as approximately $0.5(V_u/\varOmega _u)$ from the ramp. This distance should approximately correspond to the foot length of a structured shock. Figure 5 compares the original profile with a denoised one. The latter is obtained by applying Daubechies 10 wavelet transform, reducing the smaller scale 9 levels, and applying the inverse transform (WaveShrink procedure of the WaveLab 850 package).

Figure 5. The original magnetic field magnitude truncated to $2^{13}$ points (red) and the denoised magnetic field. The denoising is done by applying Daubechies 10 wavelet transform, reducing the smaller scale 9 levels, and applying the inverse transform.

An overshoot is added as follows:

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

(3.8)\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}

For the run below the following parameters are chosen: $a=1$, $w_l=0.3(V_u/\varOmega _u)$, $w_r=1.87(V_u/\varOmega _u)$, $x_l=0.1(V_u/\varOmega _u)$, $x_r=0.15(V_u/\varOmega _u)$. Because of the overshoot, the cross-shock potential has to be modified, so the downstream value $s_d=0.2$ while ${s_{\max }=0.35}$.

Figure 6 shows the reduced distribution function $f(x,v_x)$, in the format similar to figure 3. It has been shown that an overshoot enhances ion reflection (Gedalin et al. Reference Gedalin, Dimmock, Russell, Pogorelov and Roytershteyn2023a). This enhancement is easily seen comparing figures 3 and 6. Ion reflection is strong because of the large $v_T/v_u=0.17$ (Sharma & Gedalin Reference Sharma and Gedalin2023). The much larger fraction of the reflected ions in the structured shock reverses the behaviour of the total ion pressure $p_{xx}$ across the shock from decreasing, as required by the conservation laws, to increasing, which is impossible in a stable shock. Although figure 7 shows this behaviour for one set of parameters but it has been found that no variation of parameters can reduce the far downstream $p_{xx}$ to below the upstream value. Thus, the observed profile (figures 1 and 5) cannot belong to a planar stationary shock. Namely, rippling should be explored as the next viable possibility.

Figure 6. The reduced distribution function $f(x,v_x)=\int f(x,\boldsymbol {v})\,{\rm d} v_y\,{\rm d} v_z$ obtained by tracing 80 000 ions through the structured shock profile. The black line shows the magnetic field magnitude.

Figure 7. Ion $p_{xx}$ throughout the shock. Black line: the non-structured shock. Red line: the planar stationary structured shock.

3.3. A rippled shock: simulations

A self-consistent two-dimensional (2-D) hybrid kinetic simulation was performed in order to provide a reference configuration for the modelling effort. In the hybrid simulation model, the ions are treated kinetically, while electrons are modelled as a massless fluid with a prescribed equation of state (e.g. Winske et al. Reference Winske, Karimabadi, Le, Omidi, Roytershteyn and Stanier2023). The simulation was performed in a 2-D domain of size $L_x \times L_y = (1024 \times 256)(c/\omega _{pi})$, covered by a uniform grid with $8192 \times 2048$ cells. The upstream parameters are $\theta _{Bn}=115^\circ$ and $\beta _e= \beta _i = 6$, where the distribution function for the ions is a drifting Maxwellian. The upstream magnetic field is in the $x$–$y$ plane. Note that the $y$-coordinate in the simulation corresponds to the $z$-coordinate in the data and test particle analyses. The plasma is injected from $x=0$ of a 2-D simulation domain with the speed $V_0 = 7 V_A$. Reflecting boundary conditions are used at $x=L_x$, and the shock is formed by the interaction between the incoming and reflecting flows. In the simulation frame of reference, the shock propagates in the negative $x$ direction with the average speed $V_{sh} \approx 3.35 V_A$, resulting in total upstream flow in the NIF $V_u \approx 10.35 V_A$. The simulation was performed using a version of the H3D code (Karimabadi et al. Reference Karimabadi, Vu, Krauss-Varban and Omelchenko2006) adapted for shock simulations. An adiabatic equation of state with adiabatic index $\gamma =5/3$ was used for the electrons. Ion distribution is sampled by computational particles with a uniform and constant statistical weight, such that in the upstream region the number of particles per cell is $N_{{\rm ppc}} = 100$. The time step used in the simulation is $\delta t \varOmega _u = 1.25 \times 10^{-3}$. In the discussion below, coordinates are normalized to upstream ion inertial length $c/\omega _{pi}$, while time is normalized to $\varOmega _u$.

The basic properties of the shock are summarized in figure 8, which shows $x$-profiles of the magnetic field components $B_x(x',y=\text {const}, t\varOmega _u=125)/B_u$, $B_y(x',y=\text {const}, t\varOmega _u=125)/B_u$, normal component of the velocity $V_x'(x',y=\text {const}, t\varOmega _u=125)/V_u$ and density $\rho (x',y=\text {const}, t\varOmega _u=125)/\rho _u$, taken at various values of $y=\text {const}$ (left column). Here, $V_x'=V_x+V_{sh}$ is the velocity in NIF. There is substantial variance of the above variables. The right column shows the $y$-averaged profiles

(3.9)\begin{equation} \langle X\rangle_y(x',t)=\frac{1}{L_y}\int {{\rm d}y}\kern0.07em X(x',y,t), \end{equation}

where $X=B_x,B_y,V_x,\rho$, for various moments in the time interval $125< t\varOmega _u <220$, and $x' = x + V_u(t-t_0)$ is the NIF coordinate. The $y$-average profiles are stationary in NIF.

Figure 8. Profiles of (top to bottom) $B_x$, $B_y$, $\rho$ and $V_x$ in the hybrid simulation. The left column shows $y=\text {const}$ cuts taken at each cell, i.e. there are 2048 cuts, and a single instance of time $t\varOmega _{u}=125$. In each panel, each blue line is the single profile, e.g. $B_y(x, y=const, t=\text {const})$, while the black line corresponds to the $y$-average. The vertical dash lines mark the region used in the analysis of shock front perturbations and distribution function below. The right column shows $\langle B_x\rangle _y(x',t)$, $\langle B_y\rangle _y(x',t)$, $\langle V_x'\rangle _y(x',t)$ and $\langle \rho \rangle _y(x',t)$, for various moments in the time interval $125< t\varOmega _u <220$. Here, $x' = x + V_u(t-t_0)$ is the NIF coordinate, with $t_0\varOmega _u = 125$, and $V_x=V_x'+V_{sh}$.

The structure of the shock front perturbations is further illustrated in figure 9, which shows a section of the simulation domain of size $(60 \times 60) c/\omega _{pi}$ near the shock front. Several distinct types of perturbations are present, with different wavelengths and polarization properties. We identify the shock rippling by the enhancements at the shock front, shown in the figure by red arrows. The wavelength of these fluctuations is approximately $25 c/\omega _{pi}$. These perturbations propagate in the negative $y$-direction, with the velocity $V_y \approx -6 V_A$. Note that the fronts of the rippling wave are not perpendicular to the shock front but are oblique. One of these fronts is along a red arrow in the figure. Figure 10 illustrates the frequency spectrum of $B_y$ fluctuations along the shock front in the NIF. The spectrum is obtained by performing the fast Fourier transform of $B_y$ collected in a narrow region around $x' = 593.5 c/\omega _{pi}$ (see figure 8). The spectrum shows the existence of perturbations with wavenumbers in the range $k_y \lesssim 0.6$, propagating roughly with the same phase speed in the negative $y$ direction. Further, there appears to exist a second class of weaker fluctuations, propagating in the opposite direction with a somewhat smaller phase speed.

Figure 9. Two-dimensional profiles of (clockwise) $B_x$, $B_y$, $B_z$, $\rho$ and pressure tensor $\hat p_{xx}$ in a hybrid simulation at $t\varOmega _u=125$.

Figure 10. Frequency–wavenumber spectrum of $B_y$ fluctuations along the shock front.

The reduced ion distribution function collected in the same narrow region of $x$ is shown in figure 11. The same rippling as in figure 9 is seen in the reduced distribution function and is marked by red arrows. The rippling seems to be closely related to the non-locality of the reflection process is also seen (Gedalin Reference Gedalin2023).

Figure 11. Reduced distribution function at the shock front (left to right): $f(y,v_x)$, $f(y,v_y)$, $f(y,v_z)$. The rightmost panel shows $B_y(y)$ at the same $x$-location. The distributions are shown in the NIF.

Despite the spatial and time variability of the shock structure, the shock maintains stable average fluxes of the conserved quantities, as expected. This is shown in figure 12, which shows fluxes of mass $F_m$, total momentum $F_p$ and total energy $F_\epsilon$ in the simulation, normalized to their corresponding upstream values. Omitting resistivity terms, the conserved fluxes take the following form in the hybrid model:

(3.10)\begin{gather} F_m = n_i V_i, \end{gather}

(3.11)\begin{gather}F_p=\hat p_{i,xx} + m_i n_i V_{ix}^2 + p_e +\frac{1}{4{\rm \pi}} \left ( \frac{B^2}{2} - B_x^2\right ), \end{gather}

(3.12)\begin{gather} F_\epsilon = \chi_{i,x}+(\boldsymbol{V}_i\boldsymbol{\cdot}\hat{p}_i)_x+\frac{1}{2}m_i n_i V_{i}^2V_{ix} \nonumber\\ \quad +\frac{1}{2} \mathrm{Tr} (\hat p_i)V_{ix}+\frac{\gamma}{\gamma-1} p_eV_{ex}+\frac{c}{4{\rm \pi}} (\boldsymbol{E}\times\boldsymbol{B})_x, \end{gather}

where

(3.13)\begin{equation} \hat p_{i,ab}(x',y,t)=m_p\int (v_a-V_a)(v_b-V_b)f(\boldsymbol{v},x',y,t)\,{\rm d}^3\boldsymbol{v}, \end{equation}

is the ion kinetic pressure tensor in the co-moving frame, and

(3.14)\begin{equation} \chi_{i,x}= \frac{1}{2}m_p\int (\boldsymbol{v}-\boldsymbol{V})^2(v_x-V_x)f(\boldsymbol{v},x',y,t)\,{\rm d}^3\boldsymbol{v}, \end{equation}

is the ion heat flux. The fluxes in (3.10)–(3.12) are further $y$-averaged. We observe that time-averaged fluxes of mass and momentum are conserved to approximately 0.5 % or better between the upstream and downstream regions. The average energy flux shows somewhat greater variation, a little over 1 %, across the region included in figure 12. This is likely due to a combination of numerical dissipation in the simulation and slight time variability of the upstream conditions (which could be traced to numerical dissipation in the upstream region).

Figure 12. Fluxes of (top to bottom) mass, momentum and energy in the hybrid simulation. In each panel, the coloured lines show flux computed at an instance of time as indicated by the colour bar. The black lines show time-averaged fluxes. Each flux is normalized by its upstream value.

3.4. A rippled shock: modelling

We model a rippled shock front assuming that the shock is approximately stationary in the frame moving with the ripples along the shock front with the ripple velocity $V_u(V_{ry},V_{rz})$. We assume that the residual time dependence seen in the simulations is of secondary importance and may be ignored when modelling the shock structure. In what follows we exploit the principles proposed by Gedalin & Ganushkina (Reference Gedalin and Ganushkina2022). In the rippling frame, the upstream plasma velocity is $V_u(1,- V_{ry},-V_{rz})$, so that the electric field would be

(3.15)\begin{gather} E^{(0)}_x={-}V_{ry}\sin\theta_{Bn}-K_E\frac{{\rm d} B_z}{{\rm d}\kern0.07em x} \end{gather}

(3.16)\begin{gather}E^{(0)}_y=\sin\theta_{Bn}+V_{rz}\cos\theta_{Bn} \end{gather}

(3.17)\begin{gather}E^{(0)}_z={-}V_{ry}\cos\theta_{Bn}, \end{gather}

if the shock were planar, which is denoted by the superscript ${(0)}$. We introduce the scalar and vector potentials

(3.18a,b)\begin{gather} A^{(0)}_y,\quad \frac{\partial A^{(0)}_y}{\partial x}=B^{(0)}_z(x), \end{gather}

(3.19a–c)\begin{gather}A^{(0)}_z,\quad \frac{\partial A^{(0)}_z}{\partial y}=B^{(0)}_x,\quad -\frac{\partial A^{(0)}_z}{\partial x} =B^{(0)}_y, \end{gather}

(3.20)\begin{gather}\phi=V_{ry}\sin\theta_{Bn} x+K_EB^{(0)}_z-(\sin\theta_{Bn}+V_{rz}\cos\theta_{Bn})y+V_{ry}\cos\theta_{Bn} z, \end{gather}

where $\boldsymbol {B}^{(0)}$ is given by (3.1)–(3.3). Rippling is added by replacing the dependence on $x$ with the dependence on

(3.21)\begin{gather} X=x-\psi \end{gather}

(3.22)\begin{gather}\psi=a_r\cos(k_xx+k_yy+k_zz)G \end{gather}

(3.23)\begin{gather}G=g(x,x_L,w_L,x_R,w_R). \end{gather}

The shape (3.22) is a generalization of waveform proposed by Gedalin & Ganushkina (Reference Gedalin and Ganushkina2022) and is inspired by the above simulation showing that the front of the wave is not perpendicular to the shock front. Thus, rippling is a nonlinear wave propagating at an angle to the shock front. This wave damps into the downstream region and seems to not propagate into the upstream region. This behaviour is modelled using the localization function $g$. The replacement $x\rightarrow X$ is to be done in the vector potential and scalar potential. Then

(3.24)\begin{gather} B_x=\frac{\partial A_z}{\partial y} - \frac{\partial A_y}{\partial z}=\cos\theta -\psi_y B_y^{(0)}+\psi_zB^{(0)}_z \end{gather}

(3.25)\begin{gather}B_y=\frac{\partial A_x}{\partial z} - \frac{\partial A_z}{\partial x}=B_y^{(0)}(1-\psi_x) \end{gather}

(3.26)\begin{gather}B_z=\frac{\partial A_y}{\partial x} - \frac{\partial A_x}{\partial y}= B_z^{(0)}(1-\psi_x) \end{gather}

(3.27)\begin{gather}E_x={-}\frac{\partial \phi}{\partial x}=({-}V_{ry}\cos\theta_{Bn}-K_EB^{(0)}_y)(1-\psi_x) \end{gather}

(3.28)\begin{gather}E_y={-}\frac{\partial \phi}{\partial y}=\sin\theta_{Bn}+V_{rz}\cos\theta_{Bn} +V_{yr}\cos \theta_{Bn}\psi_y+K_EB_y^{(0)}\psi_y \end{gather}

(3.29)\begin{gather}E_z={-}\frac{\partial \phi}{\partial z}={-}V_{ry}\cos\theta_{Bn} +K_EB_y^{(0)}\psi_z, \end{gather}

where $B_z^{(0)}=B_z^{(0)}(X)$, $B_y^{(0)}=B_y^{(0)}(X)$ and $\psi _\xi =\partial \psi /\partial \xi$. For the present analysis, we add rippling to the non-structured shock profile (3.1). The parameters chosen to make the model resemble the simulated shock are

(3.30a–e)\begin{gather} k_y=0,\quad k_z=2{\rm \pi} /3,\quad k_x=k_z\cos(145^\circ),\quad V_{ry}=0,\quad V_{rz}={-}0.6 \end{gather}

(3.31a–d)\begin{gather}x_L=x_R=0,\quad w_L=1/M_A,\quad w_R=8,\quad a_r=0.75/M_A . \end{gather}

The three components of the magnetic of the rippled shock model are shown in figure 13 (central part only). Figure 14 further illustrates the structure of the magnetic field of the model, showing two families of the magnetic field magnitude cuts. The model was used to trace 160 000 ions. The incident distribution was the same Maxwellian as above. The initial $x$ and $y$ were the same for all ions, while the initial $z$ were randomly distributed along the wavelength $0\leq z<3$. Figure 15 presents the results of the tracing as a 2-D distribution of the ion pressure $p_{i,xx}$, where

(3.32)\begin{equation} p_{i,xx}(x,z)=m_p\int v_x^2 f(\boldsymbol{v},x,z)\,{\rm d}^3\boldsymbol{v}. \end{equation}

The distribution function $f(\boldsymbol {v},x,z)$ is obtained by counting particles appearing in $500\times 50$ cells of the size $\Delta x=0.018$ and $\Delta y=0.06$. During the tracing $z$-coordinate of an ion may leave the stripe $0\leq z<3$. In this case, it is counted with $z$-coordinate shifted toward the stripe by adding an integer number of wavelengths. Each ion count is multiplied with the weight $|v_{0x}|$ (Gedalin, Pogorelov & Roytershteyn Reference Gedalin, Pogorelov and Roytershteyn2023b). Figure 16 shows the reduced distribution function $f(x,v_x)=\int f(\boldsymbol {v},x,z) \,{\rm d} v_y\,{\rm d} v_z\,{\rm d} z$.

Figure 13. The three components of the magnetic of the rippled shock model. Three wavelengths along the shock front are shown.

Figure 14. (a) Value of $|B|/B_u$ as a function of $x$ for various $z$. (b) Value of $|B|/B_u$ as a function of $z$ for various $x$.

Figure 15. (a) Value of $p_{i,xx}$ as a function of $x$ and $z$. (b) Value of $B_x$ as a function of $x$ and $z$.

Figure 16. The reduced distribution function $f(x,v_x)$.

The behaviour of ions is further illustrated by the cuts at $x=0$. Figure 17 shows the distribution function $f(z,v_x)$ at $x=0$. Figure 18 shows $z$-integrated distributions $f(v_x,v_y)$ and $f(v_x,v_z)$ for $x=0$. Finally, figure 19 shows the magnetic field magnitude calculated using the planar stationary momentum balance (3.6). Only the values far upstream and downstream have physical sense since in the region, where rippling is still noticeable, (3.6) is not valid and should be replaced with a more complicated version (Gedalin & Ganushkina Reference Gedalin and Ganushkina2022). The downstream value of $|B|/B_u$ is rather close to $B_d/B_u$ chosen in the model.

Figure 17. The distribution function $f(z,v_x)$ at $x=0$.

Figure 18. (a) Value of $f(v_x,v_y)$ at $x=0$. (b) Value of $f(v_x,v_z)$ at $x=0$.

Figure 19. The magnetic field magnitude calculated using the planar stationary momentum balance (3.6). The derived downstream magnetic field is consistent with the chosen model magnetic field. The conversion to the downstream value starts rather quickly, less than three ion upstream convective gyroradii.