Implications of weak rippling of the shock ramp on the pattern of the electromagnetic field and ion distributions

Abstract

Collisionless shocks undergo structural changes with the increase of Mach number. Observations and numerical simulations indicate development of time-dependent rippling. It is not known at present what causes the rippling. However, effects of such rippling on the field pattern and ion motion and distributions can be studied without precise knowledge of the causes and detailed shape. It is shown that deviations of the normal component of the magnetic field from the constant value indicate certain spatial dependence of the rippling. Deviations of the motional electric field from the constant value indicate time dependence. It is argued that whistler waves should propagate towards upstream and downstream regions from the rippled ramp. It is shown that the downstream pattern of the fields and ion distributions should follow the rippling pattern, while collisionless relaxation should be faster than in the stationary planar case.

1 Introduction

Collisionless shocks are among the most fundamental strongly nonlinear phenomena in space plasmas. The unfading interest regarding collisionless shocks and the ongoing research during the last six decades is determined by the fact the these shocks are the most efficient accelerators of charged particles in plasmas (Axford, Leer & Skadron 1977; Krymskii 1977; Bell 1978; Blandford & Ostriker 1978; Vasilev, Toptygin & Chirkov 1978; Toptyghin 1980; Jokipii 1982; Drury 1983; Blandford & Eichler 1987). The highest energies are achieved at supernova remnant (SNR) shocks (Vink 2020). Acceleration processes depend on the shock structure, which makes understanding of the latter the key issue in shock physics. All information about the SNR shocks is obtained remotely via the electromagnetic emission coming from the particles which are heated and accelerated by the shocks. In situ observations of collisionless shocks are possible only in the heliosphere. Supernova remnant shocks are believed to be high-Mach number shocks (Reynolds 2004; Vink 2004b, a; Jones 2011; Raymond 2018). Heliospheric shocks with Mach numbers possibly approaching the Mach numbers of SNR shocks are observed at outer planets (Masters et al. 2013; Sulaiman et al. 2015; Madanian et al. 2021). Most shock observations have been and are being performed at the Earth bow shock. As a result of these observations and theory development, it seems that at present the structure of low-Mach-number shocks is thoroughly studied observationally and understood rather well (Greenstadt et al. 1980; Russell et al. 1982; Mellott & Greenstadt 1984; Jones & Ellison 1987; Gosling, Winske & Thomsen 1988; Farris, Russell & Thomsen 1993; Gedalin 1996b; Balikhin et al. 2008; Gedalin, Friedman & Balikhin 2015; Gedalin et al. 2022) (see, however, Wilson et al. (2017)). With the increase of Mach number the shock front undergoes structural changes, developing rippling and time dependence (Bale et al. 2005; Moullard et al. 2006; Lobzin et al. 2008; Krasnoselskikh et al. 2013; Burgess et al. 2016; Hao et al. 2016). These changes do not appear abruptly when the Mach number exceeds some critical value but gradually become more and more pronounced when the Mach number increases (Ofman & Gedalin 2013). Rippling and time dependence (reformation), as well as generation of the propagating upstream whistlers (Wilson et al. 2009; Hull et al. 2012; Ramírez Vélez et al. 2012; Wilson et al. 2012,  2017) may be interrelated (Burgess et al. 2016; Gingell et al. 2017; Umeda & Daicho 2018; Omidi et al. 2021). At present, it is not quite clear what causes rippling. One of the plausible explanations is an instability of the waves propagating along the shock surface (Lowe & Burgess 2003; Burgess & Scholer 2007; Burgess et al. 2016). Such surface modulations should affect the processes at the shock front, among them ion reflection (Johlander et al. 2016,  2018). In this paper we study the implications of weak time-dependent rippling on the shock ramp and the adjacent upstream and downstream regions. We propose an analytical model of the magnetic and electric fields inside a weakly rippled shock transition layer, and analyse the consequences which can be verified with numerical simulations and applied to observations to estimate the rippling parameters.

2 Weak non-stationary rippling: a model

In the absence of a good theory we model rippling as spatial and temporal dependence localized within the ramp. We start with a monotonic magnetic profile of a low-Mach-number shock ramp which we model using the following expressions (Gedalin et al. 2015):

(2.1)\begin{gather} \frac{B_z}{B_u\sin\theta}=\frac{R+1}{2}+\frac{R-1}{2}\tanh\frac{3x}{D}, \end{gather}

(2.2)\begin{gather}B_x=B_u\cos\theta, \end{gather}

(2.3)\begin{gather}B_y=k_B\frac{{\rm d}B_z}{{\rm d} x}, \end{gather}

(2.4a,b)\begin{gather}E_y=V_uB_u\sin\theta, \quad E_z=0, \end{gather}

(2.5)\begin{gather}E_x={-}\frac{{\rm d}\phi_{{\rm NIF}}}{{\rm d} x}={-}k_E \frac{{\rm d}B_z}{{\rm d} x}. \end{gather}

Here, subscript $u$ refers to the upstream region, $\theta$ is the angle between the shock normal and the upstream magnetic field vector, $B_u$ is the upstream magnetic field magnitude, $E_y$ is the motional electric field, $E_x$ is the cross-shock electric field, $B_y$ is the non-coplanar component of the magnetic field, $D$ is the ramp width and $R$ is the ratio of the downstream to upstream $B_z$. The shock normal is along the $x$-direction and the non-coplanarity direction is $y$. The analysis is done in the normal incidence frame (NIF), where the upstream plasma flow is along the shock normal. The coefficients $k_E,k_B$ are obtained from

(2.6)\begin{gather} -\int E_x\,{\rm d}x=\phi_{{\rm NIF}}\equiv s_{{\rm NIF}}(m_pV_u^2/2e), \end{gather}

(2.7)\begin{gather}-\int (E_x+V_u\tan\theta B_y)\,{{\rm d} x}=\phi_{HT}\equiv s_{HT}(m_pV_u^2/2e), \end{gather}

where $\phi _{{\rm NIF}}$ and $\phi _{HT}$ are the cross-shock potentials in the NIF and the de Hoffman–Teller frame (HT), respectively. In the HT frame the upstream plasma flow is along the upstream magnetic field. The relations for $B_y$ and $E_x$ are approximations derived from the two-fluid plasma model for low-Mach-number shocks, both subcritical and supercritical (Goodrich & Scudder 1984; Gosling et al. 1988; Gedalin 1996b). Within this approximation $k_B=c\cos \theta /M\omega _{pi}$ is the whistler length. Here, $M=V_u/V_A$ is the Alfvén Mach number, $c/\omega _{pi}$ is the ion inertial length, $\omega _p=\sqrt {4{\rm \pi} n_ue^2/m_p}$, $V_A=B_u/\sqrt {4{\rm \pi} n_um_p}$ is the Alfvén speed, $n_u$ is the upstream proton number density and $m_p$ is the proton mass. For simplicity, the plasma is assumed to consist of protons and electrons only. For brevity, in what follows we write down expressions for the fields using $c\equiv 1$. The speed of light may be easily restored at the end using dimension arguments. Let us introduce the vector and scalar potentials, as follows:

(2.8)\begin{gather} A_y=A(x)+B_u\cos\theta z-V_uB_u\sin\theta t, \end{gather}

(2.9)\begin{gather}A_z=k_B\frac{{\rm d}A}{{\rm d} x}, \end{gather}

(2.10)\begin{gather}\phi=k_E\frac{{\rm d}A}{{\rm d} x}, \end{gather}

so that

(2.11)\begin{gather} B_z=\frac{\partial A_y}{\partial x}, \end{gather}

(2.12)\begin{gather}B_x=\frac{\partial A_y}{\partial z}=B_u\cos\theta, \end{gather}

(2.13)\begin{gather}B_y={-}\frac{\partial A_z}{\partial x}, \end{gather}

(2.14)\begin{gather}E_x={-}\frac{\partial\phi}{\partial x}, \end{gather}

(2.15)\begin{gather}E_y={-}\frac{\partial A_y}{\partial t}=V_uB_u\sin\theta, \end{gather}

(2.16)\begin{gather}E_z={-}\frac{\partial A_z}{\partial t}=0. \end{gather}

Let us now introduce rippling, as follows. Let $X=x+f$, $f(x,y,z,t)=a\psi (y,z,t)g(x)$, $g(x\rightarrow \pm \infty )=0$, $({\rm d}g/{{\rm d} x})(x\rightarrow \pm \infty )=0$, where $a$ is the amplitude (dimensions of length), while $\psi$ and $g$ are dimensionless. Consider a vector potential and a scalar potential

(2.17)\begin{gather} A_y=A(X)+B_u\cos\theta z-V_uB_u\sin\theta t, \end{gather}

(2.18)\begin{gather}A_z=k_B\frac{\partial A}{\partial X}, \end{gather}

(2.19)\begin{gather}\phi=k_E\frac{\partial A}{\partial X}. \end{gather}

Note that it is always possible to choose the gauge where one of $\boldsymbol {A}$ components vanishes. The fields are now

(2.20)\begin{gather} B_z=B_u\sin\theta B(1+f_x), \end{gather}

(2.21)\begin{gather}B_x=B_u\sin\theta\left(\cot\theta -Bf_z+k_B B_Xf_y\right), \end{gather}

(2.22)\begin{gather}B_y={-}k_BB_u\sin\theta B_X (1+f_x), \end{gather}

(2.23)\begin{gather}E_x={-}k_EB_u\sin\theta B_X(1+f_x), \end{gather}

(2.24)\begin{gather}E_y=V_uB_u\sin\theta-B_u\sin\theta Bf_t-k_EB_u\sin\theta B_X f_y, \end{gather}

(2.25)\begin{gather}E_z={-}B_u\sin\theta k_BB_X f_t-k_E B_u\sin\theta B_Xf_z. \end{gather}

The details of the derivation and the definitions of $B_X$, $f_t$, $f_x$, $f_y$ and $f_z$ are given in the Appendix (A).

3 Implications for the fields in the ramp

In what follows we analyse what could be the signatures of the weak rippling in observations and simulations. In the first order on derivatives $f_t,f_x,f_y,f_z,f_{xx}$ one has

(3.1)\begin{gather} \frac{\partial B_z}{\partial x}=B_u\sin\theta \left(B_X(1+2f_x)+ B f_{xx}\right), \end{gather}

(3.2)\begin{gather}\frac{\partial B_z}{\partial y}=B_u\sin\theta B_Xf_y, \end{gather}

(3.3)\begin{gather}\frac{\partial B_x}{\partial y}=\frac{\partial B_x}{\partial z}=0, \end{gather}

(3.4)\begin{gather}\frac{\partial B_y}{\partial x}={-}k_BB_u\sin\theta\left( B_{XX} (1+2f_x)+ B_X f_{xx}\right), \end{gather}

(3.5)\begin{gather}\frac{\partial B_y}{\partial z}={-}k_BB_u\sin\theta B_{XX} f_z. \end{gather}

The most notable distinction from the stationary planar shock is that $B_x$ and $E_y$ are no longer constant throughout the ramp, and $E_z\ne 0$. If $k_B$ is also small, the major deviations are

(3.6)\begin{gather} \delta B_x\approx{-}B_u\sin\theta Bf_z, \end{gather}

(3.7)\begin{gather}\delta E_y\approx{-}B_u\sin\theta Bf_t-k_EB_u\sin\theta B_X f_y, \end{gather}

(3.8)\begin{gather}\delta E_z\approx{-}k_E B_u\sin\theta B_Xf_z. \end{gather}

As an example, consider ripples localized within the ramp and propagating along the shock front, of the form $\psi =\sin (k_yy+k_zz-\omega t)$. Then there is no phase difference between $\delta B_x$, $\delta E_y$ and $\delta E_z$, and

(3.9)\begin{gather} \delta E_y={-}\frac{\omega}{k_z} \delta B_x + \frac{k_z}{k_y}\delta E_z, \end{gather}

(3.10)\begin{gather}\frac{\delta E_z}{\delta B_x}=k_E\frac{B_X}{B}. \end{gather}

Figure 1 illustrates the effect of rippling on the components of the fields. Figure 1(a,b) provides a visual comparison of the main magnetic component $B_z$ for a stationary shock and its rippled counterpart. Figure 1(c,d) shows the field components for which the effect is especially pronounced. For this visualization the following parameters were used: $M=2.5$, $\theta =65^\circ$, $B_d/B_u=2.2$, $s_{{\rm NIF}}=0.5$, $s_{HT}=0.1$, $D=c/\omega _{pi}$, $k_zV_u/\varOmega _u=2{\rm \pi}$, $k_yV_u/\varOmega _u={\rm \pi} /2$, $\omega =\varOmega _u$, $a=0.3(c/\omega _{pi})$, where $\varOmega _u=eB_u/m_pc$ and $(B_d/B_u)^2=R^2\sin ^2\theta +\cos ^2\theta$. The localizing function is $g(x)=\cosh ^{-2}(x/D)$. The rippling amplitude, $a$ and the wavelength along the $z$ direction $2{\rm \pi} /k_z$ are taken to be similar to what was found numerically by Ofman & Gedalin (2013), albeit for slightly higher Mach numbers. The profiles are shown for $y=0$ and $t=0$. Note that the rippling amplitude is rather small but the effect is quite noticeable, especially in the components $B_x$ and $E_y$. Thus, even if rippling may be difficult to recognize by the main magnetic component or the magnetic field magnitude, the two mentioned components easily disclose non-planarity and/or time dependence. For the chosen model of rippling $B_x$ deviates from the constant value mainly because of the spatial dependence on $z$, while $E_y$ deviates from the constant value due to the temporal dependence and the spatial dependence on $y$. Figure 1(c,d) suggest that the rippling parameters can be estimated from observations by comparing the variations of the components of the fields. For example, $\delta E_y+({\omega }/{k_z}) \delta B_x$ and $\delta B_x$ seem to not overlap, which suggests that minimization of $\int (\delta E_y+ \lambda \delta B_x)\delta B_x \,{{\rm d} x}$, where $\lambda$ is a variable parameter, may provide an estimate of $\omega /k_z$.

Figure 1. A visual comparison of the main magnetic component (a) $B_z$ for a stationary planar shock, (b) $B_z$ for a rippled shock, (c) $B_x$ for a rippled shock and (d) $E_y$ for a rippled shock. Parameters are given in text.

There is an overshoot with $\max (|B|/B_u)=2.73$, and the magnetic field magnitude drops to below the upstream value, $\min (|B|/B_u)=0.76$. The normal component of the magnetic field varies in the range $-0.99\leq B_x/B_u<1.84$ while without rippling one has $B_x/B_u=0.42$. For all practical purposes the local normal can be defined as the direction of $\boldsymbol {\nabla }|\boldsymbol {B}|$ in the region where the magnitude of the gradient is maximum. For the chosen model rippled profile the deviations of this direction from $\hat {x}$ reach values $> 30^\circ$, as can be seen from figure 2. Thus, for a weakly rippled shock the observational determination of the shock normal using magnetic coplanarity or minimum variance might have a $\pm 30^\circ$ error depending on the spacecraft trajectory across the shock.

Figure 2. Close-up on the magnetic field magnitude illustrating changes of the local normal direction along the rippled shock front.

Shocks which exhibit rippling typically have overshoots. An overshoot is not included in (2.1) and in the visualization. The generality of the expressions, however, is limited only by the model (2.17)–(2.19). An overshoot can be easily incorporated in the model profile, for example, as in Gedalin, Pogorelov & Roytershteyn (2021).

4 Implications for upstream and downstream waves close to the ramp

The rippled ramp is acting as a boundary which is perturbed according to $\psi =\sin (k_yy+k_zz-\omega t)$. These boundary perturbations should generate waves propagating towards upstream and downstream. The only low-frequency electromagnetic wave which can propagate towards upstream is the whistler wave which has the Doppler shifted dispersion relation

(4.1)\begin{equation} \tilde{\omega}=\omega-k_xV_u=\frac{\sqrt{k_x^2+k_y^2+k_z^2}(k_x\cos\theta+k_z\sin\theta)}{\omega_{pi} }v_A,\end{equation}

where $\tilde {\omega }$ is the whistler frequency in the frame of the upstream plasma flow and $k_x<0$. The relation (4.1) determines $k_x$. Note that for $k_y=k_z=0$ and $\omega =0$ this relation reduces to the phase-standing whistler for a planar stationary shock. In this interpretation the time-dependent rippled shock surface is the source of the whistlers propagating into the upstream region. The inverse should be also true: if whistlers are constantly escaping from the shock front they should leave a corresponding imprint at the front itself. Thus, the rippling pattern at the shock front and the whistler pattern in the upstream region can be expected to be mutually consistent. For the upstream whistler $\cos \theta _{\boldsymbol {k},\hat {n}}=|k_x|/k$ and $\cos \theta _{\boldsymbol {k},\boldsymbol {B}_u}=(k_x\cos \theta +k_z\sin \theta )/k$, where $\theta _{\boldsymbol {k},\hat {n}}$ and $\theta _{\boldsymbol {k},\boldsymbol {B}_u}$ are the angles between the propagation direction of the whistler and the shock normal and the upstream magnetic field, respectively. For the parameters chosen for the above visualization $k_xV_u/\varOmega _u\approx -3.86$, $\theta _{\boldsymbol {k},\hat {n}}\approx 57^\circ$, $\theta _{\boldsymbol {k},\boldsymbol {B}_u}\approx 57^\circ$.

In the downstream region the dispersion relation for the whistler changes accordingly,

(4.2)\begin{equation} \omega-k_{d,x}V_{d,x} - k_zV_{d,z}=\frac{\sqrt{k_{d,x}^2+k_y^2+k_z^2}(k_{d,x}\cos\theta_d+k_z\sin\theta_d)}{\omega_{pi,d} }v_{A,d}, \end{equation}

where $V_{d,x}$ and $V_{d,z}$ are the components of the downstream flow velocity, $\theta _d$ is the angle between the shock normal and the downstream magnetic field, $\omega _{pi,d}$ is the ion plasma frequency calculated with the downstream ion density and $v_{A,d}$ is the Alfvén speed in the downstream region. Note that in the downstream region waves propagate from the ramp into the downstream and the corresponding $k_{d,x}>0$. In numerical simulations capable of resolving whistlers two sets of fronts would be observed diverging from the ramp (Yuan et al. 2009; Riquelme & Spitkovsky 2011). If whistler waves are not resolved properly, small pieces of such diverging fronts may be still observed. Diverging waves should remove energy from the ramp, thus providing an additional channel of the redistribution of the energy of the directed flow of the incident ions. The amplitude of the diverging waves and the amplitude of the rippling depend on the mechanism which causes rippling and further sustains it. This mechanism is not known at present and is the subject of intensive studies (Lowe & Burgess 2003; Burgess & Scholer 2007; Johlander et al. 2016,  2018; Omidi et al. 2021).

5 Two-fluid hydrodynamics within the rippled shock

The two-fluid approach was used in attempts to describe the shock front of a laminar (low Mach number, low $\beta$) oblique shock (see, e.g. Gedalin 1998). Although a shock-like profile was not obtained, some useful estimates of the scales were derived. Here we outline a semiquantitative extension of the two-fluid description with the above modelled rippling. The two-fluid model has to be adapted separately to the different conditions inside the ramp and downstream of the ramp. There are no changes in the upstream region in comparison with the standard description (Gedalin 1998). The ramp width is substantially smaller than the ion convective gyroradius in supercritical and even laminar subcritical shocks (Russell et al. 1982; Mellott & Greenstadt 1984; Farris et al. 1993; Newbury & Russell 1996; Bale et al. 2005; Hobara et al. 2010; Krasnoselskikh et al. 2013). Therefore, it is more appropriate to treat the ions kinetically. The collisionless Vlasov equation which simply states that the distribution function is constant along the particle trajectory, $f_i(\boldsymbol {r}_i,\boldsymbol {v}_i,t)=f_0(\boldsymbol {r}_0,\boldsymbol {v}_0,t_0)$, where $\boldsymbol {r}_i$ and $\boldsymbol {v}_i$ are the solutions of the equations of motion

(5.1a,b)\begin{equation} \frac{{\rm d}\boldsymbol{r}_i}{{\rm d}t}=\boldsymbol{v}_i, \quad \frac{{\rm d}\boldsymbol{v}_i}{{\rm d}t}=\frac{e}{m_p}\left(\boldsymbol{E}+\frac{\boldsymbol{v}_i}{c}\times \boldsymbol{B}\right), \end{equation}

with the initial conditions $\boldsymbol {r}_i(t=t_0)=\boldsymbol {r}_0$, $\boldsymbol {v}_i(t=t_0)=\boldsymbol {v}_0$. The equations of motion inside the ramp are not integrable even in the stationary planar case. It was shown that a reasonable solution can be derived in the following approximation: (a) the ramp is narrow; and (b) the ratio of the upstream thermal speed to the flow speed is small, $v_{iT}/V_u\ll 1$ (Gedalin 1997,  2021). Here $v_{iT}=\sqrt {T_u/m_p}$, $T_u$ being the temperature of the incident ion distribution. In this approximation

(5.2a,b)\begin{gather} V_{i,x}=\sqrt{V_u^2-2e\phi_{{\rm NIF}}/m_m}, \quad V_{i,y}=V_{i,z}=0, \end{gather}

(5.3)\begin{gather}n_i=\frac{n_uV_u}{\sqrt{V_u^2-2e\phi_{{\rm NIF}}/m_p}}, \end{gather}

where $n_i$ is the ion number density and $\boldsymbol {V}_{i}$ is the bulk (hydrodynamical) velocity of the ions. The main effect is the deceleration by the cross-shock electric field as given by (5.2a,b)–(5.3).

The electron velocity $\boldsymbol {V}_e$ can be obtained from

(5.4)\begin{equation} \frac{4{\rm \pi} en}{c} (\boldsymbol{v}-\boldsymbol{V}_e)=\boldsymbol{\nabla}\times \boldsymbol{B}, \end{equation}

where quasineutrality $n_i=n_e=n$ is assumed. Neglecting ion velocity in the $y$ and $z$ directions inside the ramp one has

(5.5)\begin{gather} V_{ex}=V_{i,x}-\frac{B_u\sin\theta}{4{\rm \pi} ne}\left(B_Xf_y+k_BB_{XX}f_z\right), \end{gather}

(5.6)\begin{gather}V_{ey}=\frac{B_u\sin\theta}{4{\rm \pi} ne}\left(B_X(1+2f_x)+Bf_{xx}\right), \end{gather}

(5.7)\begin{gather}V_{ez}=\frac{B_u\sin\theta}{4{\rm \pi} ne}\left(k_BB_{XX}(1+2f_x)+k_BB_Xf_{xx}\right), \end{gather}

where $V_x=\sqrt {V_u^2-2e\phi /m}$ and we restricted ourselves with the lowest and first order only. In the lowest order the electron velocity is

(5.8)\begin{gather} V^{(0)}_{ex}=V_{i,x}, \end{gather}

(5.9)\begin{gather}V^{(0)}_{ey}=\frac{cB_u\sin\theta}{4{\rm \pi} ne}B_X, \end{gather}

(5.10)\begin{gather}V^{(0)}_{ez}=\frac{cB_u\sin\theta}{4{\rm \pi} ne}k_BB_{XX}, \end{gather}

as in the stationary planar case. In the approximation of massless cold electrons one has $\boldsymbol {E}+\boldsymbol {V}_e\times \boldsymbol {B}/c=0$ which eventually gives

(5.11)\begin{gather} \frac{e}{m_p}\left(\frac{{\rm d}\phi_{{\rm NIF}}}{{\rm d}X}\right)\frac{n_uV_u}{\sqrt{V_u^2-2e\phi_{{\rm NIF}}/m_p}}=\frac{1}{8{\rm \pi} m_p}\frac{{\rm d}}{{\rm d}X}B^2, \end{gather}

(5.12)\begin{gather}V_u-V_{i,x}=\frac{B^2-B_u^2}{8{\rm \pi} n_uV_um_p}, \end{gather}

(5.13)\begin{gather}\frac{V_{i,x}}{V_u}=1-\frac{b^2-1}{2M^2}, \quad b=B/B_u, \end{gather}

(5.14)\begin{gather}1-\frac{n_u}{n_i}=\frac{b^2-1}{2M^2} \rightarrow 1+\frac{1}{2M^2}=\frac{V_{i,x}}{V_u}+\frac{b^2}{2M^2}. \end{gather}

This expressions are identical to the expressions obtained for the stationary planar shock with the only replacement $x\rightarrow X$ (Gedalin 2021). The above means that in the lowest-order approximation the ramp structure and the ion velocity remain the same as in the stationary planar case, only the position of the ramp edges depend on $y,z,t$. This has important implications for the ion motion and the ion distribution just behind the ramp.

The main correction to the electron velocity due to the rippling is that $V_{ex}\ne V_x$, which probably may be observable in spacecraft measurements. In a quasi-perpendicular shock, where $k_B$ is small, this velocity difference is directly related to the spatial variations in the $y$ direction. For simplicity, we write down the expression $\boldsymbol {E}+\boldsymbol {V}_e\times \boldsymbol {B}/c=0$ for $\cos \theta =0$ and $k_B=0$ as

(5.15)\begin{gather} -k_E B_X(1+f_x)+\frac{B_u}{4{\rm \pi} ne}B\left(B_X(1+3f_x)+Bf_{xx}\right)=0, \end{gather}

(5.16)\begin{gather}V_u- Bf_t-k_E B_X f_y-B\left(V_x-\frac{B_u}{4{\rm \pi} ne}\left(B_Xf_y+k_BB_{XX}f_z\right)\right)=0, \end{gather}

(5.17)\begin{gather}-k_E B_Xf_z -\frac{B_u}{4{\rm \pi} ne}BB_Xf_z=0. \end{gather}

In the lowest order we have

(5.18a,b)\begin{equation} -k_E +\frac{B_u}{4{\rm \pi} ne}B=0, \quad V_u-V_{i,x}B=0. \end{equation}

Both mean $n/B=\text {const.}$, as should occur in a perpendicular shock. The first-order corrections add the following constraints:

(5.19)\begin{gather} -k_E B_Xf_x+\frac{B_u}{4{\rm \pi} ne}B\left(3B_Xf_x)+Bf_{xx}\right)=0, \end{gather}

(5.20)\begin{gather}- Bf_t-k_E B_X f_y+B\frac{B_u}{4{\rm \pi} ne}\left(B_Xf_y+k_BB_{XX}f_z\right)=0. \end{gather}

For $f=g(x)\sin (k_yy+k_zz-\omega t)$ the equations take the form

(5.21)\begin{gather} -k_E B_Xg_x+\frac{B_u}{4{\rm \pi} ne}B\left(3B_Xg_x)+Bg_{xx}\right)=0, \end{gather}

(5.22)\begin{gather}\omega B-k_yk_E B_X +B\frac{B_u}{4{\rm \pi} ne}\left(k_yB_X+k_zk_BB_{XX}\right)=0. \end{gather}

Note that in (5.22) all variables depend only on $X$ and all derivatives are with respect to $X$. It is tempting to interpret (5.22) as a consistence condition for $k_y,k_z,k_E,k_B$ and $B=(X)$, while (5.21) may be interpreted as an equation for viable $g(x)$. However, at this stage using (5.21) and (5.22) to place restrictions on the model parameters is premature, given the number of approximations made to separate the first-order corrections. These two equations only show that there are no gross inconsistencies in the proposed model of the rippling.

Two-fluid hydrodynamics in the downstream region requires taking into account the non-gyrotropy of the ion distribution and its slow gyrotropization. It is more convenient to replace the two-fluid approach with the conservation laws together with the collisionless relaxation principles (Gedalin et al. 2015).

6 Implications for downstream collisionless relaxation

Upon crossing the ramp ions begin to gyrate. In a stationary planar shock the total downstream ion pressure $p_{ij}$ is a function of the distance from the ramp $L$. The total pressure includes the dynamic pressure $nmV_iV_j$ and the kinetic pressure $P_{ij}$, $p_{ij}=nmV_iV_j+P_{ij}$, where $V_i$ is the bulk flow velocity. As a result of the kinematic collisionless relaxation, the kinetic pressure $P_{ij}$ gradually gyrotropizes and further isotropizes, while the dynamic pressure tensor reduces to three components only: $nmV_x^2$; $nmV_z^2$ and $nmV_xV_z$. In a stationary planar shock the total downstream ion pressure depends on the distance $x_d$ from the downstream edge of the ramp, $p_{ij}=p_{ij}(x_d)$. In the approximation of a small amplitude rippling this dependence may be replaced with the dependence on $x_n=x_d+a\psi (y,z,t)$. In a more general way, the ion distribution function becomes dependent on $\psi$. Therefore, all other moments, such as the bulk velocity vector, are also functions of $\psi$. For $\psi =\sin (k_yy+k_zz-\omega t)$ this would mean a spatially periodic pattern propagating along the shock front. Such a pattern should be easily observed in simulations but is difficult to identify even with four-spacecraft measurements. One immediate implication of the new dependence is that weak spatial and/or temporal averaging results in smearing out the peak values of the moments. In particular, the minimum value of averaged $p_{xx}$ is larger than the minimum value of $p_{xx}$ in the case if the shock were stationary and planar. Accordingly, the maximum value of averaged $p_{xx}$ decreases. The conservation laws read

(6.1)\begin{gather} n_t+\sum_{j=x,y,z} (nV_j)_j=0, \end{gather}

(6.2)\begin{gather}(nmV_i)_t+\sum_{j=x,y,z}(p_{ij}+\varPi_{ij})_j=0, \end{gather}

(6.3)\begin{gather}\varPi_{ij}=\frac{1}{8{\rm \pi}} (B^2\delta_{ij}-2B_iB_j). \end{gather}

Averaging over $y,z,t$ we arrive at the conservation laws in the form

(6.4)\begin{gather} \langle nV_x\rangle=\text{const.}, \end{gather}

(6.5)\begin{gather}\langle p_{ix}\rangle+\langle\varPi_{ix}\rangle=\text{const.}, \end{gather}

where $\langle \cdots \rangle$ means averaging. Smearing out the peak values of the pressure means a reduction of the amplitude of the downstream magnetic field oscillations. In this way rippling enhances kinematic collisionless relaxation. Figure 3 illustrates the dispersion of ion trajectories caused by the rippling. For this purpose figure 3(b) shows $x$ versus $v_x$ for an ion moving from upstream with the velocity of the flow, $\boldsymbol {v}_{\text {initial}}=(V_u,0,0)$, in the shock without rippling with the above chosen profile. Figure 3(a) shows $x$ versus $v_x$ for ions with $\boldsymbol {v}_{\text {initial}}=(V_u,0,0)$ starting at randomly chosen $0< z< V_u/\varOmega _u$ and the same $x_{\text {initial}}$. For figure 3(a) the shock is rippled with the parameters mentioned above but $k_y=0$ and $\omega =0$. The blue dotted line shows the magnetic field magnitudes corresponding to the positions $x$ of all ions, independently of $y$ and $z$.

Figure 3. Trajectories, $x$ versus $v_x$, for ions entering the shock with the velocity of the upstream flow at different initial $z$: (a) the rippled shock; (b) no rippling. The blue dotted line shows the magnetic field magnitudes corresponding to the positions $x$ of all ions.

Since the incident ions cross the ramp in different positions because of rippling, it can be expected that the downstream heating parameters, averaged over $y$, $z$ and $t$, would differ from those which are achieved in a stationary planar shock with the same parameters. In order to compare the parallel and perpendicular heating in a rippled shock with its stationary counterpart, 40 000 initially Maxwellian distributed ions with $\beta _i=0.2$ were traced across the shock in both cases and parallel and perpendicular temperatures were calculated well downstream of the transition region. The normalized upstream temperature in both cases is $T_u/m_pV_u^2=0.016$. The distributions behind the shock are strongly anisotropic. Without rippling the temperatures are $T_{d,\parallel }/m_pV_u^2= 0.016$ and $T_{d,\perp }/m_pV_u^2= 0.178$. No parallel heating occurs. With rippling we obtained $T_{d,\parallel }/m_pV_u^2= 0.018$, $T_{d,\perp }/m_pV_u^2= 0.182$, which means weak parallel heating and lower anisotropy.

Corrections to the averaged conservation laws may be obtained by replacing

(6.6)\begin{equation} \partial_i=a\psi_i \partial_x, \quad i=t,y,z \end{equation}

which gives

(6.7)\begin{gather} \partial_x\left(nV_x+\sum_{i=t,y,z} a\psi_i nV_i\right)=0, \end{gather}

(6.8)\begin{gather}\partial_x\left[a\psi_t nmV_i +p_{xx}+\varPi_{xx} + \sum_{j=y,z} a\psi_j(p_{ij}+\varPi_{ij})\right]=0. \end{gather}

Thus, for example, $nV_x$ and $p_{xx}+\varPi _{xx}$ are no longer constant throughout the shock but the deviations are proportional to the small amplitude of the rippling.

The above analysis was done in the lowest-order approximation which describes the rippling as a simple shift of the position of the ramp of the kind $\Delta x=a\psi (y,z,t)$. Beyond this approximation, the maximum magnetic field and the cross-ramp potential also change with the shift, which further affects the ion motion and can be expected to cause more efficient relaxation. As a result, the spatial scale of the gyrotropization and isotropization downstream of the ramp should be smaller than the corresponding scale in a stationary planar shock. In the present study we model rippling with a monochromatic wave. Even multispacecraft observations provide information only about a small part of the shock surface. Simulations (Umeda & Daicho 2018; Omidi et al. 2021) show that rippling is only approximately monochromatic. Finite width of the spectrum would further enhance gyrotropization and isotropization by adding randomness in the gyrophases of the ion which are mixed at a fixed spatial position behind the ramp.

7 Discussion

At present, there is no theory of the rippling development, and we do not know what the rippling parameters should be and how they depend on the shock parameters, such as Mach number and the angle between the shock normal and the upstream magnetic field. Therefore, our study focused on the implications of a time-dependent rippling for the pattern of the electric and magnetic fields inside the ramp and in the upstream and downstream regions adjacent to the ramp, as well as for the ion motion and distributions in the rippled structure. It is clear that the functional dependence of the vector potential on the coordinates and time may be different from the model adopted in this study. Accordingly, the field profiles may differ in detail from those shown in the figures. However, the general conclusions about the deviations from the stationary planar fields are rather independent of the particular shape. Upstream and downstream whistlers diverging from the ramp may not propagate far from the shock, and possibly only traces of such whistlers with the footprints at the rippled ramp would be observed. The downstream pattern should be induced by the ion distributions following the rippling pattern, and the relaxation to thermal equilibrium should be faster because of the enhanced mixing. Waves which cross the shock and instabilities in the foreshock, foot and downstream region would affect the ion motion and also change the observed fields. Yet, we expect that the main effect would be produced by the macroscopic fields of the rippled shock front. All numerical illustrations in this study were done with a rather modest rippling and for shock parameters which essentially excluded ion reflection. In supercritical shocks, reflected ions play an important role in the formation of downstream distributions and contribute significantly to ion heating. Ion reflection will be also affected by the rippling. The effect may be expected to be stronger since a reflected ion crosses the ramp up to three times (Gedalin 1996a,  2016; Balikhin & Gedalin 2022). Detailed study of the ion motion and distributions in a supercritical rippled shock, together with the dependence on the rippling parameters, will be performed later.

8 Conclusions

We modelled rippling of the shock front as a monochromatic wave propagating along the shock front but localized within the ramp. Such rippling causes similar patterns in the fields inside and around the ramp. The most prominent observable changes of the fields inside the ramp are variations of the normal component of the magnetic field and of the motional component of the electric field. The rippling should cause whistlers diverging from the ramp. The whistlers should propagate obliquely to the shock normal and to the ambient magnetic field. The downstream magnetic field pattern should roughly follow the pattern of the rippling. The amplitude of the magnetic field oscillations should be smaller than in the stationary planar shock with the same parameters. Gyrotropization and isotropization should occur at smaller scales than what would happen in a stationary planar shock with the same shock parameters. Our conclusions can be verified with numerical simulations and used to determine the rippling parameters from observations. The objective of this study was to show that even weak rippling has clear effects on the observable fields and particle distributions in the shock front. We defer for future studies analysis of more realistic rippling in supercritical shocks with significant ion reflection.