Ion Acceleration Efficiency at the Earth's Bow Shock: Observations and Simulation Results

In this study, we use data from the four MMS spacecraft and primarily from MMS2. Ion moments and distribution functions are from the Fast Plasma Investigation Dual Ion Spectrometer (FPI-DIS) (Pollock et al. 2016) onboard MMS. FPI-DIS can sample the full ion distribution function every 150 ms. The energy range of the ion measurements used here is from 2–18 keV, and 28 keV during different times. Magnetic field data are from the FluxGate Magnetometer instrument, which provides magnetic field vector measurements at a resolution of 128 Hz (Russell et al. 2016). All upstream solar-wind parameters are from a set of spacecraft positioned upstream of MMS at the first Lagrange point. The solar-wind data are time shifted to the bow shock and provided by the OMNI database with a time resolution of 1 minute (King & Papitashvili 2005).

We perform a statistical study using encounters with the bow shock by the MMS spacecraft. The events are selected between 2017 October and 2018 April during a period when MMS crosses the bow shock every orbit. The orbit of MMS during the time period starts out crossing the bow shock on the dusk-side flank and ends in the dawn-side flank while covering the subsolar region in between. This time interval is also close in time to solar minimum, which reduces the risk of energetic ions from the Sun or geomagnetic storms. The MMS events are selected from the Scientist In The Loop (SITL) reports, which are documents that prioritize burst data downlink for MMS. We pick out events from the SITL reports where both the upstream (unshocked) and downstream (shocked) plasmas were clearly observed in the same time interval. We are interested in how upstream shock parameters influence ion acceleration at the bow shock. Therefore, we only select events with stable upstream OMNI conditions for 10 minutes around the event. We only use events where the magnetic field direction standard deviation is <8° and the maximum deviation is <25°. We also only consider events where the solar-wind Alfvén Mach number has a standard deviation <3 and a maximum deviation <6. These standard deviations can be seen as a rough estimate of error bars of these parameters. We further disregard eight events where the OMNI data clearly disagree with the upstream solar wind observed by MMS. Lastly, we only include events where FPI-DIS covers up to >20 times the solar-wind energy. In the end we end up with 154 bow-shock crossing events by MMS with data in burst resolution. The positions of all used shock crossings in Geocentric Ecliptic Coordinates (GSE) are shown in Figures 1(a)–(c). The shock angle determination is described below. The data cover the entire day side of the bow shock and a wide variety of solar-wind conditions.

Zoom In Zoom Out Reset image size

In order to determine the Mach numbers and shock angle of the bow shock for all events we need to determine the local bow-shock normal vector $\hat{n}$ for each event. It is in principle possible to determine $\hat{n}$ from local measurements of the plasma (e.g., Abraham-Shrauner 1972). However, this can be deceptive when energetic ions are present since the upstream plasma may have large upstream B fluctuations that make a local determination of $\hat{n}$ and θ_Bn highly inaccurate (e.g., Battarbee et al. 2020). Obtaining $\hat{n}$ from timing analysis (Schwartz 1998) is not possible because of the relatively small size of the spacecraft tetrahedron formation, which is ∼20 km in most of the cases studied here. Instead we determine the local bow-shock normal by using the bow-shock model by Farris et al. (1991). We fit the bow-shock model through the position of MMS at the time of the event and obtain $\hat{n}$ from this fitted model (Schwartz 1998). This method is possible since all MMS events are in direct contact with the bow shock and the shock parameters we obtain are the average upstream conditions and not local parameters that are sensitive to upstream fluctuations and small-scale structures within the shock. When we have obtained $\hat{n}$, we can calculate the shock-normal angle θ_Bn between the model shock normal and the solar-wind upstream magnetic field obtained from the OMNI database. Another parameter of interest is how far out on the flank of the bow shock MMS is positioned. We define the angle ϕ as the angle between the MMS position and the Sun–Earth line seen from the center of the Earth at the time of the event. Then ϕ = 0 at the subsolar point of the bow shock and large values of ϕ means that MMS is on either the dawn- or dusk-side flank.

We also calculate the Mach number of the shock from the solar-wind parameters provided by OMNI. We calculate both the Alfvén Mach, M_A , and the fast magnetosonic Mach, M_ms, numbers in the shock frame assuming a zero shock speed relative to the spacecraft. The Alfvén Mach number of the shock is

$$\begin{eqnarray}&&{M}_{A}=\displaystyle \frac{{{\boldsymbol{V}}}_{u}\cdot \hat{n}}{{v}_{A}}=\displaystyle \frac{{V}_{u}}{{v}_{A}}\cos {\theta }_{{Vn}},\end{eqnarray} \tag{ 1 }$$

where V _u is the upstream flow velocity, v_A is the upstream Alfvén speed, and θ_Vn is the angle between $\hat{n}$ and V _u . We also calculate the fast magnetosonic Mach number

$$\begin{eqnarray}&&{M}_{\mathrm{ms}}=\displaystyle \frac{{{\boldsymbol{V}}}_{u}\cdot \hat{n}}{{v}_{\mathrm{ms}}({\theta }_{\mathrm{Bn}})},\end{eqnarray} \tag{ 2 }$$

where the v_ms(θ) is the fast magnetosonic speed for a propagation angle θ to B . We assume that the upstream electron temperature is the same as the ion temperature T_i so that the sound speed ${c}_{s}^{2}=4{k}_{B}{T}_{i}/{m}_{i}$, where m_i is the proton mass and k_B is the Boltzmann constant.

The shock angle and Mach number distributions of the MMS shock crossings are shown in Figures 1(d)–(e). Compared to the typical solar wind, quasiparallel shock geometries are somewhat overrepresented, likely due to selection bias in the SITL reports. The Mach numbers have rather typical solar-wind values with M_A ∼ 5–15 and M_ms ∼ 3–8, although with some limitations on high Mach number by the selection condition that FPI-DIS covers energies greater than 20 times the solar-wind energy. This data set has a good coverage of quasiparallel and quasiperpendicular shock regions as well as low to high Mach numbers.

Next, we quantify the presence of high-energy ions at the shock in order to compare this value to the derived shock parameters. We adopt the same definition as Caprioli & Spitkovsky (2014a) of ion acceleration efficiency ε(E₀) to be the fraction of energy density of ions with energy E_i > E₀ to the total ion energy density, measured in the downstream. We express this as

$$\begin{eqnarray}&&\varepsilon ({E}_{0})={\left\langle \displaystyle \frac{{U}_{i}({E}_{i}\gt {E}_{0})}{{U}_{i}({E}_{i}\gt 0)}\right\rangle }_{\mathrm{downstream}}\end{eqnarray} \tag{ 3 }$$

where U_i (E_i > E₀) is the ion energy density above the energy E₀ and the brackets indicate time averaging over a downstream interval. In this definition of ion acceleration efficiency, we assume that most available energy in the shock goes to the ion population, either by heating or acceleration. The benefit of this definition is that it naturally adopts a value between 0 and 1, and is less sensitive to uncertainties in upstream plasma measurements than definitions based on energy flux as used by e.g., Ellison et al. (1990). This makes the definition of ε based on energy density well suited for our study with many shock crossings and somewhat uncertain upstream conditions. The ion energy density can be expressed as

$$\begin{eqnarray}&&{U}_{i}({E}_{i}\gt {E}_{0})=4\pi \sqrt{\displaystyle \frac{2}{{m}_{i}^{3}}}{\int }_{{E}_{0}}^{{E}_{\max }}{\rm{d}}{E}_{i}\sqrt{{E}_{i}^{3}}{f}_{i}({E}_{i})\end{eqnarray} \tag{ 4 }$$

where f_i (E_i ) is the spherical mean of the ion phase-space density as a function of energy. Like Caprioli & Spitkovsky (2014a) we investigate ion acceleration efficiency above an energy E₀ = 10E_sw, where ${E}_{{\rm{sw}}}={m}_{i}{V}_{{\rm{sw}}}^{2}/2$ is the solar-wind bulk energy and V_sw the solar-wind bulk speed. We refer to these ions as energetic. We will however also be investigating ion acceleration above 5E_sw, which we refer to as suprathermal (Caprioli et al. 2015). To be able to calculate the ion acceleration efficiency we need to know at what parts of the time intervals MMS is in the downstream. We manually select the downstream times for all events. We set the downstream intervals from where the ion distribution is significantly thermalized and where the plasma density is increased. The time averaging in ε is then done in these intervals.

To truly represent acceleration efficiency, U_i should be measured in the downstream frame of reference. Therefore, we calculate f_i (E_i ) in the local downstream plasma frame. This is done by calculating the average phase-space densities in spherical shells centered on the plasma flow velocity. This allows for a direct comparison of acceleration efficiency at shock crossings with different downstream flow speeds.

Next, we look closer at three shock crossings by MMS from the data set. The events are shown in Figure 2. The manually determined downstream intervals are marked in gray. The first event shown in Figures 2(a)–(b) shows a quasiperpendicular shock. As typical for quasiperpendicular shocks, the transition from upstream to downstream is sharp. As can be seen, there are practically no ions with energy E > 10E_sw and therefore the ion acceleration efficiency is close to zero. Figures 2(c)–(d) shows an encounter with the quasiparallel bow shock close to the subsolar point. Here, the shock transition is much patchier and more extended. During the time interval, MMS observes three Short Large Amplitude Magnetic Structures (SLAMS) (Schwartz & Burgess 1991) seen as sharp increases in B in the upstream. The last SLAMS appears to be merging with the shock proper and the plasma transitions to downstream at this point. As we can see in Figure 2(d), there is a large density of energetic ions at this time and the acceleration efficiency is relatively high with 6% of the ion energy density being contained in energetic ions. The last shock crossing in Figures 2(e)–(f) also shows a quasiparallel shock crossing. In this event, MMS crosses the bow shock quite far out on the flank of the dawn side of the bow shock. Compared to the previous event, the Mach number M_A is also lower. Here, ions are clearly accelerated at the shock, but the energy of the ions rarely goes above 10E_sw, see Figure 2(f). Therefore, the acceleration efficiency is rather low at 0.1%. As we will see below, these events illustrate trends of ion acceleration in the statistics of all analyzed events.

Zoom In Zoom Out Reset image size

The Earth's bow shock is relatively small in astrophysical settings and has a curved shape. This means that, unless the IMF is perfectly aligned with the solar wind, a field line will only be connected to the bow shock for a limited time. Therefore, we can say that the field lines connected to the bow shock have a limited lifetime, even during stable upstream conditions. Ions that are injected and accelerated at the bow shock will eventually be convected away by the solar wind and therefore only have a limited amount of time to gain energy. We therefore try to estimate how long a field line has been connected to the bow shock when the MMS measurements are made and we call this the connection time T_c . Scholer et al. (1980) found that energetic ion events upstream of the shock are more probable for high T_c . Under most circumstances, a field line first connects to the bow shock at the tangent point where θ_Bn = 90°. The field line then traverses the bow shock with decreasing θ_Bn. To calculate T_c , we first use the same bow-shock model (Farris et al. 1991) we used to determine $\hat{n}$ to estimate the global shape of the bow shock. We then calculate the position of the tangent plane where the field lines connect to the bow shock where θ_Bn = 90°. We then get T_c from the distance between the spacecraft and the tangent plane in the direction of the solar wind. We illustrate this procedure in Figure 3, which depicts the bow shock seen from the Sun in the event shown in Figure 2(c)–(d). This method assumes a steady solar wind and we therefore only consider T_c < 600 s as this was the time we require the OMNI data to be steady, described above. The maximum achievable energy in DSA is expected to increase with T_c (Axford 1981; Caprioli & Spitkovsky 2014b). Here, we have no information about the maximum energy due to limited energy coverage (<28 keV). We instead use T_c to investigate the evolution of ion injection into DSA.

Zoom In Zoom Out Reset image size

Now, we look at all 154 selected shock crossings by MMS. Figure 4(a) shows the average phase-space density as a function of energy normalized to E_sw for all shock crossings averaged over the downstream time intervals and represented in the local downstream frame. We can see that downstream of the quasiperpendicular bow shock (shown in red colors), the ion distribution is heated and is nearly Maxwellian. There is also a clear bump in the distribution at 2–4E_sw. This is the result of ions being reflected at the shock and completing one gyration before going downstream, gaining energy in the process (Paschmann et al. 1980). Downstream of the quasiparallel shock (shown in blue colors), the ions are also heated but retain more of the solar-wind core at E_sw. A high-energy tail extending beyond 10E_sw is also present downstream of the quasiparallel shock. Next, we use these distributions to quantify the acceleration efficiency at the bow shock.

Zoom In Zoom Out Reset image size

Figure 4(b) shows a scatter plot of ε(10E_sw ), calculated from the distribution functions in Figure 4(a) as a function of θ_Bn. A clear trend is that the acceleration efficiency is close to zero for θ_Bn ≳ 50° in good agreement with (Caprioli & Spitkovsky 2014a). For shocks with θ_Bn≲50° the acceleration efficiency is significantly higher, up to ∼15 %. However, there is a large spread in ε for quasiparallel shocks. The color in Figure 4(b) indicates the Alfvén Mach number. There seems to be a trend of lower ε when M_A is low. This will be investigated further below. Another explanation for the large variations might be the dynamic and turbulent nature of quasiparallel shock waves and the fact that MMS measurements represent a snapshot at one point of the shock. There is also a slight trend for lower ε for θ_Bn < 20°. However, due to the $\sim \sin {\theta }_{{Bn}}$ probability for θ_Bn , there are few events in this range, making this result uncertain. Overall the acceleration efficiency agrees with simulation results by Caprioli & Spitkovsky (2014a).

We can clearly see in Figure 4 that ion acceleration efficiency has a strong dependence on the shock angle. Now, we further investigate how ε depends on other shock parameters. The clearest trends are found in the Mach numbers in Figures 5(a)–(b). It is clear that ε is significantly lower for lower Mach numbers M_A < 7 and M_ms < 5. Especially for M_ms where θ_Bn < 45°, the positive trend to ε seems to continue for the entire parameter range captured here. One clear trend is that ion acceleration is on average much more efficient close to the subsolar point than on the flanks of the bow shock, see Figure 5(c). This result is somewhat unexpected since, although the Mach numbers depend on ϕ, this dependence is too weak to explain the difference in ion acceleration between subsolar and flank bow shock. To further examine this, we study the ε dependence on connection time T_c normalized to the inverse ion gyrofrequency ω_ci in Figure 5(d). We see that T_c is strongly correlated with θ_Bn. In the quasiparallel regime, however, we see that an increasing connection time first leads to a somewhat higher acceleration efficiency, but that after $\sim 250{\omega }_{{ci}}^{-1}$ there is a decrease in ε. Perhaps surprisingly, we find no correlation between T_c and ϕ. This is because field lines do not have to pass by the nose of the bow shock but can be connected to flanks of the bow shock all the time.

Zoom In Zoom Out Reset image size

In order to test and validate the results from the MMS statistics, we compare the observations to a two-dimensional, hybrid-kinetic Vlasiator simulation of Earth's bow shock. Vlasiator (von Alfthan et al. 2014; Palmroth et al. 2018) is a global hybrid-Vlasov model where ions are described as distribution functions and electrons as a cold massless charge-neutralizing fluid. The simulation solves the Vlasov equation coupled with Maxwell's equations to propagate the plasma. The electric field is closed by Ohm's law E = − V _i × B + j × B /q_e n, where j is the current density, n is the number density, and q_e is the elementary charge. The main benefit to a Vlasov scheme is that it uses noise-free distribution functions directly instead of relying on integrating potentially noisy particles like in particle-in-cell simulations.

The simulation used in this study simulates near-Earth space in the meridional plane, ignoring dipole tilt. The simulation plane, which corresponds to the X–Z plane in GSE coordinates, is different from the MMS orbit, which is mainly in the X–Y plane. However, in the context of the bow shock, this simulation setup corresponds to a Parker spiral like setup (Turc et al. 2020). The solar wind in the simulation consists of a purely proton plasma and has a flow speed of 750 km s⁻¹ in the −X direction with a plasma density of 1 cm⁻³. The IMF is in the simulation plane at $5\,\mathrm{nT}\times (\hat{x}-\hat{z})/\sqrt{2}$ and is therefore southward and forms a 45° angle to the Sun–Earth line. This means that we model the quasiperpendicular and quasiparallel bow shock at the same time. The solar-wind Mach numbers are therefore M_A = 6.9, M_ms = 5.9, and ion beta β_i = 0.7. We note that the solar wind is faster than the typical solar wind, but the dimensionless parameters like Mach numbers and ion beta are typical for the solar wind. The extent of the simulation box is X ∈ [ − 48, 64]R_E and Z ∈ [ − 60, 40]R_E with a spatial resolution of 300 km. The 3D velocity space extends ±4000 km s⁻¹ in each direction with a resolution of 30 km s⁻¹. Vlasiator uses a sparse velocity grid in order to save on computational resources. This means that the simulation is only storing the distribution function where f exceeds some value ${f}_{\min }$ (Palmroth et al. 2018). In this case, ${f}_{\min }={10}^{-15}$ s³ m⁻⁶.

Figure 6(a) shows the plasma density in the simulation after a simulation time of 1361 s ($\sim 650{\omega }_{{ci}}^{-1}$). We can see the bow shock as a sharp increase in density. The quasiparallel part of the shock at Z < 0 is characterized by large fluctuations in density and magnetic field. The approximate position of the magnetopause is marked and we see a number of flux transfer events formed. One of the flux transfer events has launched a "bow wave" (Pfau-Kempf et al. 2016) toward the bow shock, which slightly distorts the shape of the bow shock. Normally in the simulation, the full 3D ion velocity distribution function is only saved in a few cells each time step. However, here we use data from a time that can be used to restart the simulation. For this time step, the ion distribution function is saved in all simulation cells, allowing for detailed analysis of ion acceleration at the bow shock.

Zoom In Zoom Out Reset image size

The Vlasiator simulation used here trades high spatial resolution for correctly modeling the global processes at the bow shock. The grid cell size in this simulation is 300 km while the upstream ion inertial length is ∼230 km. The thermal gyroradius of a solar-wind ion is ∼190 km and the gyroradius of a specularly reflected ion is on the order of 1000 km. This suggests that not all ion-kinetic physics can be expected to be fully resolved at the bow shock. Despite this, Pfau-Kempf et al. (2018) showed that even when not resolving the ion inertial length, the Vlasiator model is able to correctly model many kinetic effects such as ion reflection and upstream waves in an oblique 1D shock. In addition to this, the main strength of Vlasiator in this study is that it is a noise-free kinetic model while maintaining realistic separation of kinetic and global scales, which allows for direct comparison with observations from the Earth's bow shock.

In order to quantify ion acceleration in the simulation we need to estimate where in the simulation the bow shock is located. While assuming a bow-shock shape of a conic section (Schwartz 1998) is often the most accurate when working with spacecraft data, the bow shock in the simulation is distorted due to feedback from kinetic effects and its shape is not well described as a conic section or polynomial. The bow-shock model we are using here for estimating its position is constructed from contour points where the plasma density has a value of 2n_sw. This value is used as an approximation to when the plasma transitions from upstream to downstream since the compression ratio at all points of the shock is expected to be greater than two (Battarbee et al. 2020). The model bow-shock shape is then constructed by fitting a smooth curve to these contour points. The estimated bow-shock shape is shown as a magenta line in Figure 6(a) and the contour points used to make the fit in Figure 6(b). Important parameters such as M_A and θ_Bn are then calculated everywhere using the normal vector obtained from the model bow-shock curve. We emphasize that this postprocessing does not affect the physics of the simulation.

Figure 6(b) shows the energy density for ions with energy above 5E_sw (suprathermal ions) in a region close to and upstream of the estimated bow-shock position. We show the energy density above 5E_sw here instead of the 10E_sw that we used for the spacecraft observations due to the sparse velocity grid of Vlasiator, which we shall see below leads to the distribution function rolling off at around 10E_sw. This is a limitation in the Vlasov scheme of Vlasiator when studying high-energy particles. To get around this limitation we instead quantify the early stage of ion acceleration using 5E_sw as a lower energy limit and compare this to the observations with a similar energy limit.

We can see several interesting features in the concentration of suprathermal ions in Figure 6(b). First, it is clear that the quasiperpendicular bow shock produces practically no suprathermal ions and these ions are only found in the quasiparallel region. We can see the presence of a field-aligned beam upstream of the bow shock originating from the shock where θ_Bn ∼ 45° (Kempf et al. 2015). There seems to be no increase of suprathermal ions just downstream of where the field-aligned beam connects to the bow shock. Further out on the flank on the quasiparallel side the suprathermal ions are found close to the bow shock and extend both in the upstream and downstream. The extent of the suprathermal ions on the upstream side seems to increase further out on the flank. As the position on the flank translates to the time a field line has been connected, this indicates a DSA-like acceleration of the ions (Jones & Ellison 1991). The energy density of these ions is also highly structured in the upstream. The details of how the upstream suprathermal ion population evolves and how it is modulated by fluctuations like SLAMS-like structures in Vlasiator are left for a future study.

Next, we investigate the downstream ion distributions to quantify ion acceleration in the simulation. Therefore, we divide all simulation points in bins based on their angle to the Sun–Earth line. We use 256 equally-spaced bins in ϕ between −100° and 80^◦, where negative angles correspond to Z < 0. To better compare with the spacecraft observations which were made in direct contact with the shock, we only use simulation cells that are within a distance of 0.25R_E ≈ 1600 km to our bow-shock estimate. Lastly, we put a constraint for n>2n_sw in order to only use cells that we estimate to be downstream of the shock. The resulting average ε measured in the local downstream plasma frame is shown in Figure 6(c) along the Z position of the bow shock. It is clear from the figure that ε is higher on the quasiparallel side that there is a large spread in ε and that ε seems to decrease on the flank of the bow shock, all consistent with the observations.

The average distributions for each ϕ bin are shown in Figure 7(a). These downstream distributions show similarities to those observed by MMS in Figure 4(a) in that the high-energy tail extending further for quasiparallel shock geometries. However, as can be seen in Figure 7(a), the ion distributions appear to roll off at ∼10E_sw. This is due to the sparse velocity grid used in Vlasiator.

Zoom In Zoom Out Reset image size

One difference between the ion distributions observed by MMS and in Vlasiator is found in quasiperpendicular shock geometries. The simulation results do show a bump in the ion distributions at 2–4E_sw due to ion reflection when θ_Bn > 45°. It is, however, less clear than in the MMS observations. One explanation could be that the Mach number in the simulation is lower than most events in the MMS data set and that ion reflection is less efficient at lower Mach numbers. Another reason for the difference could be that Vlasiator does not correctly model the ion reflection at the quasiperpendicular shock, possibly due to the limited spatial resolution and that the electron pressure gradient term is omitted in Ohm's law. However, we estimate the electron pressure gradient term would contribute at most 10% of the shock-normal electric field at the quasiperpendicular shock. We also observe clear evidence of near-specular ion reflection at this part of the shock and that ∼10% of the solar-wind ions are reflected (c.f. Leroy et al. 1982). We shall further investigate ion distributions at the quasiperpendicular shock in a future study, as well as effects of spatial resolution.

The ion acceleration efficiency above 5E_sw as a function of θ_Bn is shown in Figure 7(b). The simulation results are shown on top of the MMS observations. To better compare the measurements to the simulation, we show the acceleration efficiency measured by MMS using only the part of the distribution with 5E_sw < E < 10E_sw. Like in Figure 4(b), we see that ion acceleration for this energy limit is much higher on the quasiparallel parts of the bow shock. The Vlasiator results show a qualitative match with the low Mach number shock crossings by MMS. Vlasiator also shows a large spread in acceleration efficiency in the quasiparallel part, further highlighting the dynamic nature of ion acceleration in this region of the bow shock. The simulation seems to show that ion acceleration starts to become efficient below θ_Bn ≲ 30° rather than θ_Bn ≲ 50° for the MMS observations. We also find no decrease in the acceleration efficiency for θ_Bn < 20°, suggesting that the simulation does not fully capture the θ_Bn-dependence on the ion acceleration efficiency. We note that the observed θ_Bn-dependence is more similar to that found in simulations by Caprioli & Spitkovsky (2014a) than what we find in the Vlasiator simulation. Except for this, ion acceleration efficiency in the Vlasiator simulation shows good quantitative agreement with observations of the bow shock by MMS.

Even though we are using only one time step of the simulation, the general curved shape of the bow shock allows us to investigate the time evolution of the suprathermal ions at the shock. Because of the two-dimensional simulation setup, determining the time a field line has been connected to the bow shock is somewhat simplified from a true 3D case. To determine the connection time T_c , we find the field line, which is connected to the tangent point where θ_Bn = 90° and determine T_c from the solar-wind speed and the distance along X between this line and the model bow shock. We estimate that the connection time is 120 s (250ω_ci ) at the subsolar point and 370 s (770ω_ci ) on the outermost flank in Figure 6. The bow-shock shape in the simulation stabilizes at a time before 800 s while the time step here is 1361 s, so the travel time estimates are not affected by simulation initialization and the large-scale bow shock shape changes.

It is clear from Figure 6(c) that the ion acceleration at the shock is less efficient on the far quasiparallel flank of the bow shock, similar to the MMS observations. In a 2D setup like this, the connection time T_c increases with Sun–Earth angle ϕ on the −Z flank. This would thus lead to the interpretation that, at very large values of T_c , the acceleration efficiency would decrease. Terasawa (1979) showed that in a curved shock, particle injection changes over connection time T_c as a function of local acceleration efficiency (in their model, magnetic compression and effective shock velocity). Like them, we also see a decrease in magnetic compression ratio at the quasiparallel flank, explaining how large values of T_c can in fact be connected with low-acceleration efficiency , and how T_c in itself is insufficient to estimate particle energization. Similarly, in the MMS observations, we see acceleration efficiency initially increase with T_c , then decreasing after $\sim 250{\omega }_{{ci}}^{-1}$ (see Figure 5(d)). This decrease happens despite the connection time having no correlation with the angle ϕ due to the 3D nature of the bow shock. These results, being in agreement, suggest that the relevant acceleration timescale to energies 5–10E_sw is much less than the field-line connection timescale.