3.1. Wavelet-Fourier transform

Since instabilities typically are dominated by one or a few of the most rapidly growing unstable modes, our first task is to transform the fluctuating fields into local plane-wave modes using a combined WFT. Since the boundary conditions in the transverse plane $(y,z)$ are periodic, the fluctuating electric field $\delta \boldsymbol {E}(x,y,z)$ and fluctuating magnetic field $\delta \boldsymbol {B}(x,y,z)$ are Fourier transformed over these two directions to obtain the complex Fourier coefficients as a function of $(x,k_y,k_z)$, yielding the transformed fields $\widetilde {\delta \boldsymbol {E}} (x,k_y,k_z)$ and $\widetilde {\delta \boldsymbol {B}} (x,k_y,k_z)$.

Next, since the steady-state shock physics leads to significant non-periodic variations in the normal direction, we employ the wavelet transform (Torrence & Compo Reference Torrence and Compo1998) along the shock normal, $\hat {\boldsymbol {n}}=\hat {x}$,

(3.3)\begin{equation} \delta \hat{\boldsymbol{E}}(k_x,k_y,k_z;x) \equiv W_\psi \{\widetilde{\delta \boldsymbol{E}}(x,k_y,k_z)\} = \sqrt{|k_x|} \int_{-\infty}^{\infty} {{\rm d} x}^\prime \psi^*_{\sigma_0}[(x^\prime-x) k_x] \widetilde{\delta\boldsymbol{E}}(x',k_y,k_z), \end{equation}

with the complex Morlet wavelet

(3.4)\begin{equation} \psi_{\sigma_0}(x)={\rm e}^{{\rm i} k_x x} \,{\rm e}^{-({1}/{2})({x k_x}/{\sigma_0})^2} {\rm \pi}^{{1}/{4}}\sqrt{\frac{k_x}{\sigma_0}},\end{equation}

where $\sigma _0$ is a parameter that allows a trade off between resolution in position and wavenumber (Najmi & Sadowsky Reference Najmi and Sadowsky1997) and $\widetilde {\delta \boldsymbol {E}}(x,k_y,k_z)$ is the Fourier transformed fluctuating electric field. At higher $\sigma _0$, the wavenumber is better constrained at the cost of spatial resolution. It should be noted that we normalize the result at each wavenumber in order to maintain the relative energy of each mode, as suggested in Torrence & Compo (Reference Torrence and Compo1998). By combining the Fourier transform of the perturbed fields in the transverse plane with the wavelet transform in the normal direction, we are able to determine the approximate wavevector of unstable fluctuations locally at position $x$, enabling us to explore the properties of the unstable modes as they vary in the normal direction through the shock.

The same WFT is applied to obtain the local wavevectors of the perturbed magnetic field as a function of the position along the shock normal, $\delta \hat {\boldsymbol {B}}(k_x,k_y,k_z;x)$. Figure 5 shows that the unstable mode amplitudes peak in the shock transition region, rather than in the upstream or downstream regions. In this figure, we plot the wavelet transform of the dominant transverse Fourier mode $(k_y d_{i,0},k_z d_{i,0})=(-0.52,-1.05)$ at position $x_0/d_{i,0}=39.875$, as shown in figure 6(c). Here, $|\delta \hat {\boldsymbol {B}}(x_0=39.875d_{i,0})|$ assumes a maximum at $(k_x d_{i,0},k_y d_{i,0},k_z d_{i,0})= (1.62,-0.52,-1.05)$. There exists non-zero amplitude of $|\delta \hat {\boldsymbol {B}}(x_0=39.875d_{i,0})|$ far upstream at low $k_x$, although this may be due to the inability of the wavelet transform to localize waves with small wavenumber. Most critically, we note that the dominant transverse Fourier mode $(k_y d_{i,0},k_z d_{i,0})=(-0.52,-1.05)$ matches the observed ripple structure seen in figure 4. We will further motivate our selection of this wave mode in § 3.2, explain how to measure frequency using WFT in § 3.3 and compare our measured wave with dispersion relations from gyrotropic linear theory in § 3.4 to show that we are able to isolate and identify the wave present in this simulation.

Figure 5. (a) Wavelet transform of the magnetic field fluctuations, $\delta \hat {\boldsymbol {B}}$, of an oblique shock with a shock velocity of $M_A \approx 7.88$ and a shock normal of $\theta _{B_n} = 45^\circ$ with fixed $k_{y,0} = -0.52$ and $k_{z,0} = -1.05$ (i.e. $\delta \hat {\boldsymbol {B}}(x;k_x,k_{y,0},k_{z,0})$) at time $t=20 \varOmega _{i,0}^{-1}$. A vertical black line indicates the position $x/d_{i,0}=39.875$ at which we determine the dominant wave mode. We see larger values for $|\delta \hat {\boldsymbol {B}}|$ in the ramp and overshoot of the shock. (b) Compressible magnetic field component, $\overline {B}_z$, for reference.

Figure 6. Projections of $|\delta \hat {\boldsymbol{B}}(k_{x},k_{y},k_{z})|$ (a–c) and $|\delta \hat {\boldsymbol {B}}(k_{\|},k_{\perp 1},k_{\perp 2})|$ (d–f) in the shock ramp at $x_0/d_{i,0}=39.875$ using hexagonal binning due to the non-uniformity of data in FACs. There is a dominant wave mode in the FAC system with $(k_{\perp 1}d_{i,0},k_{\perp 2}d_{i,0},k_{\|}d_{i,0}) = (1.97,0.00,0.34)$, corresponding to $(k_{x}d_{i,0},k_{y}d_{i,0},k_{z}d_{i,0}) = (1.62,-0.52,-1.05)$ in the simulation coordinates, along with its conjugate mode due to the reality condition. There is little power off the $k_{\perp 2}$ axis. The ‘plus’ structure in panel (e) arises due to our projection of a grid in simulation coordinates $(x,y,z)$ onto our FAC system.

3.2. Magnetic FAC system

We construct a local magnetic FAC system at $x_0/d_{i,0} =39.875$, defined by the orthonormal basis $(\hat {e}_{\perp 1},\hat {e}_{\perp 2},\hat {e}_{\parallel })$, where

(3.5)\begin{gather} \hat{e}_{{\parallel}} = \frac{\overline{\boldsymbol{B}}(x_0)}{|\overline{\boldsymbol{B}}(x_0)|}, \end{gather}

(3.6)\begin{gather}\hat{e}_{{\perp} 2} = \frac{\overline{\boldsymbol{B}}(x_0) \times \hat{\boldsymbol{n}}} {|\overline{\boldsymbol{B}}(x_0) \times \hat{\boldsymbol{n}}|} =\frac{\overline{\boldsymbol{B}}(x_0) \times \hat{x}} {|\overline{\boldsymbol{B}}(x_0) \times \hat{x}|}, \end{gather}

(3.7)\begin{gather}\hat{e}_{{\perp} 1} = \frac{\hat{e}_{{\perp} 2} \times \hat{e}_{{\parallel}}} {|\hat{e}_{{\perp} 2} \times \hat{e}_{{\parallel}}|}. \end{gather}

At $x_0$, we measure the averaged magnetic field $\overline {\boldsymbol {B}}(x_0)/B_0 = (0.707, 0.698, 1.390)$, which notably has a substantial non-zero component in the $\hat {y}$ direction.

Figure 6 shows a singular dominant WFT wave mode (with conjugate pair) around $(k_x d_{i,0},k_y d_{i,0},k_z d_{i,0}) = (1.61,-0.52,-1.05)$, which corresponds to $(k_{\perp 1}d_{i,0},k_{\perp 2}d_{i,0}, k_{\|}d_{i,0}) = (1.97,0.00,0.34)$ in our FAC system. The localization of amplitude around $k_{\perp 2}=0$ suggests that our FAC system is a natural coordinate system for describing the unstable mode underlying the shock ripple. Note that the transverse Fourier wavevector of this dominant unstable mode $(k_y d_{i,0},k_z d_{i,0})=(-0.52,-1.05)$ is consistent with the dominant plane-wave structure seen in the total magnetic field seen in figure 4.

3.3. Frequency estimation using Faraday's law

Since the dominant wave mode in the FAC system has $k_{\perp 2}=0$, we take the wavevector to lie in the $(\hat {e}_{\perp 1},\hat {e}_{\parallel })$ plane, such that $\boldsymbol {k} = k_{\perp 1} \hat {e}_{\perp 1} + k_\parallel \hat {e}_{\parallel }$. Assuming local plane-wave modes that vary as $\exp ({\rm i} \boldsymbol {k} \boldsymbol{\cdot} \boldsymbol {x} - {\rm i} \omega t)$, we Fourier transform Faraday's law in time and space to relate $\omega \boldsymbol {B} = \boldsymbol {k}\times \boldsymbol {E}$. We convert these equations into a dimensionless normalization for the complex WFT coefficients using $\boldsymbol {B}' = \delta \hat {\boldsymbol {B}}/B_0$, $\omega ' = \omega /\varOmega _i$ and $\boldsymbol {E}' = \delta \hat {\boldsymbol {E}}/(v_A B_0)$, where $B_0$ is the magnitude of the upstream magnetic field and we use a dimensionless wavevector $\boldsymbol {k} d_i$, noting that $d_i = v_A/\varOmega _i$ should be the local value of the ion inertial length when comparing with linear wave modes. Thus, the components of Faraday's law may be expressed in this dimensionless normalization by

(3.8)\begin{gather} \omega' B'_{{\perp} 1} =-k_\parallel d_i E'_{{\perp} 2}, \end{gather}

(3.9)\begin{gather}\omega' B'_{{\perp} 2} =k_\parallel d_i E'_{{\perp} 1} - k_{{\perp} 1} d_i E'_\parallel, \end{gather}

(3.10)\begin{gather}\omega' B'_{{\parallel}} = k_{{\perp} 1} d_i E'_{{\perp} 2}. \end{gather}

The complex WFT coefficients for the dimensionless fluctuating electric field $\boldsymbol {E}'$ and magnetic field $\boldsymbol {B}'$ from the simulation for the dominant wave mode are presented in table 1. We can use (3.8), (3.9) and (3.10) to obtain three separate determinations of the real frequency $\omega$, also presented in table 1. Note that, since the original fields $\boldsymbol {E}$ and $\boldsymbol {B}$ before the WFT transform are real, the complex WFT coefficients must satisfy the reality conditions $\hat {\boldsymbol{E}}(\boldsymbol{k})=\hat {\boldsymbol{E}}^*(-\boldsymbol{k})$ and $\hat {\boldsymbol{B}}(\boldsymbol{k})=\hat {\boldsymbol{B}}^*(-\boldsymbol{k})$, so combining the coefficients for a given $\boldsymbol{k}$ and its conjugate $-\boldsymbol{k}$ yields strictly real fields, and therefore we obtain a real value for $\omega$ from each equation.

Table 1. Dimensionless complex WFT coefficients for the fluctuating electric field $\boldsymbol {E}'$ and magnetic field $\boldsymbol {B}'$ for the local plane-wave mode $(k_{\perp 1}d_{i,0},k_{\perp 2}d_{i,0}, k_{\|}d_{i,0}) = (1.97,0.00,0.34)$, along with determinations of the real frequency $\omega$ from each of (3.8), (3.9) and (3.10) normalized to the upstream cyclotron frequency, $\varOmega _{i,0}$.

3.4. Comparison with linear wave modes

To determine the linear wave mode associated with the dominant wave vector identified through our WFT analysis at $x/d_{i,0}=39.875$, we will compare the real frequency computed from Faraday's law with solutions for the frequency from the Vlasov–Maxwell linear dispersion relation for that wave vector using the PLUME dispersion relation solver (Klein & Howes Reference Klein and Howes2015). To do so, we must employ the plasma conditions locally at $x/d_{i,0}=39.875$ in the frame of the plasma. Relative to the upstream values of magnetic field, density and ion and electron temperatures, the local parameters have values $B/B_0=1.709$, $n_i/n_0=n_e/n_0=2.045$, $T_i/T_{i0}=3.510$ and $T_e/T_{e0}=1.611$. Computing the local dimensionless parameters yields $\beta _i=2.240$, $T_i/T_e=2.179$, $k_\parallel d^{({\rm loc})}_i=0.2378$ and $k_{\perp 1} d^{({\rm loc})}_i=1.371$, where the ratio of the local ion inertial length to its upstream value is $d^{({\rm loc})}_i/d_{i,0}=(n_0/n_i)^{1/2}=0.699$. Similarly, the frequency from (3.9) must be converted to a normalization relative to the local ion cyclotron frequency $\varOmega ^{({\rm loc})}_i/\varOmega _{i,0} =(B/B_0)$, yielding a wave frequency $\omega /\varOmega ^{({\rm loc})}_i=0.404$.

The solutions from the PLUME solver for scans over $k_{\|} d^{({\rm loc})}_i$ and $k_{\perp 1} d^{({\rm loc})}_i$ are presented in figure 7. We compare these PLUME solutions with empirical dispersion relations from Klein et al. (Reference Klein, Howes, TenBarge, Bale, Chen and Salem2012) and Howes, Klein & TenBarge (Reference Howes, Klein and TenBarge2014) for a kinetic Alfvén wave valid in the limit $k_\parallel d_i \ll 1$

(3.11)\begin{equation} \frac{\omega}{\varOmega_i}=k_{\|}d_i\sqrt{1+\frac{(k_{{\perp}} d_i)^2}{1+2/\beta_i(1+T_e/T_i)}}, \end{equation}

and from the MHD dispersion relation for the fast and slow magnetosonic waves valid in the limit $k d_i \ll 1$

(3.12)\begin{align} \frac{\omega}{\varOmega_i} = k d_i \left[\frac{1+\beta_i(1+T_e/T_i)\pm\sqrt{[1+\beta_i(1+T_e/T_i)]^2-4\beta_i (1+T_e/T_i)(k_{\|}^2/k^2)}}{2} \right],\end{align}

where the upper (lower) sign corresponds to the fast (slow) magnetosonic wave.

Figure 7. Comparisons between the measured frequency and wavelength of the ripple on the surface of the shock using empirical dispersion relations (Klein et al. Reference Klein, Howes, TenBarge, Bale, Chen and Salem2012; Howes et al. Reference Howes, Klein and TenBarge2014) (dotted lines), and dispersion relation for using PLUME, a Vlasov–Maxwell linear dispersion relation solver (Klein & Howes Reference Klein and Howes2015) (solid lines). Blue points are measured wavelength and frequency of the dominant wave mode. The plotted dispersion relations are kinetic Alfvèn wave (black), fast magnetosonic wave (red, lowest frequency), first three ion Bernstein modes (red, higher frequencies) and slow magnetosonic wave (green). Frequency is normalized to the local ion cyclotron frequency $\varOmega ^{({\rm loc})}_i$ and wavelength is normalized to the local ion inertial length $d^{({\rm loc})}_i$. (a) Dispersion relation as a function of $k_{\|} d^{({\rm loc})}_i$ at fixed $k_{\perp } d^{({\rm loc})}_i= 1.371$. (b) Dispersion relation as a function of $k_{\perp } d^{({\rm loc})}_i$ at fixed $k_{\|} d^{({\rm loc})}_i = 0.2378$.

In figure 7(a,b), we plot the PLUME solutions for the fast magnetosonic wave (lowest frequency solid red curve) along with the three lowest ion Bernstein modes (the $n=1$ through $n=3$ modes are the successively higher frequency solid red lines). In panel (b), the MHD fast wave dispersion relation (3.12) (dotted red) agrees well with the fast magnetosonic wave for up to $k_\perp d^{({\rm loc})}_i=0.5$ and $\omega /\varOmega ^{({\rm loc})}_i=1$, at which point the fast magnetosonic wave undergoes a mode conversion to the $n=1$ ion Bernstein mode. Note that, as $k_\perp d^{({\rm loc})}_i$ increases, the fluid solution for the fast magnetosonic wave (3.12) corresponds to the regions of increasingly higher $n$ ion Bernstein modes where those modes have a positive perpendicular group velocity $\partial \omega /\partial k_\perp$, as has been found in previous kinetic studies of the fast mode at nearly perpendicular propagation (Swanson Reference Swanson1989; Stix Reference Stix1992; Li & Habbal Reference Li and Habbal2001; Howes Reference Howes2009). We also plot the PLUME solutions for the kinetic Alfvèn wave (solid black), which show good agreement with the empirical kinetic Alfvèn wave solution (3.11) (dotted black), and we plot the MHD solution for the slow magnetosonic wave (dotted green).

To identify the wave mode responsible for the rippling of the shock front, we plot on figure 7 the wave frequency $\omega /\varOmega ^{({\rm loc})}_i=0.404$ calculated from (3.9) (blue point). It is clear from this comparison that the Alfvèn mode agrees with our frequency computation from the dominant unstable wave vector. We argue that this Alfvénic nature of the unstable fluctuation identifies the instability responsible for the shock ripple as the AIC for four reasons. First, the AIC instability launches waves on the Alfvén branch (Alfvén waves, kinetic Alfvén waves, ion cyclotron waves) (Tajima, Mima & Dawson Reference Tajima, Mima and Dawson1977), which we have just shown are present in this simulation. Second, figure 8 shows that we have a partial ring distribution, which matches the form of a perturbed distribution that was used to derive the AIC instability in Otani (Reference Otani1988). Third, figure 8 also shows temperature anisotropy $T_{\perp }/T_{\|}>1$ in the incoming beam (if a bi-Maxwellian distribution is fit to the multiple populations present) that is required for the AIC instability. This result is consistent with previous studies (Winske & Quest Reference Winske and Quest1988; McKean et al. Reference McKean, Omidi and Krauss-Varban1995; Burgess et al. Reference Burgess, Hellinger, Gingell and Trávníček2016). Fourth, our measured temperature anisotropy, $T_{\perp }/T_{\|}=2.01$, is above the marginal stability threshold established in Hellinger et al. (Reference Hellinger, Trávnıček, Kasper and Lazarus2006), $T_{\perp }/T_{\|}>1.67$, for the AIC.

Figure 8. Distribution of ions in the ramp ($x/d_{i,0} = 39.875$) of a $M_A \approx 7.88$ and $\theta _{B_n} = 45^\circ$ from a three-dimensional dHybridR simulation. The projection onto $v_x$,$v_z$ (b) shows three distinct populations of ions: the stream on the bottom, the population of having been reflected once in the middle and the doubly reflected ions on top.

One should note, computing $\omega$ with (3.8), (3.9) and (3.10) using measurements of $\boldsymbol {E}^\prime$ and $\boldsymbol {B}^\prime$ at some fixed $\boldsymbol {k}$ will not necessarily give similar values of $\omega$. This can occur as there may be multiple wave modes with different polarization may be superimposed. Figure 9 shows the polarization of a kinetic Alfvén wave and kinetic fast magnetosonic wave using the same parameters as those measured at $x/d_{i,0} = 39.875$ in the simulation. For $(k_{\perp 1} d^{({\rm loc})}_i,k_{\perp 2} d^{({\rm loc})}_i ,k_{\|} d^{({\rm loc})}_i)= (1.37,0.00,0.238)$, the amplitudes of $E_{\|}$ and $E_{\perp 2}$ are very small relative to $E_{\perp 1}$ for the kinetic Alfvén wave, suggesting that (3.9) will give the best estimate of the wave frequency, while (3.8) and (3.10) are worse choices for computing the frequency of an Alfvén wave with this wave vector. Thus, we choose the frequency computed using (3.9), $\omega = 0.69 \varOmega _i$, as the frequency of the dominant wave mode underlying the rippling of the shock surface.

Figure 9. Magnitude of the components of the eigenfunction response of a kinetic Alfvèn wave (a,b) and kinetic fast magnetosonic wave for the electric and magnetic fields (c,d) in a homogeneous plasma computed using PLUME, normalized to $|E_{\perp,1}|$, using local parameters at $x/d_{i,0} = 39.875$, as discussed in § 3.4. The vertical black line corresponds to the perpendicular wavenumber of the dominant wave mode discussed in § 3.4.

We close with a brief discussion of the validity of using the Vlasov–Maxwell linear dispersion relation, derived for spatially homogeneous plasma conditions with no bulk plasma flow, to interpret the results of kinetic numerical simulations of a collisionless shock in which the plasma is strongly spatially inhomogeneous with a rapid flow of the plasma through the approximately steady-state structure of the shock transition. The study of how kinetic instabilities arise and saturate in the presence of other dynamical plasma processes – such as plasma turbulence, magnetic reconnection and collisionless shocks – represents an important area of research at the frontier of plasma physics. For example, a recent analysis of spacecraft measurements from the Solar Orbiter mission (Müller et al. Reference Müller2020) argues that the complex interaction between turbulent fluctuations and kinetic instabilities ultimately regulates the proton-scale energetics of the solar wind (Opie et al. Reference Opie, Verscharen, Chen, Owen and Isenberg2023). In the case of the collisionless shock investigated here, we identify the wavevector of the dominant fluctuation in the shock-rest frame, and then use Faraday's law to calculate the frequency of the plasma response in that frame of reference. It is important the emphasize that, in the shock-rest frame, the non-Maxwellian ion velocity distributions shown in figure 8 are relatively steady in time, and therefore it seems plausible that a sufficiently rapidly growing instability (with a growth rate of order the ion cyclotron period) will be able to arise locally from the free energy in the unstable ion velocity distribution. Even in the presence of gradients of the plasma density, plasma temperature and magnetic field, we still expect the linear response of the plasma to give rise to the wave-like fluctuations appropriate for the plasma parameters at that position. An unanswered question, that is beyond the scope of this investigation, is how the super-Alfvénic flow of the plasma through shock transition competes with the growth of a kinetic instability, but the kinetic numerical simulation results presented here clearly show that, for the shock parameters $M_A \simeq 7.88$ and $\theta _{Bn}=45^\circ$, a kinetic instability is indeed able to give rise to shock rippling even in the presence of the bulk flow.

4.1. Steady-state shock energization

We present our field–particle correlation analysis of the rate of ion energization due to the steady-state shock dynamics using the averaged correlations $\overline {C}_{E_j}(\boldsymbol {v})$ for each component of the averaged electric field $\overline {E}_j$ in figure 10. The correlation is computed over an integration volume with normal thickness $\Delta x=0.25 d_{i,0}$, centred at $x_0/d_{i,0}=39.875$, and area covering the full transverse extent of the domain $L_y\times L_z=12 d_{i,0} \times 12 d_{i,0}$. We separate the energization by each of the three components of the averaged electric field in the simulation coordinates, $(\overline {E}_x,\overline {E}_y,\overline {E}_z)$, and present three views of the three-dimensional velocity space using integrations along each of the three velocity space dimensions, for a total of nine panels. For example, the reduction along the $v_z$ dimension for the $\overline {E}_y$ averaged correlation is given by $\overline {C}_{E_y}(v_x,v_y) = \int {\rm d}v_z \overline {C}_{E_y}(v_x,v_y,v_z)$. In addition, we plot the same three reduced views of the ion velocity distribution $f_i(\boldsymbol {v})$ in the top row, using a logarithmic colour map to highlight the different small populations of reflected particles.

Figure 10. (a–c) Total velocity distribution function of ions $f_i(x,\boldsymbol {v})$ at the transition from the foot to the ramp of the shock at $x_0/d_{i,0} = 39.875$. (d–l) Average energization from the field–particle correlation $\overline {C}_{E_j}(x_0,\boldsymbol {v})$, where we integrate the three-dimensional velocity space over the third velocity-space coordinate in each column. All quantities are computed in the shock-rest frame and normalized to $n_0$, the upstream particle density.

Note that the net energization by each of the three field components is most clearly visualized along the reverse diagonal of the nine panel correlation plots, panels (f,h,j). The dependence of the correlation (4.7) on $\partial \overline {f}/\partial v_j$, combined with the fact that this derivative must pass through zero along the $v_j$ axis, implies that the correlation with $\overline {E}_j$ consists of positive and negative regions along the $v_j$ axis. Using the red–white–blue colour map, it can be difficult to assess by eye the net rate of energization. But upon the integration over the $v_j$ dimension – given by panels (f,h,j) on the reverse diagonal – the net energization is easily observed. The quantitative value of the net rate of energization by the $\overline {E}_j$ over the integration volume, $\langle \,j_{j,s} \overline {E}_j \rangle$, is also plotted at the top of the reverse diagonal panels.

Furthermore, to assist in interpretation of the velocity-space signatures, we note that when a bipolar signature is observed (consisting of adjacent blue and red regions), the weighting by $v_j^2$ in the correlation means that the colour further from the origin (at higher velocity magnitude $|v_j|$) usually dominates. To standardize our terminology for bipolar signatures, we will state the colour with the lower magnitude of velocity $|v_j|$ first. Therefore, a blue–red bipolar signature generally denotes net particle energization, and a red–blue bipolar signature typically indicates a net loss of particle energy. We can interpret bipolar signatures as particles being accelerated/decelerated by the $E_j$ component of the electric field from the ‘blue’ region to the ‘red’ region in a collisionless system.

The contribution from averaged cross-shock electric field $\overline {E}_x$, evaluated at $x_0/d_{i,0}=39.875$, to the rate of change of phase-space energy density due to the steady-state shock energization is given by the velocity-space signatures in figure 10, panels (d–f). Panels (d,e) show the deceleration of the incoming ion beam in the $(v_x,v_y)$ and $(v_x,v_z)$ projections, with a red–blue bipolar signature indicating a net loss of energy by the incoming beam population. This net loss of ion energy by $\overline {E}_x$ is clearly shown in panel (f) as the dark blue region at $(v_z,v_y) \simeq (0,0)$. Also visible in panel (f) is the small population of reflected ions heading back upstream ($v_x>0$) at $v_y/v_{ti} \simeq -5$ that gain energy from the cross-shock electric field, and a subsequent loss of energy by $\overline {E}_x$ as those reflected particles return downstream at $v_y/v_{ti} \simeq -9$ and $3\lesssim v_z/v_{ti} \lesssim 7$.

The primary mechanism of ion acceleration at quasi-perpendicular collisionless shocks is mediated by the motional electric field, $\overline {E}_y$, with the three viewpoints of $\overline {C}_{E_y}$ shown in panels (g–i). In panel (g), the dominant velocity-space signature is the blue–red ‘crescent’ indicating acceleration of the reflected ion population by the motional electric field (Juno et al. Reference Juno, Howes, TenBarge, Wilson, Spitkovsky, Caprioli, Klein and Hakim2021, Reference Juno, Brown, Howes, Haggerty and TenBarge2022), a mechanism known as shock-drift acceleration (Paschmann et al. Reference Paschmann, Sckopke, Bame and Gosling1982; Sckopke et al. Reference Sckopke, Paschmann, Bame, Gosling and Russell1983; Anagnostopoulos & Kaliabetsos Reference Anagnostopoulos and Kaliabetsos1994; Anagnostopoulos et al. Reference Anagnostopoulos, Rigas, Sarris and Krimigis1998; Ball & Melrose Reference Ball and Melrose2001; Anagnostopoulos, Tenentes & Vassiliadis Reference Anagnostopoulos, Tenentes and Vassiliadis2009; Park et al. Reference Park, Ren, Workman and Blackman2013). The same blue–red crescent signature is visible from a different viewpoint in panel (i). The net ion energization due to shock-drift acceleration by the motional electric field is most easily viewed in panel (h), where the red diagonal feature is an additional velocity-space signature of shock-drift acceleration (Juno et al. Reference Juno, Brown, Howes, Haggerty and TenBarge2022) and is perpendicular local magnetic field direction at this point $x/d_{i,0}=39.875$, where the shock transitions from the foot to the ramp. Integration of $\overline {C}_{E_y}$ over velocity space shows that the net rate of positive ion energization, $\langle j_y \overline {E}_y \rangle$, is dominated by the acceleration of the reflected ion population by the averaged motional electric field. A key result of our field–particle correlation analysis of the steady-state shock energization is the velocity-space signature of shock-drift acceleration, shown in panels (g)–(i), together indicating the acceleration of the reflected ion population by upstream the motional electric field for a $\theta _{Bn}=45^\circ$ shock.

The steady-state shock energization due to the self-consistently generated $\overline {E}_z$ is shown in panels (j)–(l). In panel (k), we show the energization of three separate populations of ions: (i) a loss of energy for the incoming ion beam at $v_z \simeq 0$; (ii) a small population of ions that have been reflected once in the range $2 \lesssim v_z/v_{ti} \lesssim 12$; and (iii) a smaller population of ions that have been reflected a second time from the shock (and may be escaping upstream, as discussed in Juno et al. Reference Juno, Brown, Howes, Haggerty and TenBarge2022) at $12 \lesssim v_z/v_{ti} \lesssim 18$. The net loss of energy by the incoming beam and net gain of energy by the two populations of reflected ions by $\overline {C}_{E_z}$ is clearly seen in panel (j). Note, however, that the net rate of energization by $\overline {E}_z$ is more than a factor of 30 smaller than that by $\overline {E}_x$ and $\overline {E}_y$.

Using this analysis, we can analyse the energization of each population of particles separately. With the field–particle correlation technique, we can clearly separate the net loss of energy of the incoming ion beam from the energization of the reflected ion population. A fluid description of particle energization providing only $\boldsymbol {j}_i \boldsymbol{\cdot} \boldsymbol {E}$, a quantity integrated over velocity space, loses this clear separation of energization mechanisms operating on distinct populations of ions in velocity space.

4.2. Instability energization

We present our field–particle correlation analysis of the rate of ion energization due to the kinetic instabilities driven at the shock using the fluctuating correlations $\widetilde {C}_{E_j}$ for each component of the fluctuating electric field $\delta E_j$ in figure 11. These results are presented using the same format as used for the steady-state shock energization in figure 10, enabling direct comparisons of features and amplitudes of the energization rates.

Figure 11. (a–c) Total velocity distribution function of ions $f_i(x,\boldsymbol {v})$ at the transition from the foot to the ramp of the shock at $x_0/d_{i,0} = 39.875$. (d–l) Instability energization from the fluctuating correlation $\widetilde {C}_{E_j}(x_0,\boldsymbol {v})$, where we integrate the three-dimensional velocity space over the third velocity-space coordinate in each column. All quantities are computed in shock-rest frame.

First, we note that the net rates of energization due to the three fluctuating electric field components $\delta E_j$ is more than two orders of magnitude smaller than the rate of energization by the steady-state physics of the shock transition. Second, the rate of energy transfer from all three components of $\widetilde {C}_{E_j}$ is positive, meaning that the instability-generated fluctuations are being damped at this position, giving their energy to the particles. In the next section, we will explore how the net rate of energization by instabilities varies as a function of the normal distance through the shock.

The key new qualitative feature of interest in these fluctuations correlations $\widetilde {C}_{E_j}$ is the appearance of a ‘tripolar’ signature, such as that arising from the fluctuating normal component of the electric field $\delta E_x$ in panels (d,e). We do not observe any such tripolar features due to the steady-state shock physics in figure 10. Such a tripolar feature is a consequence on the non-uniformity of the electric field across the transverse domain, e.g. shock ripple, as we explain below. This tripolar feature indicates increasing phase-space energy at velocities flanking the average velocity of the incoming stream.

The tripolar feature in figure 10(d) can be explained by computing the correlation using sub-regions across the transverse plane of the simulation. In the top panel of figure 12, we plot the normal component of the fluctuating electric field $\delta E_x$ across the full transverse plane $(y,z)$ at $x/d_{i,0}=39.875$. We then divide this plane up into 16 subregions of size $3 d_{i,0} \times 3 d_{i,0}$ (white lines) and compute the fluctuating correlation $\widetilde {C}_{E_x}(v_x,v_y)$ in each of these sub-regions, plotted in the lower panel. The second column of the lower panel illustrates how the sum of each of these sub-regions combines to produce the tripolar signature. Numbering the four elements of the second column from top to bottom, the first and third elements have $\delta E_x>0$ and yield a blue–red bipolar (positive energization) signature, while the second element has $\delta E_x<0$ and yields a red–blue bipolar (negative energization) signature. The fourth element has a mixture of $\delta E_x>0$ and $\delta E_x<0$ within the sub-region, and yields a red–blue–red tripolar signature. Critically, the $v_x$ value of the zero crossing of the bipolar signatures in the first three elements shifts slightly to right in first and third elements and slightly to the left in the second element. (The vertical dotted green line, at the average velocity of the incoming ion stream in the shock-rest frame assists in visualization of this small shift.) Therefore, their sum would yield a net positive (red) value at the lowest $|v_x|$ values, a net negative (blue) at intermediate $|v_x|$ values and a positive (red) value at the highest $|v_x|$ values – their sum would then be a red–blue–red tripolar signature, with a net positive rate of ion energization. Thus, the qualitative appearance of a tripolar signature indicates energization by the fluctuating electric field (across the transverse plane) associated with an instability, and the velocity-space-integrated rate of energization yields the net ion energization rate by that instability.

Figure 12. (a) Fluctuating cross-shock electric field $\delta E_x$ at $x_0/d_{i,0} = 39.875$, plotted across the transverse plane, with subregion boundaries indicated (white lines). (b) Fluctuating correlation $\widetilde {C}_{E_x}(v_x,v_y)$, generated using spatial subregions of size $3 d_{i,0} \times 3 d_{i,0}$ in the $(y,z)$ plane, with the same range $\Delta x=0.25 d_{i,0}$ centred at $x_0/d_{i,0} = 39.875$ as used in figures 10 and 11. The dotted green vertical line on each panel indicates the average velocity of the incoming ion stream in the shock-rest frame.

Using our separation of the steady-state shock physics from the instability physics, we find that the instability driving the rippling of the shock surface in the ramp accounts for $\lesssim 1$% of the total ion energization. Despite this small amplitude of energy transfer, we are able to isolate the dominant wave mode causing the shock ripple, distinguish shock energization from instability energization using the velocity-space signatures and quantitatively analyse the rate of energy transfer to the ions in velocity space using the field–particle correlation technique. With these methods, we can analyse the potentially different energization of distinct populations of particles, even for weak mechanisms.

4.3. Steady-state shock and instability energization through the shock

The profile of energization of ions throughout the shock due to the averaged fields and the fluctuating fields is shown in figure 13, where the total rate of energization can be computed from integrating the field–particle correlation over velocity space using

(4.9)\begin{equation} \int d^3\boldsymbol{v} C_{\boldsymbol{E}} (\boldsymbol{v}) = \boldsymbol{j} \boldsymbol{\cdot} \boldsymbol{E} . \end{equation}

Thus, we can compute energization due to the steady-state shock physics and to the instability physics separately by computing $\int d^3\boldsymbol {v} \widetilde {C}_{\boldsymbol {E}} (\boldsymbol {v})$ and $\int d^3\boldsymbol {v} \overline {C}_{\boldsymbol {E}} (\boldsymbol {v})$, respectively.

Figure 13. Comparison of energy transfer between the particles and transverse-plane averaged fields(a) as well as between that particles and the fluctuating fields (b) and kinetic energy/thermal energy of ions (c) near the shock transition region in the shock-rest frame. The vertical black line indicates the position that the instability isolation method was applied in § 3. Quantities are normalized as specified in § 2 and by $n_0$, the upstream particle density.

The $\int d^3\boldsymbol {v} \overline {C}_{E_x}(x) (\boldsymbol {v})$ term (top panel, red dotted) shows deceleration of the ions (loss of energy) by the cross-shock electric field throughout the shock transition region, $44 \gtrsim x/d_{i,0} \gtrsim 38$. We see net acceleration of ions (gain of energy) throughout the same shock transition region by $\overline {E}_y$ in the $\int d^3\boldsymbol {v} \overline {C}_{E_y}(x) (\boldsymbol {v})$ term (top panel, green dashed), which accelerates the reflected ions through the shock-drift acceleration mechanism, as demonstrated by our field–particle correlation analysis in figure 10(h). The work done by the $\overline {E}_z$ component of the electric field on the ions is generally negligible relative to that due to the other components.

The relative rate of energization of the ions due to the instability-driven fluctuating fields is small compared with the steady-state shock energization throughout the entirety of the shock transition region, with amplitudes often much less than 10 % of the steady-state shock energization. The ion kinetic instabilities that generate the electromagnetic field fluctuations observed in this simulation arise by tapping free energy in unstable ion velocity distributions caused by the steady-state shock physics in the shock transition. Therefore, the condition that the velocity-integrated total fluctuating correlation is negative, $\int d^3\boldsymbol {v} \widetilde {C}_{\boldsymbol {E}}(\boldsymbol {v})<0$, indicates clearly the regions where these instabilities arise. As shown in figure 13(b), the ions drive these instabilities over the range $37 \lesssim x/d_{i,0} \lesssim 39.5$, which corresponds to the ramp and overshoot regions of the shock. These instability-driven fluctuations appear to propagate a short distance upstream into the foot of the shock, before they are fully damped out upon reaching $x/d_{i,0} \gtrsim 42$. Thus, at the position $x_0/d_{i,0}=39.875$ analysed in our field–particle correlation analysis, the unstable electromagnetic fluctuations are being damped, leading to a small, but positive, net energization of the ions at that position, as shown in figure 11. Unstable electromagnetic fluctuations may also be swept downstream to $x/d_{i,0} \lesssim 37$ with the bulk plasma flow, appearing as turbulent fluctuations downstream that we expect will eventually be damped and ultimately lead to heating of the downstream plasma.