Isolation and phase-space energization analysis of the instabilities in collisionless shocks

Abstract

We analyse the generation of kinetic instabilities and their effect on the energization of ions in non-relativistic, oblique collisionless shocks using a 3D-3V (three spatial with three velocity components) simulation by dHybridR, a hybrid particle-in-cell code. At sufficiently high Mach number, quasi-perpendicular and oblique shocks can experience rippling of the shock surface caused by kinetic instabilities arising from free energy in the ion velocity distribution due to the combination of the incoming ion beam and the population of ions reflected at the shock front. To understand the role of the ripple on particle energization, we devise a new instability isolation method to identify the unstable modes underlying the ripple and interpret the results in terms of the governing kinetic instability. We generate velocity-space signatures using the field–particle correlation technique to look at energy transfer in phase space from the isolated instability driving the shock ripple, providing a viewpoint on the different dynamics of distinct populations of ions in phase space. Together, the field–particle correlation technique and our new instability isolation method provide a unique viewpoint on the different dynamics of distinct populations of ions in phase space and allow us to completely characterize the energetics of the collisionless shock under investigation.

1. Introduction

Collisionless shocks above a critical Mach number must invoke energy dissipation mechanisms other than resistivity to enable a steady-state structure in the transition region (Leroy et al. 1981; Wu et al. 1984). Such supercritical shocks (Edmiston & Kennel 1984; Kennel, Edmiston & Hada 1985) can form when a large obstacle is put into a plasma flow that is travelling significantly faster than the fast magnetosonic velocity in the frame of the obstructing object. This disturbance creates a plasma wave that steepens until a shock forms, creating a transition from supersonic flow upstream to subsonic flow downstream. Collisionless shocks impact the dynamics and energetics of a plasma system in several ways, such as the ability to accelerate particles to high energy (Caprioli, Zhang & Spitkovsky 2018), and the partitioning of energy between electrons and ions (Savoini & Lembege 2010), which is allowed due to the lack of particle collisions to drive ion and electron temperatures to equilibrium.

For quasi-perpendicular and oblique shocks, a fraction of incoming ions are reflected at the shock transition and travel back upstream a distance that is typically of order one gyroradius, creating an unstable distribution of incoming and reflected ions that can generate unstable electromagnetic fluctuations which mediate the dissipation of energy and aid in forming the structure of the shock transition. This reflection plays a significant role in creating the foot/ramp/overshoot structure seen in the transverse (to the shock normal direction) magnetic field that is typical of collisionless shocks (Leroy et al. 1981, 1982; Balogh & Treumann 2013; Burgess & Scholer 2015). A gradual increase in the transverse magnetic field occurs in the foot, with a more rapid increase in the ramp of the shock. Beyond this, the overshoot increases the transverse component of the magnetic field above its eventual downstream asymptotic value.

For shocks with a Mach number and shock normal angle in a particular regime, the surface of the shock starts to ripple. This occurs as kinetic instabilities are driven by the unstable distribution created by the combination of the incoming and reflected populations of ions. This has been observed in simulation (Winske & Quest 1988; McKean, Omidi & Krauss-Varban 1995; Lowe & Burgess 2003; Burgess et al. 2016) and in situ with spacecraft (Johlander et al. 2016, 2018). Shock rippling has been studied previously using two-dimensional (Winske & Quest 1988; McKean et al. 1995; Lowe & Burgess 2003) and three-dimensional simulations (Burgess et al. 2016). Studies that simulate only two spatial dimensions potentially risk suppressing some of the important degrees of freedom of the system (Burgess & Scholer 2007; Zacharegkas et al. 2022). Here, we focus on the ion-scale fluctuations caused by Alfvénic modes. The shock ripple may also impact other unstable fluctuations arising in the shock transition, for example the rippling may interact with whistler waves that propagate upstream and downstream (Gedalin & Ganushkina 2022).

Modelling this instability with the theory of linear waves in a homogeneous plasma has been difficult, as the instability that best matches the observations depends on the position relative to the shock transition. In the shock foot, the instability appears to be driven by streaming instabilities as the reflected population has a significant component of its velocity perpendicular to the velocity of the incoming stream, and this has been proposed to be due to either the modified two stream instability (Winske & Quest 1988) or the ion Weibel instability (Burgess et al. 2016). In the shock ramp, previous studies have concluded that the instability corresponds to the Alfvén ion cyclotron instability (AIC), as the waves generated locally tend to travel at the local Alfvén velocity and the plasma meets the temperature anisotropy threshold, factors that are also required by the AIC instability (Gary et al. 1976; Winske & Quest 1988; Lowe & Burgess 2003; Klein & Howes 2015). Agreement between linear theory with the AIC instability has been found (McKean et al. 1995), but the analysis methods used make it difficult to isolate instabilities in a restricted time or space domain. This difficulty is further enhanced in higher energy shocks in which the instabilities excited can cause non-stationary, non-steady-state behaviour such as shock reformation because of how rapidly these instabilities evolve in space and time. Furthermore, modelling this instability is also complicated by the non-homogeneous nature of the plasma through the shock transition.

Previous work has compared the linear theory of kinetic instabilities in homogeneous plasmas with two-dimensional shock simulations (Winske & Quest 1988; McKean et al. 1995; Burgess & Scholer 2007). These studies tracked the shock as a function of time to measure the growth of any expected instabilities, or Fourier transformed a block of the simulation in space and time, which loses the locality of the dynamics of the shock along the shock normal direction and the locality of the plasma parameters. While these attempts have been successful at examining the dynamics in the ramp, to our knowledge, no direct comparisons between linear theory and the instabilities in the foot have been made. Only quantitative similarities between properties of both the instabilities in the foot and linear instabilities that are suspected to be present have been found (Burgess et al. 2016). This motivates the need for development of a more versatile method for comparing simulations with linear theory, as such a tool will have applications beyond analysis of the instabilities causing shock ripple.

The purpose of this paper is to develop tools for analysing how kinetic instabilities affect the dynamics and energetics of a collisionless shock and to demonstrate these tools on the main instability causing the surface of the shock to ripple. We develop a method that can locally identify the kinetic instability (or instabilities) present and characterize its properties from simulation data by analysing the ‘fluctuating fields’, as defined in § 3. We show that these fluctuations are consistent with the properties of wave modes from the linear Vlasov–Maxwell dispersion relation based on the local plasma parameters at that position, as shown in figure 7. Furthermore, we use this definition of the fluctuating fields to compute the energization due to the fluctuations arising from kinetic instabilities in the shock transition. Thus, we devise a means to separate the approximately steady-state bulk shock energization of particles due to the shock transition from the energization due to the instability-driven fluctuations, with the aim to quantify and explain how the electromagnetic fluctuations arising from kinetic instabilities affect the energization of particles through the shock transition.

In this paper, we will present a novel instability isolation method and apply it to a simulation of an oblique collisionless shock. We use this method in conjunction with the field–particle correlation technique to produce the velocity-space signatures of particle energization in phase space due to an isolated instability. We use these methods to investigate the kinetic instability responsible for the shock ripple and assess the impact of the instability on the ion energization in the shock. These methods allow us to separate the energization of particles due to the steady-state physics of the shock from that due to the kinetic instabilities arising in the shock transition. In § 2, we describe the set-up and details of the 3D-3V dHybridR simulation of an oblique collisionless shock. Section 3 presents the instability isolation method devised to analyse the properties of the modes driven unstable by kinetic instabilities in the shock transition. We analyse the particle energization in the ramp of this simulated shock, separating the energization due to the bulk steady-state shock physics (a transverse-plane average of the fields through the shock) from the energization due to the kinetic instabilities responsible for the shock ripple in § 4. We conclude in § 5 by discussing new avenues of investigation made possible by these techniques to understand the dynamics and energetics of other non-stationary shocks.

2. Simulation set-up

We present fully three-dimensional hybrid simulations using the massively parallel code dHybridR (Gargaté et al. 2007; Haggerty & Caprioli 2019). We produce a shock by sending a supersonic flow in the $-\hat {x}$ direction towards a reflecting wall at $x=0$. The resulting colliding flows generate a shock that propagates in the $+\hat {x}$ direction in the simulation domain. We maintain periodic boundary conditions in the $\hat {y}$ and $\hat {z}$ directions. We use $1000$ particles per cell in the initial state of each cell and at the upstream boundary where particles are continuously injected. We normalize the magnetic fields and number density to their values in the upstream region, $B_0$ and $n_0$. Lengths are scaled to the upstream ion inertial length $d_{i,0} \equiv c / \omega _{pi,0}$, where $\omega _{pi,0} = \sqrt {4 {\rm \pi}n_0 e^2/m_i}$ is the ion plasma frequency using the upstream density with c as the speed of light, e as the elementary charge, and $m_i$ as the ion mass. Time is scaled to the inverse ion gyrofrequency $\varOmega _{i,0}^{-1} \equiv c m_i / e B_0$ based on the upstream magnetic field. Velocity is normalized using to the upstream Alfvén velocity, $v_{A,0} = B_0/\sqrt {4 {\rm \pi}n_0 m_i} = d_{i,0} \varOmega _{i,0}$. Electric fields are normalized to $B_0 v_{A,0} / c$. The ratio of the Alfvén velocity to the speed of light is $v_{A,0}/c = 1/125$. The plane normal to the shock propagation direction, the transverse plane, has an area of $L_y \times L_z = 12 d_{i,0} \times 12 d_{i,0}$. The length of the simulation along the direction of the shock normal is $L_x = 98.75 d_{i,0}$. The simulation is divided into square cells with sides of length $0.25 d_{i,0}$. We set the initial state of the inflowing plasma to have beta values $\beta _i = \beta _e = 1$, where the species $s$ plasma beta is given by $\beta _s = 8 {\rm \pi}n_s T_{s,0}/B_0^2$ and temperatures are expressed in units of energy (absorbing the Boltzmann constant into $T_{s,0}$, the upstream ion temperature). The upstream ion thermal velocity is defined by $v_{ti,0}= \sqrt {2 T_{i,0}/m_i}$. We inject particles at the upstream boundary with velocity $U_{inj} = -6 v_{A,0} \hat {x}$ in the frame of the simulation box and impose an initial magnetic field with shock normal angle $\theta _{B_n} = 45^\circ$ with $\boldsymbol {B}_0 = (1/\sqrt {2})B_0 \hat {x} + (1/\sqrt {2})B_0 \hat {z}$.

Figure 1 shows the magnetic field structure, ion density and bulk flow velocity at $\varOmega _{i,0} t = 20$. The incoming flow in the shock-rest frame (as inferred from the simulation) has an Alfvén Mach number $M_A=7.88$, and the shock does not experience any shock reformation or significant ‘breathing’, i.e. the shock velocity is stable throughout the simulation. We track the shock in the simulation frame as discussed in Appendix B, and Lorentz transform the electromagnetic fields and ion velocity distributions to the shock-rest frame. The electric and magnetic fields averaged across the transverse plane (in the shock-rest frame) are shown in figure 2. At $\varOmega _{i,0} t = 20$, the shock has reached a position $x/d_{i,0} \simeq 40$. The compressible component of the magnetic field, $B_z$, has a compression ratio of $r \approx 3.5$. This ratio is in approximate agreement with the predicted compression ratio of $r = 3.62$ using MHD (magnetohydrodynamics) Rankine–Hugoniot jump conditions (see Appendix B). In this parameter regime and at this time in the simulation, MHD and kinetic shocks are in good agreement, but simulations of kinetic shocks can deviate (Bret 2020; Caprioli, Haggerty & Blasi 2020; Haggerty & Caprioli 2020). We observe ions being reflected in the phase-space plot $f_i(x,v_x)$, shown in figure 3, creating the instability discussed in § 3. The distribution of ions evolves rapidly along the $\hat {x}$ direction near the shock transition region, creating a small range in $x$ where the distribution is unstable, as discussed in § 4.3.

Figure 1. Two-dimensional slice of magnetic fields and fluid moments of a dHybridR simulation with shock velocity of $M_A=7.88$ (in the downstream rest frame) and $\theta _{B_n} = 45^\circ$ shock at $\varOmega _{i,0} t = 20$. At this time, the shock is at $x / d_{i,0} \simeq 40$.

Figure 2. Average (over the transverse plane) magnetic fields (a) and electric fields (b) from the hybrid simulation at $\varOmega _{i,0} t = 20$. The black curve shows the transverse magnetic field jump $B_{z,2}/B_{z,1}=3.66$ predicted by the MHD Rankine–Hugoniot jump conditions for a collisionless shock with $M_A=7.88$, $\theta _{B_n} = 45^\circ$, and $\beta =2$.

Figure 3. Ion distribution function $f_i(x,v_x)$, integrated over $v_y$ and $v_z$, in the shock-rest frame. Ions are reflected back upstream a distance of order one $d_{i,0}$.

There is rippling in the shock at $x / d_{i,0} = 39.875$. This location corresponds to the most upstream extent of the ramp. Figure 4 shows a two-dimensional slice of the rippling over the transverse plane. We see a periodic rippling in the magnetic field, indicating the presence of an instability launching a wave that is rippling the location of the shock front across the transverse direction. The width of the ripple is estimated to be $\Delta x /d_{i,0} \simeq 1$ from figure 4. There is no variation in the transverse plane in the upstream region, showing that the unstable ion velocity distribution arising in the shock transition induces instabilities which produce a variation of the fields in transverse directions.

Figure 4. Two-dimensional slice of the total magnetic field, $|\boldsymbol {B}(x_0,y,z)|$, at $x_0/ d_{i,0} = 39.875$ (a) and the total electric field, $|\boldsymbol {E}(x_0,y,z)|$ (b) over the transverse plane. There is rippling of the shock with $k_y d_{i,0} = -0.52$ and $k_z d_{i,0} = -1.05$. Line slices of the compressible magnetic field component, $B_z(x,y_i,z_0)$ (c) and $E_x(x,y_i,z_0)$ (d) as a function of $x$ with a fixed value of $z_0/d_{i,0} = 0.125$ for a set of discrete values of $y_i/d_{i,0}$.

3. Instability isolation method

To isolate the unstable fluctuations from the variation of the electromagnetic fields due to the steady-state physics of the shock transition in our collisionless shock simulations, we first separate the average of these fields over the transverse plane of the shock (i.e. along the shock face) from the fluctuations across that plane, where the boundary conditions in the plane are periodic. We define the transverse-plane-averaged fields (hereafter denoted the averaged fields for brevity) by

(3.1)\begin{equation} \overline{\boldsymbol{E}}(x) \equiv \frac{1}{L_y L_z}\int_0^{L_y}\,{{\rm d} y} \int_0^{L_z}\,{\rm d}z \boldsymbol{E}(x,y,z).\end{equation}

The fluctuating fields are then computed by

(3.2)\begin{equation} \delta \boldsymbol{E}(x,y,z) = \boldsymbol{E}(x,y,z)-\overline{\boldsymbol{E}}(x),\end{equation}

where $\boldsymbol {E}(x,y,z)$ is the total electric field from the simulation in the shock-rest frame. The motivation for this separation of the average and fluctuating fields is that the variation of the fields and flows associated with compression through the shock transition, in the absence of instabilities or other time variation in the shock-rest frame, is inherently one-dimensional along the shock normal¹. We denote this one-dimensional physics of the averaged fields the steady-state shock physics, and the dynamics and energetics of the unstable modes are the instability physics. Note that the one-dimensional steady-state shock physics is only approximately steady state – for supercritical shocks, the structure through the shock transition does undergo slight oscillations in time in the shock-rest frame (sometimes denoted ‘breathing’ of the shock), oscillations which increase in amplitude as the Mach number increases. Upon reaching or exceeding the second critical Mach number, $M_{2c,A} \approx 15.1$ (Krasnoselskikh et al. 2002; Oka et al. 2006), these oscillations lead to the phenomenon of shock reformation, but as this simulation has a significantly lower Alfvén Mach number than this second critical threshold, $M_A< M_{2c,A}$, we expect no reformation in this simulation.

It should be noted that despite our separation of the steady-state shock physics from the instability physics using (3.1) and (3.2), these phenomena are causally related. That is, the steady-state shock physics, which has no inherent variation in the transverse plane, generates the unstable ion velocity distributions that give rise to the instabilities, so the energetics of the steady-state shock physics is not decoupled from the instability physics, but rather drives those instabilities.

The basic steps of the instability isolation method, using the electromagnetic fields in the shock-rest frame, are:

1. (i) Compute the averaged electric field $\overline {\boldsymbol {E}}(x)$ using (3.1) and the fluctuating electric field $\delta \boldsymbol {E}(x,y,z)$ using (3.2).

2. (ii) Fourier transform the fluctuating fields over the transverse plane $(y,z)$ to obtain $\widetilde {\delta \boldsymbol {E}} (x,k_y,k_z)$.

3. (iii) Perform a wavelet transform along the shock normal direction to obtain the wavelet-Fourier transform (WFT) $\delta \hat {\boldsymbol {E}} (k_x,k_y,k_z;x)$.

4. (iv) Use the same procedure as in (i)–(iii) above on the magnetic field to obtain the WFT of the magnetic field, $\delta \hat {\boldsymbol {B}} (k_x,k_y,k_z;x)$.

5. (v) At the position $x=x_0$ along the normal direction, compute the local averaged total magnetic field, $\overline {\boldsymbol {B}}(x_0)$.

6. (vi) Generate a local magnetic field-aligned coordinate (FAC) system $(\hat {\boldsymbol {e}}_{\perp 1},\hat {\boldsymbol {e}}_{\perp 1},\hat {\boldsymbol {e}}_{\parallel })$ using the normal direction $\hat {\boldsymbol {n}}=\hat {x}$ and $\overline {\boldsymbol {B}}(x_0)$.

7. (vii) Rotate the fluctuating WFT fields $\delta \hat {\boldsymbol {E}} (k_x,k_y,k_z;x_0)$ and $\delta \hat {\boldsymbol {B}} (k_x,k_y,k_z;x_0)$ into the FAC system.

8. (viii) To estimate the frequency of an unstable (local) plane-wave mode given by $(k_{\perp 1},k_{\perp 2},k_\parallel ;x_0)$ from the WFT transform, use Faraday's law in the FAC coordinate system.

9. (ix) Finally, compare the estimated frequency of the unstable mode from Faraday's law to the frequencies of the different wave modes from a linear dispersion relation solver.

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 1998) 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 1997) 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 (1998). 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 2015). 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. (2012) and Howes, Klein & TenBarge (2014) 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. 2012; Howes et al. 2014) (dotted lines), and dispersion relation for using PLUME, a Vlasov–Maxwell linear dispersion relation solver (Klein & Howes 2015) (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 1989; Stix 1992; Li & Habbal 2001; Howes 2009). 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 1977), 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 (1988). 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 1988; McKean et al. 1995; Burgess et al. 2016). Fourth, our measured temperature anisotropy, $T_{\perp }/T_{\|}=2.01$, is above the marginal stability threshold established in Hellinger et al. (2006), $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. 2020) argues that the complex interaction between turbulent fluctuations and kinetic instabilities ultimately regulates the proton-scale energetics of the solar wind (Opie et al. 2023). 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. Particle energization

The field–particle correlation technique allows us to understand the energetics in the 3D-3V phase of a kinetic plasma by analysing the transfer of energy between fields and particles (Klein & Howes 2016; Howes 2017; Howes, Klein & Li 2017; Klein, Howes & TenBarge 2017). It has been used successfully to explore energization of particles in kinetic turbulence simulations (Klein et al. 2017; Howes, McCubbin & Klein 2018; Li et al. 2019; Horvath, Howes & McCubbin 2020; Klein et al. 2020) and in situ with spacecraft (Chen, Klein & Howes 2019; Afshari et al. 2021), as well as electron energization in simulations of collisionless magnetic reconnection (McCubbin, Howes & TenBarge 2022), simulations of 1D-2V (one spatial with two velocity components) perpendicular collisionless shocks (Juno et al. 2021), and electron acceleration by Alfvén waves in laboratory plasmas (Schroeder et al. 2021).

The rate of change of phase-space energy density of species $s$, $w_s(\boldsymbol {r},\boldsymbol {v},t) \equiv m_s v^2 f_s(\boldsymbol {r},\boldsymbol {v},t)/2$, in a collisionless plasma can be determined by multiplying the Vlasov equation by $m_s v^2/2$, yielding

(4.1)\begin{equation} \frac{\partial w_s}{\partial t} =- \boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} w_s - q_s\frac{v^2}{2} \boldsymbol{E} \boldsymbol{\cdot} \frac{\partial f_s}{\partial \boldsymbol{v}} - \frac{q_s}{c}\frac{v^2}{2} (\boldsymbol{v} \times \boldsymbol{B}) \boldsymbol{\cdot} \frac{\partial f_s}{\partial \boldsymbol{v}}. \end{equation}

On the right-hand side of (4.1), only the middle term with the electric field leads to a net change of energy of the particles (Klein & Howes 2016; Howes et al. 2017; Klein et al. 2017), so we define the field–particle correlation by

(4.2)\begin{equation} C_{\boldsymbol{E}} (\boldsymbol{v},t) \equiv \left\langle -q_s \frac{v^2}{2} \frac{\partial f_s (\boldsymbol{r},\boldsymbol{v},t)}{\partial \boldsymbol{v}} \boldsymbol{\cdot} \boldsymbol{E}(\boldsymbol{r},t) \right\rangle,\end{equation}

where the angled brackets denote that the correlation is computed using either an average over a correlation interval in time or an integration over a volume in space. Note that integration of the correlation over velocity space of (4.2) yields the rate of work done on species $s$ by the electric field, $\int d^3\boldsymbol {v} C_{\boldsymbol {E}} (\boldsymbol {v}) = \langle \,\boldsymbol {j}_s \boldsymbol{\cdot} \boldsymbol {E} \rangle$. Since we are interested in determining the impact of each particular field component on the energization of particles, we may split the equation by each directional component $j$

(4.3)\begin{equation} C_{E_j} (\boldsymbol{v}) \equiv \left\langle -q \frac{v_{j}^2}{2} \frac{\partial f (\boldsymbol{r},\boldsymbol{v},t)}{\partial v_j} E_j(\boldsymbol{r},t) \right\rangle , \end{equation}

where the other contributions to $v^2$ in (4.2) contribute zero net energization upon integration over velocity space (Juno et al. 2021).

The field–particle correlation technique generates velocity-space signatures that contain both quantitative and qualitative features used here to further understand particle energization in collisionless shocks. The technique can be employed to probe the interactions between the electric field and the different populations of particles, e.g. the incoming beam of particles or the reflected particles at a supercritical shock. Previous approaches to understanding particle energization in weakly collisional plasmas involve either tracking individual particles in phase space (Caprioli & Spitkovsky 2014) or reducing the fundamental ‘3D-3V’ behaviour of a kinetic plasma to the fluid quantity, such as $\boldsymbol {j} \boldsymbol{\cdot} \boldsymbol {E}$ (Gershman et al. 2017). These previous approaches have notable limitations in developing a full understanding of the energy transport in phase space: by following either only individual particles or velocity-space-integrated quantities, it is difficult to assess the interactions between the fields and distinct populations of particles in velocity space. The field–particle correlation technique is especially useful in supercritical shocks, where the effect of the electric field on the small population reflected particles may comprise a significant fraction of the particle energization at the shock.

To implement the field–particle correlation technique for use with simulation data from particle-in-cell (PIC) codes, a straightforward approach would be to specify a small spatial volume and bin all of the particles within that volume into velocity bins, generating a velocity distribution function within that volume $f(\boldsymbol {v})$. We would then be able to take the velocity derivatives of that distribution function and correlate them with the components of the electric field, as given by (4.3). However, in this approach, the electric field must be averaged over the spatial integration volume, losing the spatial dependence of the electric field within the volume. A solution to this issue is to use an alternative formulation of the correlation at time $t_0$ (Chen et al. 2019)

(4.4)\begin{equation} C^\prime_{E_j} (\boldsymbol{v}) = \langle q v_j f(\boldsymbol{r},\boldsymbol{v},t_0) E_j(\boldsymbol{r},t_0) \rangle, \end{equation}

where the electric field determined at the location of each particle and the brackets indicate an integration over a spatial volume. For shock problems, we typically choose a volume over a small range $\Delta x$ in the direction normal to the shock and covering the full transverse extent of the domain $L_y \times L_z$. As presented in Chen et al. (2019), this alternative correlation $C^\prime _{E_j}$ can then be used to compute the standard correlation $C_{E_j}$ through the relation

(4.5)\begin{equation} C_{E_j} (\boldsymbol{v}) =-\frac{v_j}{2} \frac{\partial C^\prime_{E_j}(\boldsymbol{v})}{\partial v_j} + \frac{C^\prime_{E_j}(\boldsymbol{v})}{2}.\end{equation}

It is important that the correlation is computed this way to retain the variation of the electric field within the integration volume. We elaborate on the full details of the procedure of computing the field–particle correlation using data from PIC codes in Appendix A.

To separate the contribution to the particle energization by the steady-state shock physics from the energization due to the instability physics, we separate the transverse-plane averaged contribution from the fluctuating contribution, as given by (3.1) and (3.2), for both the electric field $\boldsymbol {E}(x,y,z)=\overline {\boldsymbol {E}}(x)+ \delta \boldsymbol {E}(x,y,z)$ and the ion velocity distribution function $f(\boldsymbol {r},\boldsymbol {v})=\overline {f}(x,\boldsymbol {v})+ \delta f(\boldsymbol {r},\boldsymbol {v})$. Substituting these decompositions into (4.3), we obtain

(4.6)\begin{align} C_{E_j} (\boldsymbol{v})& = \left[ \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \overline{f} (x,\boldsymbol{v})}{\partial v_j} \overline{E}_j(x) \right\rangle + \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \overline{f} (x,\boldsymbol{v})}{\partial v_j} \delta E_j(\boldsymbol{r}) \right\rangle \right. \nonumber\\ & \quad + \left. \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \delta f (\boldsymbol{r},\boldsymbol{v})}{\partial v_j} \overline{E}_j(x) \right\rangle + \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \delta f (\boldsymbol{r},\boldsymbol{v})}{\partial v_j} \delta E_j(\boldsymbol{r}) \right\rangle \right]. \end{align}

If the spatial average indicated by the angle brackets (also indicated by line over variable) denotes integration over the full transverse extent of the domain in $(y,z)$, with periodic boundary conditions in that transverse plane, then the two middle terms on the right-hand side of (4.6) contribute nothing since $\langle \delta E_j(\boldsymbol {r}) \rangle =0$ and $\langle \partial \delta f (\boldsymbol {r},\boldsymbol {v})/\partial v_j \rangle =0$ by the definitions (3.1) and (3.2). Therefore, under the transverse-plane average, we obtain a clean separation of the steady-state shock physics from the instability physics, defining the steady-state shock energization through the averaged correlation, $\overline {C}_{E_j}$, given by

(4.7)\begin{equation} \overline{C}_{E_j}(\boldsymbol{v}) = \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \overline{f} (x,\boldsymbol{v})}{\partial v_j} \overline{\boldsymbol{E}}(x) \right\rangle, \end{equation}

and the instability energization through the fluctuating correlation, $\widetilde {C}_{E_j}$, given by

(4.8)\begin{equation} \widetilde{C}_{E_j}(\boldsymbol{v}) = \left\langle -q \frac{v_{j}^2}{2} \frac{\partial \delta f (\boldsymbol{r},\boldsymbol{v})}{\partial v_j} \delta E_j(\boldsymbol{r}) \right\rangle.\end{equation}

Furthermore, when we integrate over the full transverse extent of the domain, we may simply use the full ion velocity distribution function $f(\boldsymbol {r},\boldsymbol {v})$ in the derivative for both (4.7) and (4.8), since the difference vanishes upon the spatial integration. Thus, the averaged correlation $\overline {C}_{E_j}(\,f)$ captures the transfer of energy between the particles and the averaged fields and the fluctuating correlation $\widetilde {C}_{E_j}(\,f)$ captures the transfer of energy between the particles and the fluctuating fields arising from instabilities.

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. 2021, 2022), a mechanism known as shock-drift acceleration (Paschmann et al. 1982; Sckopke et al. 1983; Anagnostopoulos & Kaliabetsos 1994; Anagnostopoulos et al. 1998; Ball & Melrose 2001; Anagnostopoulos, Tenentes & Vassiliadis 2009; Park et al. 2013). 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. 2022) 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. 2022) 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.

5. Conclusion

We present an in-depth analysis of the shock ripple produced by a kinetic instability in a three-dimensional, hybrid PIC simulation of an oblique, supercritical shock. We introduce a new approach – the instability isolation method – to separate the steady-state shock physics from the physics of kinetic instabilities that may arise in the shock transition. An extension of the field–particle correlation technique enables us to separate the impact of the shock physics on ion energization from that due to the instabilities.

The instability isolation method divides the variation along the shock normal direction due to steady-state shock physics from the variation transverse to the shock normal due to kinetic instabilities. The electromagnetic fields are separated into the averaged fields and the fluctuating fields. A WFT is used to identify the dominant local wave vectors arising from the instability that lead to rippling of the shock front. We apply Faraday's law using this wave vector and its associated electric and magnetic field components to estimate the real frequency of the unstable mode. Comparing this frequency with the eigenfrequencies of the Vlasov–Maxwell linear dispersion relation demonstrates the Alfvénic nature of the fluctuations causing the shock ripple, likely arising from the AIC instability driven by a perpendicular ion temperature anisotropy caused by compression of the plasma through the shock, in agreement with the findings of earlier studies (Winske & Quest 1988; McKean et al. 1995; Burgess & Scholer 2007; Burgess et al. 2016).

The field–particle correlation technique uses single-point measurements of the electromagnetic fields and particle velocity distributions to obtain a velocity-space signature characteristic of a given particle energization mechanism and to compute the rate of energy transfer between the fields and the particles. We describe a specific procedure for the implementation of this technique with data from PIC based kinetic simulation codes. Further, we describe a new extension of the technique to use separately the transverse-plane-averaged fields to compute the steady-state shock energization through a newly defined averaged correlation (4.7) and the fluctuating fields to compute the instability energization through a newly defined fluctuating correlation (4.8).

The results of the averaged field–particle correlation, applied to measurements at the transition from the ramp to the foot of the shock, are presented in figure 10, with nine panels that present the energization due to the three components of the averaged electric field, each reduced along each of the three dimensions of velocity space. The correlation with the averaged motional electric field $\overline {C}_{E_y}(\boldsymbol {v})$, shown in panels (g)–(i), displays the characteristic velocity-space signature of shock-drift acceleration of reflected ions at an oblique shock, dominating the total ion energization in the shock foot.

Computation of the fluctuating correlation, shown in figure 11, shows that the fluctuations causing the shock ripple – although they produce a measurable impact on the surface of the shock as seen here in simulation and in in situ observations (Johlander et al. 2016, 2018) – have a negligible impact on the ion energization, with energization rates more than two orders of magnitude smaller than the steady-state shock physics. But the fluctuating correlation does produce a new, qualitatively distinct ‘tri-polar’ signature of ion energization, as seen in panel (d). By decomposing the transverse plane into sub-regions, we demonstrate that the tri-polar signature arises from variations of the fluctuating electric field and fluctuating velocity distribution across the integration volume used to compute the distribution. Therefore, the appearance of a tri-polar signature is indicative of energization by the electric field arising from an instability, leading to a qualitative means to identify particle energization by instabilities.

Integrating the averaged and fluctuating correlations over velocity space yields the net ion energization at the spatial point of analysis. This enables us to produce a profile of the energy transfer between the fields and the ions as a function of the normal distance through the shock, presented in figure 13. This shows that the instability experiences net wave growth in the ramp and overshoot regions of the shock, with a fraction of the wave energy propagating a short distance upstream into the foot of the shock before being entirely damped. The remaining instability-driven wave energy appears to propagate downstream with the bulk plasma flow, generating the oft-observed turbulence in the downstream plasma.

Here we have developed new tools to characterize more completely the energetics of collisionless shocks in phase space and have shown their value for analysing instabilities arising in shocks and the resulting impact of the instability-driven fluctuations on particle energization. Our success comparing linear kinetic theory predictions to analyse simulation results shows that the instability isolation method can be used quantitatively to characterize the nature of instability-driven fluctuations, even for mechanisms that transfer small amounts of energy in the system. Future work will apply this technique to shocks with higher Mach number in which instabilities are believed to have a more significant impact on the particle energization, in particular controlling the partitioning of upstream kinetic energy between ions and electrons, e.g. instances of non-adiabatic electron heating often observed in sufficiently high Mach number shocks.