Fine structure and motion of the bow shock, and particle energisation mechanisms inferred from MMS observations

This study presents new observations of fine structure and motion of the bow shock formed in the solar wind, upstream of the Earth’s magnetosphere. The NASA’s MMS mission has recorded during 2 hours eleven encounters with an oscillatory shock, which moves with the speed of 4–17 km/s and has thickness of 130 km, or an ion gyroradius. The shock is formed by steepening of 1 mHz magnetosonic wave, creating compressional magnetic field and plasma density structures, which initiate a chain of cross-field current-driven instabilities that heat solar wind ions by the stochastic ExB wave energisation 5 mechanism. The theoretical ion energisation limits are confirmed by observations. We have identified the ion acceleration mechanism operating at shocks and explained double beam structures in the velocity space. The nature of this mechanism has been revealed as a stochastic resonant acceleration (SRA). The results provide for the first time a consistent picture of a chain of plasma processes that generate collisionless shocks and are responsible for particles energisation.

energisation capacity implied by this (wave) mechanism is where v ⊥j is the initial perpendicular velocity of a particle with mass m j (j = e for electrons, j = p for protons, and j = i for general ions). This energy corresponds typically to 200 keV for protons and 1 keV for electrons in shocks measured by MMS.
The E × B acceleration is similar to the wave surfing (surfatron) mechanism (Katsouleas and Dawson, 1983;Ohsawa, 30 1985; Ucer and Shapiro, 2001;Kichigin, 2013), which applied to shocks requires a wide front of coherent waves Shapiro and Ucer, 2003). The energy for particles is provided by the convection electric field. In contrast, the E × B wave mechanism can work on intermittent bursty waves in any direction and the energy is taken from wave electric fields ∼ 10 − 100 mV m −1 , much larger than the convection field ∼ 3 mV m −1 . The present mechanism relies on a stochastic condition, which requires sufficiently strong gradients of the electric field to render particle motion chaotic and facilitate heating 35 (Cole, 1976;Karney, 1979;McChesney et al., 1987;Balikhin et al., 1993;Stasiewicz et al., 2000;Vranjes and Poedts, 2010).
The threshold for stochastic heating has recently been generalised to the form (Stasiewicz, 2020) χ j (t, r) = div(E ⊥ ) and applied to electron and ion heating observed at the bow shock (Stasiewicz and Eliasson, 2020a, b). Here, ω cj = q j B/m j is the angular cyclotron frequency of particle species with charge q j , N c is the number density of excess charges, N is the 40 plasma number density, V 2 Aj = B 2 /(µ 0 N m j ), c is the speed of light. The equivalent formula on the right side of Equation (2) implies that stochastic heating requires charge non-neutrality fraction larger than the ratio of the Alfvén speed, V Aj , to the speed of light squared. The particles are magnetised (adiabatic) for |χ j | < 1, demagnetised (subject to non-adiabatic heating) for |χ j | ≳ 1, and selectively accelerated to high perpendicular velocities when |χ j | ≫ 1.
Acceleration of ions in quasi-perpendicular shocks is performed by lower hybrid (LH) waves which have energisation 45 capacity for protons limited by waves phase speed (Stasiewicz and Eliasson, 2021) which will be shown to apply also in the analysed case. T e , T p are electron and proton temperatures in energy units.
In this paper we shall determine motion and thickness of shocks, and analyse ion distribution functions, magnetic and electric field turbulence measured at quasi-perpendicular shocks, as well as particle heating mechanisms implied by these 50 measurements.

Oscillatory bow shock
On January 3, 2020 the NASA Magnetospheric Multiscale Spacecraft (MMS)  were in solar wind at 13:40 UTC at the beginning of the data period shown in Figure 1. MMS entered the quasi-perpendicular shock #1 at 13:47 UTC at position (10.8, 13.8, -1.6), or R = 17.6 R E GSE (geocentric solar ecliptic), and then moved further earthward with the speed of 1.7 km s −1 . The variations of the dynamical solar wind pressure which was about 1 nPa caused oscillatory movements of the shock front with speed of 4 − 17 km s −1 and has led to eleven shock crossings within 2 hours on a distance of 2 R E , labeled in panel (a) with #1-11. The first crossing was caused by the outward motion of the shock front with speed of 15 km s −1 followed by an earthward motion of the shock #2 eight minutes later with the speed of 17 km s −1 in the spacecraft frame. The last shock crossing #11 was at 15:49 UTC, position (9.1, 12.7, -2.0) R = 15.8 R E , with outward speed 4 km/s. The four MMS spacecraft 60 had an average separation distance of 21 km. The velocity of the shock fronts has been determined with inter-spacecraft timing of the magnetic field measured by the fluxgate magnetometer (Russell et al., 2016). The motion is outward for all odd shock numbers and earthward for all even shocks.
The ion differential particle flux shown in panel (a) is measured by the Fast Plasma Investigation experiment (FPI) (Pollock et al., 2016) in energy range 10 eV -30 keV. It exhibits the solar wind beam centred around 700 eV which becomes thermalised in the shock regions while some ions are accelerated to a few keV. Over-plotted is the energisation capacity of lower hybrid waves K LH given by Equation (3). This equation exhibits good agreement with MMS measurements in all of ca 40 quasiperpendicular shocks analysed by the authors. Panel (b) shows perpendicular and parallel ion temperatures, T i⊥ , T i∥ , which confirm the known fact that ion heating in shocks is stronger in perpendicular direction.
The electron temperature in quasi-perpendicular shocks is isotropic, T e = T e⊥ ≈ T e∥ , and obeys a specific relation, which 70 has been found recently (Stasiewicz and Eliasson, 2020a) with α = 1/3. This relation, named quasi-adiabatic, predicts a dip of T e /B where B has a maximum. It has been derived under the assumption that the perpendicular energy gain, T e⊥ ∝ B, during compressions of the magnetic field is redistributed to the parallel energy component by scattering on waves, leading to the above temperature relation. At quasi-parallel shocks 75 we observe similar relation, but with α = 2/3, which has also a theoretical justification. The isotropisation of electrons is due to scattering on high-frequency oblique electrostatic waves with a parallel electric field component (Stasiewicz and Eliasson, 2020b).
Panel (c) shows compressions of the electron number density N e , and the magnetic field B occurring at shocks. To understand the process of nonlinear steepening of the magnetosonic waves that leads to the formation of perpendicular shocks we 80 perform multiresolution frequency decomposition of the measured magnetic field B from Figure 1c with orthogonal wavelets (Daubechies, 2009). The decomposition shown in panel (d) is exact, i.e., the sum of all components gives the original signal, and the orthogonality means that the time integral of the product of any different pair of the frequency dyads is zero. The numbered dyads in this stacked plot represent baselines (zero levels) for the signal B f /14 nT at the indicated frequency f . The residual 'dc' magnetic field is shown as a black line at the bottom with the same normalisation.

85
The decomposition suggests that the oscillatory behaviour of the shock and wave steepening process are related to the ∼1 mHz wave at the bottom, which triggers cascades of compressional waves extending to 1 Hz and above. The maximum amplitude of compressions is observed in the 0.5 Hz channel, which can be associated with ion cyclotron waves. The proton cyclotron frequency is 0.1 Hz in the solar wind regions, but it goes up to 0.6 Hz in shock compressions. The Alfvén Mach number for the plasma flow is M A = V i /V A ≈ 7 in solar wind regions, and plasma beta is β e ≈ 1, β i ≈ 2. Some additional 90 diagnostic parameters for these shocks can be found elsewhere (Stasiewicz and Eliasson, 2020a).

Burst data analysis
In this section we focus our analysis on high resolution burst data measured during time 14:31:36-14:32:22 UTC which contains shock #4 of Figure 1a. ion measurements are transformed to a cartesian coordinate system in which thex axis is along the ExB direction, theẑ axis is along the magnetic field, and theŷ =ẑ ×x is along the electric field, forming the ExB reference system in velocity space We have used the electron convection electric field E = −V e × B to construct these coordinates and apply a different notation to distinguish between the ExB convection and wave induced acceleration E × B . Colour spectrograms  Figure 3. The shock ramp, identified with the B and T e profiles in panels (d) and (e), is within blue vertical lines.
Electrons are mostly in quasi-adiabatic regime, |χ e | < 1, which means that the temperature shown in panel (d) follows the quasi-adiabatic relation (4), T e ∝ B 1−α at shocks (Stasiewicz and Eliasson, 2020a). On the other hand, ions are in strongly plasma, whereas Figure 2d shows the opposite. Secondary beams are produced by the E × B acceleration, as we will show further in the text.
In Figure 2f we show the time-frequency spectrogram of the field measured by the electric field double probes instruments (Ergun et al., 2016;Lindqvist et al., 2016) with sampling rate 8192 s −1 . Over-plotted is the lower hybrid frequency f lh and the proton cyclotron frequency f cp . Lower hybrid drift waves, can be identified in the frequency range f cp − f lh . This frequency range contains also ion whistler waves, which could originate from mode conversion of lower hybrid waves on density striations (Rosenberg and Gekelman, 2001;Eliasson and Papadopoulos, 2008;Camporeale et al., 2012).
The vertical striations seen in the spectrogram (f) represent cascades of instabilities: LHD → MTS → ECD extending form f cp up to nf ce in a few kHz range. The presence of ECD instability at shocks has been reported in several papers (Wilson III et al., 2010;Breneman et al., 2013;Stasiewicz, 2020). Lower hybrid waves generated in the shock propagate upstream in panel 120 (f) and appear to be associated with particles in panel (b).
The velocity of the shock in Figure  The shock comprising the ramp and foot would have thickness 2r p embracing the whole proton orbit, which can be inferred from data presented in Figure 2, and in particular from the ion temperature in panel (d). These values agree with many other estimates of shocks thickness and motion published by other authors. However, shock thickness scalings based on ion inertial 130 length, or the hybrid gyroradius (r p r e ) 1/2 are not supported by measurements.

Ion distribution functions measured at shocks
We shall now inspect the ion distribution functions shown in Figure 3 in columns A, B, C, D and E, which correspond to events marked in Figure 2. Each picture shows 2-dimensional reduced distribution function in the reference system (v ExB , v E , v B ). Event A shows partly thermalised ions in the magnetosheath with some remaining non-gyrotropic features. The crescentlike structure in distribution A1 is characteristic for E × B acceleration which will be explained further. Event B shows ion distribution downstream of the shock peak, and event C on the upstream side of the peak. Event D is in the middle of the shock ramp, and E in the foot of the shock. All distributions are strongly non-gyrotropic, some with separated beams. Similar distributions with double beam structures have been reported by many authors (Paschmann et al., 1982;Gosling and Thomsen, 1985;Fuselier, 1994;Mazelle et al., 2003;Kucharek et al., 2004;Wilson III, 2016;Johlander et al., 2016) and interpreted usually in terms of specularly reflected ions, non-specularly reflected ions, gyrating ions, gyro-phase-bunched ions, or simply shock reflected ions.
Particularly puzzling are multiple peaks in the perpendicular plane (first row). Ions reflected from magnetic barriers could 145 acquire a different parallel velocity component V ∥ , but they are in the same electric field so they should have the same V ⊥ ≈ V ExB velocity component as the original solar wind beam, with possible modifications by temperature dependent gradient drifts. However, we observe secondary beams in all directions in the perpendicular plane, with similar parallel velocities, which appears to be in odds with standard plasma physics.

The E × B wave energisation mechanism 150
In this section we shall argue that the presented observations are consistent with the E × B acceleration (Stasiewicz and Eliasson, 2021). First, we should distinguish between the convection ExB drift ∼ 500 km s −1 used to establish the coordinate system, and the wave electric drift V ExB = E × B/B 2 at higher frequencies. For wave amplitudes of ∼ 50 mV m −1 in a magnetic field of 7 nT, the later is V ExB ∼ 7, 000 km s −1 , corresponding to the gyration speed of a 250 keV proton, which can explain acceleration of ions in quasi-parallel shocks (Stasiewicz and Eliasson, 2021;Stasiewicz et al., 2021;Stasiewicz and Kłos, 155 2022).
Usually, particles do not obey the electric drift in waves with frequencies higher than the gyrofrequency or wavelengths smaller than the gyroradius. In such situations, the effects average to zero over the wave period or wavelength. However, when the electric field gradient ∂ x E x = k x E x , for electrostatic waves with wavevector k x , exceeds the stochastic condition in Equation (2), the particle can be accelerated to the value in Equation (1) within a fraction of the gyroperiod.

160
The mechanism of E × B acceleration by electrostatic waves can be studied with the Lorentz equation (Stasiewicz andEliasson, 2020a, b, 2021). The previous model is generalised here for waves propagating in arbitrary direction in the perpendicular plane to the magnetic field B 0 = (0, 0, B 0 ). The position r and velocity v of an ion with mass m and charge q are determined by the equation mdv/dt = q(E + v × B 0 ) together with dr/dt = v. We assume that convection electric field E y convects plasma into electrostatic wave E w sin(ω D t − k · r) propagating in the (x, y) plane at angle α to the x-direction, with the Doppler shifted frequency ω D in the observer's frame. By using dimensionless variables with time normalized by ω −1 c , space by k −1 and velocity by ω c /k with ω c = qB 0 /m p being the angular ion cyclotron frequency, the normalised equations of motion for a test ion in a stationary (shock related) frame are upstream, α = 180 • , active during time 0.4f −1 cp . The particle is accelerated more than 4 times the initial gyration energy as shown in the lower panel. Acceleration increases uy while ux remains constant, which corresponds to the situation seen in Figure 3E1.
Here, Φ = Ω D t − x cos α − y sin α is the wave phase with the Doppler shifted angular frequency Ω D = ω D /ω c = Ω + χ d cos α with respect to Ω in the plasma frame. The normalised amplitudes of the E × B drift and the convection drift are, respectively Please note that χ w = χ p represents here the stochastic wave parameter given by Equation (2). By setting χ d = 0 we obtain 175 equations in the plasma frame of reference.
The most efficient energisation occurs on the acceleration lane (Stasiewicz and Eliasson, 2021), which corresponds to u 0 = Ω that matches particle velocity with phase speed of waves. The initial conditions for the presented here solutions are chosen in such a way that the gyration velocity v 0 at t = 0 is aligned with the (k x , k y , 0) vector, or alternatively with the phase velocity of waves, so that in the plasma frame we have where u 0 = v 0 k/ω c = kr c , and r c = v 0 /ω c is the gyroradius.
Generally the equations have chaotic solutions, because for χ w > 1 the solutions are very sensitive for initial conditions and have positive Lyapunov exponent (Balikhin et al., 1993;Stasiewicz et al., 2000). They are representative of deterministic chaos. The parameters of these equations are: Ω, χ w , χ d , u 0 , α which can be varied to fit particular physical conditions. The frequency Ω = 25 is below the lower hybrid frequency and the ratio χ w /χ d = E w /E y = 8 is realistic.  Figure 5. The same as in Figure 4, but without the convection electric field, χ d = 0. The waves are active a longer time during 0.75f −1 cp and the propagation direction of waves is reversed, α = 0. Acceleration is the same, but in the negative uy direction. At time t = −1 the proton in Figure 4 has initial gyration energy K 0 = u 2 0 and is drifting earthward with the convection speed u x = χ d = 5. During time t = 0 − 0.4 it experiences a burst of lower hybrid waves with amplitude χ w = 40. We see that 190 the gyration energy K = u 2 x + u 2 y has increased more than 4 times after a couple of wave periods. The acceleration is in the u y direction, while u x is constant, which could correspond to Figure 3E1.
The convection electric field E y is an essential element in the shock surfing (surfatron or SSA) acceleration by waves in front of shocks Shapiro and Ucer, 2003) and in shock drift acceleration (SDA) models based on the magnetic gradient drift. The situation observed in Figure 2b, where particles are accelerated along the convection electric field 195 E y suggests that we may have the surfatron case here.  Figure 2b we see ions accelerated in the negative v E direction, which is most likely due to downstream propagating waves as in case of Figure 5. For waves propagating at α = ±90 • the acceleration is in the v x , or equivalently v ExB direction, which is also observed in measurements. By changing the wave propagation angle α and the amplitude of waves χ w we can reproduce any secondary ion peak which can be found in Figure 3, row (1). A free gyration 205 after acceleration would produce crescent-like structures seen in most distributions.
It can also be seen that the duration of wave activity is not an essential factor. In Figure 6 we show the electric field E x seen by the particle along the trajectory made in Figure 5. The energisation occurs only during short coherence time after t = 0 by means of v x E x > 0. The work done by the electric field on the v x component is transferred to the v y component by the Lorentz force v x B z . The mechanism is inherently bursty and works only during short coherence times of a few wave periods. After 210 decoherence, the waves do not affect particles anymore, as can be seen in Figures 4-6. We can conclude that stochastic particle energisation by waves is performed in a sequence of coherent resonant interactions. This leads to the concept of stochastic resonant acceleration (SRA) as a complementary description of the E × B wave mechanism.
The stochastic condition in Equation (2) is necessary for energisation of particles. When χ w < 2 no significant acceleration can be produced by Equations (5)-(9), irrespectively of the values of other parameters. The convection electric field E y plays 215 no role in the E × B energisation, which could be anticipated. Indeed, transformation between the plasma frame of reference where E y = 0, χ d = 0 and the shock fixed frame with the convection electric field cannot involve E y in particle energisation because both are equivalent inertial systems.
Secondary beams in the perpendicular plane, such as seen in Figures 2b and 3 are commonly observed in front of quasiperpendicular shocks and have been usually described as shock-reflected ions. In shock reflection scenarios it has been usually 220 assumed that the E x field which makes the cross-shock potential is responsible also for the reflection. Contrary to this popular belief, strong E x field does not reflect ions upstream, but accelerates them in the y-direction through the E × B mechanism in a stochastic resonant way.

Conclusions
This research provides confirmation of the plasma heating/acceleration scenario in shocks outlined in earlier publications 225 (Stasiewicz, 2020;Stasiewicz andEliasson, 2020a, b, 2021;Stasiewicz et al., 2021;Stasiewicz and Kłos, 2022). Shocks oscillatory movements and development of compressions appear to be related to 1 mHz magnetosonic wave in Figure 1d which leads to cascades of instabilities, including rapid growth of ion cyclotron waves as seen in magnetic waveforms in Figure 1d and in the electric spectrogram in Figure 2f. The instability progresses to waves around the lower hybrid frequency (∼ 10 Hz) and further up to a few kHz, generating a cascade of instabilities LHD → MTS → ECD mentioned in the Introduction. The significance of these cross-field current-driven instabilities for heating of the solar wind plasma has been advocated earlier by many authors (Forslund et al., 1972;Lashmore-Davies and Martin, 1973;Yamada and Owens, 1977;Wu et al., 1983;Zhou et al., 1983;Winske et al., 1985;Drake et al., 1983;Gary, 1993;Daughton, 2003;Muschietti and Lembége, 2017).
Using only the fundamental Lorentz equation we have identified the E × B wave mechanism which explains how lower hybrid waves accelerate ions to velocities of 800 km s −1 , as can be seen in Figure 2. We have shown that stochastic particle 235 energisation by waves occurs in a series of coherent resonant interactions. The nature of this mechanism can be described as a stochastic resonant acceleration (SRA). The model is also capable of explaining multi-beam ion distributions measured at shocks and shown in Figure 3. The observed values for acceleration agree with the limits given by Equations (3) and (1), computed for lower hybrid waves (Stasiewicz andEliasson, 2020a, b, 2021).
Energisation of particles depends on interaction time with waves. Particles are convected rapidly across perpendicular shocks 240 with thickness of 100 km, but can spend considerably longer times in a spatially extended turbulence (a few R E ) of quasiparallel shocks. The short interaction time in quasi-perpendicular shocks limits the ion acceleration to a few keV or velocities v p < 1, 000 km s −1 as can be seen in Figures 1a, 2a,b,c and 3. Waves involved in acceleration are in the frequency range f cp − f lh . The longer interaction time with higher amplitude waves at higher frequencies in quasi-parallel shocks makes it possible to accelerate protons to velocities v p ∼ 7, 000 km s −1 that correspond to energies of 250 keV (Stasiewicz et al., 2021;245 Stasiewicz and Kłos, 2022).
Using exceptional quality, multipoint measurements of MMS we have made exact determinations of the shock ramp thickness which is about 100 km, while the ramp and foot combined have thickness of 2 gyroradii that embraces the whole ion cyclotron orbit, or 200 km. We have also pointed out that high perpendicular ion temperatures measured in front of shocks are mainly the result of secondary beams produced by the wave acceleration process.

250
Code and data availability. The data underlying this article are available to the public through the MMS Science Data Center at the Laboratory for Atmospheric and Space Physics (LASP), University of Colorado, Boulder: https://lasp.colorado.edu/mms/sdc/public. The data have been processed with the IRFU-Matlab analysis package available at https://github.com/irfu/irfu-matlab.
Author contributions. KS made data analysis and wrote the paper. ZK contributed to the text.
Competing interests. The authors declare that they have no conflict of interest.