Electron Acceleration at Earth's Bow Shock Due to Stochastic Shock Drift Acceleration

We use the Magnetospheric Multiscale mission (MMS) to study electron acceleration at Earth's quasi‐perpendicular bow shock to address the long‐standing electron injection problem. The observations are compared to the predictions of the stochastic shock drift acceleration (SSDA) theory. Recent studies based on SSDA predict electron distribution being a power law with a cutoff energy that scales with upstream parameters. This scaling law has been successfully tested for a single Earth's bow shock crossing by MMS. Here we extend this study and test the prediction of the scaling law for seven MMS Earth's bow shock crossings with different upstream parameters. A goodness‐of‐fit test shows good agreement between observations and SSDA theoretical predictions, thus supporting SSDA as one of the most promising candidates for solving the electron injection problem.


Introduction
Collisionless shock waves are ubiquitous throughout our universe and an essential source of cosmic ray generation.It has been shown by studying radio and X-ray synchrotron emissions that supernova remnant (SNR) shocks can produce ultra-relativistic electrons (Reynolds, 2008).To explain these observations, the diffusive shock acceleration (DSA) mechanism has become widely accepted as the standard theory for electron (and ion) acceleration up to ultra-relativistic cosmic ray energies (Blandford & Eichler, 1987;Drury, 1983).However, this theory requires a seed population of pre-accelerated electrons up to mildly relativistic energies to efficiently provide further energization (Amano & Hoshino, 2022;Balogh & Treumann, 2013a;Levinson, 1992;Longair, 1981).This is known as the electron injection problem, and much effort has been put down during the recent decades to resolve it.
The stochastic shock drift acceleration (SSDA) theory is currently one of the most promising candidates for solving the long-standing electron injection problem (Amano & Hoshino, 2022).SSDA builds on the common shock drift acceleration (SDA) theory, described in Wu (1984), by introducing an additional stochastic component where accelerated electrons are strongly pitch-angle scattered in the shock transition layer (STL).The scattering agent allows electrons to undergo the SDA process longer and therefore gain more energy.Theoretical studies show that this theory can produce a power law spectrum for suprathermal electrons (Katou & Amano, 2019), supported by observational evidence in a single case (Amano et al., 2020).injection efficiency compared to for example, SNR shocks (Dresing et al., 2016).The less efficient heliospheric shocks still provide excellent opportunities for testing the SSDA theory based on in situ observations.Recent work by Dresing et al. (2022) shows that interplanetary (coronal) shock waves can energize electrons up to ∼1 MeV.Which is around 10 times higher than typical electron energies observed at Earth's bow shock.
The important parameter controlling the injection is the SSDA cut-off energy, E c .By combining the theory described in Katou and Amano (2019) and recent in situ observations at Earth's bow shock by Amano et al. (2020), a theoretical scaling law for the SSDA cut-off energy is obtained where u 0 is the upstream solar wind speed normal to the shock surface in the shock frame, θ Bn the angle between the local shock normal and upstream magnetic field, D μμ the pitch-angle diffusion coefficient (scattering rate) and Ω ce the electron cyclotron frequency.At large-scale SNR shocks (having u 0 ≈ 3,000 km/s), SSDA generally produces a cut-off energy above the injection threshold, while at the weaker heliospheric shocks (u 0 ≈ 400 km/s), SSDA generally produces a cut-off energy below the injection threshold.Hence, DSA is not expected at heliospheric shocks unless θ Bn is close to 90°.Despite the large amount of shock data available over the years, only one case study has been performed at Earth's bow shock, providing observational evidence of the SSDA theory (Amano et al., 2020).Therefore, it is crucial to extend this study to multiple shock observations under different solar wind conditions to establish the SSDA validity and investigate its relationship with upstream driving conditions.
In this paper, we present new observational evidence supporting the SSDA theory.We use the Magnetospheric Multiscale (MMS) spacecraft to study seven shock crossings of Earth's bow shock.The cut-off energy, E c , and pitch-angle diffusion rate, D μμ , are determined for each shock crossing and used to test the scaling law predicted by Amano et al. (2020).

Methodology
Seven shocks are selected from the shock database in Lalti et al. (2022).The shocks are chosen such that they cover different ranges in the upstream parameters; shock angle, θ Bn , and Mach number, M A .All seven shocks display accelerated electrons in the vicinity of the shock ramp/foot.In all of them, energetic electrons are observed with varying intensity and energies up to or above 60 keV.
We use data obtained from the Fast Plasma Investigation (FPI) (Pollock et al., 2016) and Fly's Eye Electron Proton Spectrometer (FEEPS) (Blake et al., 2016) for the distribution function, density, and solar wind bulk flow.Magnetic field data is obtained from the Flux Gate Magnetometer (Russell et al., 2016) and Search Coil Magnetometer (Le Contel et al., 2016).
To test the SSDA scaling law for the seven shock crossings, the four parameters in Equation 1 need to be obtained with as great precision as possible.We use methods described in Paschmann and Daly (2000) to determine u 0 and θ Bn where the local shock normal is determined as the average of the three mixed methods also described in Paschmann and Daly (2000).These methods require reliable upstream and downstream plasma density, velocity, and magnetic field measurements.These quantities are obtained by averaging their values over multiple pairs of upstream and downstream time intervals of varying lengths (45-60 s), similar to that described in Trotta et al. (2022).To ensure the data is taken sufficiently far away from the STL, the interval pairs are chosen outside of the selected burst-data time interval of each shock.
To determine the cut-off energy, E c , and normalized pitch-angle scattering rate, D μμ /Ω ce , we follow the methods outlined in Amano et al. (2020).To increase the count statistics, the data is averaged over MMS1-3.MMS4 electron data was either missing or highly contaminated and was therefore omitted.Figure 1 shows two of the seven shock crossings investigated in this study.According to the SSDA theory, pitch-angle scattering of electrons in the STL causes an exponential increase in energetic electron phase space density (PSD).Therefore, the scattering rate can be determined from the e-folding time, τ(E), obtained by fitting the energetic electron intensity profiles across the shock.The time intervals marked with yellow in Figure 1 indicate the intervals used when calculating the pitch-angle scattering rate.The power law in the PSD is predicted to form when the PSD starts to saturate.Therefore, the time intervals used for the estimates of the power law and cutoff energy are chosen directly downstream of the region used for scattering rate estimates and are marked with green in Figure 1.
The PSD profiles, measured by MMS, are fitted to an analytical expression f(E) ∝ E p exp( E/E c ), where the spectral index (p) and cut-off (E c ) are used as fitting parameters.The fitting expression describes a power-law with exponential cut-off and the fitting data is restricted to data points above 1 keV which is well above the thermal energy.The resulting fits are shown in Figure 2a and the determined cut-off energies in Table 1.
The pitch-angle scattering (diffusion) rates as a function of energy are shown in Figure 2b along with their respective theoretical threshold.Under the assumption of wave amplitudes much smaller than the background magnetic field, quasi-linear theory can be used to obtain the scattering rate and theoretical threshold Figure 1.Shock crossings number 2 and 3 (see Table 1) showing Magnetospheric Multiscale mission burst data of (a) magnetic fields, (b) high energy electron spectra (c) low energy electron spectra (d) ion spectra (e) ion velocities (f) ion densities.The yellow-marked time intervals are used to determine the pitch angle diffusion rates (see Figure 2b) from the phase space density-increase at the shocks following the method outlined in Amano et al. ( 2020).The green-marked time intervals are used for the determination of the cut-off energies. and respectively where V sh is the shock velocity, E sh = m e u 2 0 /2 cos 2 θ Bn and η is a numerical factor between 0 and 1 determined by the fractional thickness of the STL that characterizes the fraction of the time the accelerated electrons are observed.See Katou and Amano (2019) and Amano et al. (2020) for a detailed description.The theory states that a scattering rate above the threshold can efficiently trap electrons within the acceleration region, resulting in a power-law distribution.The distribution cut-off is expected to occur as the scattering rate goes below the threshold.Therefore, the pitch-angle scattering rate used in the scaling law is the one evaluated at the cut-off and is determined as the average of the data points having their standard deviation below the theoretical threshold.These points are marked with their respective error bars in Figure 2b and the determined scattering rate is documented for each shock in Table 1.
Having determined the four parameters constituting the SSDA scaling law depicted in Equation 1, we can define the ratio of the observed cutoff energy derived from the distribution function fits versus the theoretically predicted cutoff energy derived from Equation 1 in which we use the observed parameter values ). (4) This way of writing the scaling law enables us to perform a χ 2 -test providing a goodness-of-fit of the theory (Cowan, 1998) where σ X,k is the standard deviation of the observed X-value, X k , for each (k) shock.The X-parameter as defined in Equation 4 is a function of four individual parameters.Similar to the scaling law, it has a strong dependence on the shock angle θ Bn , and small changes in the angle close to 90°can give large deviations.Also, the scattering rate error (standard deviation) has the same order of magnitude for many shocks (see Table 1).Therefore, to estimate X k and its uncertainties σ X,k we perform Monte Carlo simulations, generating sample distributions of X k using the measured quantities on the right-hand side of Equation 4. We assume that each of the four measured quantities can be described by a normal distribution constructed from their observationally determined value (mean) and standard error.Two example distributions are illustrated in Figures 2c and 2d.
The X k -distributions are then generated by drawing random values from these four normal distributions and utilizing Equation 4. The X k and σ X k used in Equation 5 are then obtained as the mean and standard error of the X kdistributions, see the red and green dashed vertical lines in Figures 2c and 2d.

Results
Table 1 contains the time stamp, Alfvén Mach number, scaling law parameters, and calculated X-values for all seven shock crossings considered in this study.Figure 3a visually represents the observed X-values and their standard error.The theoretical prediction (X = 1) is marked by the red horizontal line.Using Equation 5, we obtain a χ 2 -value (goodness-of-fit) of χ 2 = 0.55, suggesting good agreement between our observations and the theory (Cowan, 1998).
Figure 3b shows the normalized pitch angle scattering rate versus Alfvénic Mach number and θ Bn .In addition to the seven shocks studied, we include the shock studied in Amano et al. (2020).We observe a clear decreasing trend in the scattering rate for increasing Alfvén Mach number and/or θ Bn being closer to 90°.The reasons for these dependencies are discussed in Section 4.
The cut-off energy, E c , is observed to vary between 9 and 24 keV.The spectral index (p) does not vary much among the shocks.On average, we find a spectral index p = 3.4 ± 0.2 consistent with previous work by Oka et al. (2018) and Gosling et al. (1989) at Earth's bow shock.

Discussion
Out of the seven shock crossings investigated, we identified two types of shock crossings.They are distinguished by the amount of energetic electron flux registered at the high-energy channels (47-515 keV).Figure 1 illustrates two examples of the different shock types where energetic electron flux can be seen in the panels (b1) and (b2).Shocks number 1 and 2 display similar features to the case study shown in Amano et al. (2020), with enhanced energetic electron flux observed over a large portion of the STL and the flux being well above the noise levels of the energy channels.On the contrary, shocks number 3-7 have their enhanced electron flux closer to the noise level and much more localized to the vicinity of the shock ramp.The reason for this difference is currently unknown.Shocks 1 and 2 appear to have the greatest Mach numbers but it's not a significant difference.However, one major difference between the two groups (1-2 vs. 3-7) is a much higher solar wind speed for the first group (1-2) compared to all other shocks, see Table 1.
The normalized χ 2 -value obtained by Equation 5 lies safely below 1 and corresponds to a p-value of 0.8.Thus suggesting that our seven shocks follow the SSDA scaling law quite well (Cowan, 1998).This is illustrated in Figure 3a where all except one shock have the theoretical X-value (X = 1) within their standard error limit (shock number 5 having it on the edge).
An important part of this study is the estimation of the standard errors.As described in the above section, the uncertainty in u 0 and θ Bn are estimated by performing multiple calculations based on different intervals while the standard error on the cut-off energies, E c , is determined from the data-fits in Figure 2a.
The scattering rate error is the most difficult to estimate.It is calculated from several quantities all carrying some degree of uncertainty, see Equation 2. However, we consider the most dominating errors to originate from the intensity profile fits (see Amano et al. (2020)), determining the e-folding time τ(E), and the uncertainty in the shock velocity (V sh ) determination (Paschmann & Daly, 2000).Therefore, the standard errors shown in Figure 2b are obtained assuming contribution from only these two parameters.Another factor complicating the pitch-angle scattering rate calculations is the effects of the data compression (Barrie et al., 2017) present at the highest energy channel(s) of the FPI-instrument (Pollock et al., 2016).Data compression effects typically occur when the plasma density (count rate (Barrie et al., 2017)) is large which is often where the SSDA power-law is expected to form.Consequently, these data points cannot be fully trusted and must be removed from the fitting data.Unfortunately, these energy levels (around 10-27 keV) are where the cut-off typically appears and are therefore responsible for the relatively large errors on the pitch-angle scattering rate, see Table 1.
When performing the Monte Carlo simulations, all four scaling law parameters are assumed to be normally distributed with their estimated error as standard deviation.Due to this assumption, a small fraction of the drawn Monte Carlo samples might give an unphysical result and are therefore removed.However, this removal is estimated to have a small impact on the X-distributions (Figures 2c and 2d) for all shocks.When calculating the χ 2 , we use the average and standard errors determined from the X k -distributions instead of the obtained X-values using the formula in Equation 4. Note that these two values differs slightly.We consider determining the mean and standard error from the sampled X k -distributions to be more accurate rather than obtaining them from the scaling law formula, Equation 4. Acquiring the standard deviations for the X-values, obtained directly from the scaling law, requires standard deviations of the four measured quantities (E c , u 0 , θ Bn , and D μμ ) to be small, relative to their measured values, which is not the case for all our shocks.
The pitch-angle scattering/diffusion rate is essential for the SSDA theory and determines the cut-off energy for energetic electrons.Figure 3b shows that smaller pitch-angle scattering rates are observed at higher Alfvénic Mach numbers M A and/or θ Bn closer to 90°.Due to the relatively small data set (8 shocks), observational bias cannot be completely ruled out.Note that the threshold wave power (required to sustain the acceleration) has a strong dependence on the Alfvén Mach number and shock obliquity, ∝(M A / cosθ Bn ) 2 (Katou & Amano, 2019).Previous work by Katou and Amano (2019) shows that this threshold decreases with increasing M A / cos θ Bn .In other words, for a given strength of wave power, the efficiency of SSDA extends for larger M A / cos θ Bn and therefore a smaller scattering rate is required to sustain the process.The data in Figure 3b appears to support this claim and the theoretical result stated in Katou and Amano (2019).Also, an important pitch-angle scattering agent is the presence of upstream whistler waves (Amano et al., 2020;Amano & Hoshino, 2022), previous work on whistler waves predicts the importance of a critical whistler Mach number(s) (Balogh & Treumann, 2013b;Krasnoselskikh et al., 2002;Lalti et al., 2020;Oka et al., 2017).For M A < M wh , whistler waves can propagate upstream of the quasi-perpendicular shocks to form standing wavefronts.All our shocks have an Alfvénic Mach number larger than (or similar to) its respective (linear) whistler Mach number, M wh = 0.5 ̅̅̅̅̅̅̅̅̅̅̅̅ ̅ m i / m e √ |cos θ Bn |.Whether this whistler-super-criticality can affect the electron pitch angle scattering rate at the shock and contribute to the decreasing trend observed in Figure 3b is left for future investigation.
The data presented in this work supports the SSDA mechanism and cements it as one of the most promising candidates for a solution to the injection problem.However, some open questions remain to be answered regarding suprathermal electron acceleration up to mildly relativistic energies.Recent observations show that other electron acceleration processes can contribute to the solution of the injection problem by combining different acceleration mechanisms (Lindberg et al., 2023).Also, the important role of the pitch-angle scattering agent is still poorly understood.The strong emissions of whistler waves at Earth's bow shock can sufficiently pitch-angle scatter electrons with energies ≥1 keV via cyclotron resonance (Amano et al., 2020).However, at frequencies above 100 Hz, the whistler wave power is generally not high enough at Earth's bow shock to sufficiently scatter electrons with energies much less than 1 keV, and hence other wave-particle interaction(s) is needed to explain the pitch angle scattering of electrons at thermal energies (10-100 eV).Electrostatic waves are currently under investigation and could be the source of scattering at thermal energies (Kamaletdinov et al., 2022).

Figure 2 .
Figure 2. (a) Magnetospheric Multiscale mission electron distribution functions with fits.(b) Pitch angle scattering rates and theoretical thresholds (c-d) X-parameter distributions obtained from Monte Carlo simulations using Equation 4 for shock numbers 2 and 3 respectively.The red dashed line indicates the average X-value used in χ 2 -fit and the green dashed lines indicate the standard deviation.

Figure 3 .
Figure 3. (a) Observed X-parameter, defined in Equation 4. The theoretical prediction obtained using the stochastic shock drift acceleration scaling law in Equation 1 is X = 1, red line.(b) Scatter plot of pitch angle scattering rate versus Alfvén Mach number and shock angle, θ Bn .

Table 1
The Calculated Parameters of the Seven Shocks Studied in This Paper