Helium in the Earth's foreshock: a global Vlasiator survey

The foreshock is a region of space upstream of the Earth's bow shock extending along the interplanetary magnetic field (IMF). It is permeated by shock-reflected ions and electrons, low-frequency waves, and various plasma transients. We investigate the extent of the He²⁺ foreshock using Vlasiator, a global hybrid-Vlasov simulation. We perform the first numerical global survey of the helium foreshock and interpret some historical foreshock observations in a global context.

The foreshock edge is populated by both proton and helium field-aligned beams, with the proton foreshock extending slightly further into the solar wind than the helium foreshock and both extending well beyond the ultra-low frequency (ULF) wave foreshock. We compare our simulation results with Magnetosphere Multiscale (MMS) Hot Plasma Composition Analyzer (HPCA) measurements, showing how the gradient of suprathermal ion densities at the foreshock crossing can vary between events. Our analysis suggests that the IMF cone angle and the associated shock obliquity gradient can play a role in explaining this differing behaviour.

We also investigate wave–ion interactions with wavelet analysis and show that the dynamics and heating of He²⁺ must result from proton-driven ULF waves. Enhancements in ion agyrotropy are found in relation to, for example, the ion foreshock boundary, the ULF foreshock boundary, and specular reflection of ions at the bow shock. We show that specular reflection can describe many of the foreshock ion velocity distribution function (VDF) enhancements. Wave–wave interactions deep in the foreshock cause de-coherence of wavefronts, allowing He²⁺ to be scattered less than protons.

1 Introduction

The Earth's bow shock forms due to the interaction of the supermagnetosonic solar wind with our planet's magnetic field. As in other heliospheric shocks, solar wind particles interacting with the shock undergo a variety of processes, including reflection and acceleration. Upstream of the bow shock, in regions where plasma is magnetically connected to the shock, the reflected particles form a region called the foreshock. It is a very complex environment, populated by a variety of suprathermal ion distributions (Thomsen, 1985; Fuselier, 1995; Wilson, 2016), waves (Hoppe et al., 1981; Blanco-Cano et al., 2009; Wilson, 2016), and nonlinear transient structures (Kajdič et al., 2017; Blanco-Cano et al., 2018). The edges of the foreshock are magnetically connected to quasi-perpendicular regions of the Earth's bow shock (where the angle between the shock normal and the magnetic field θ_Bn ≳ 45^∘), whereas the central region of the foreshock is magnetically connected to the quasi-parallel bow shock (where θ_Bn ≲ 45^∘).

Most studies of the foreshock have concentrated on studying the proton dynamics and properties of ultra-low frequency (ULF) waves. Suprathermal ion distributions in the foreshock include field-aligned ion beams (FABs), gyrating distributions and hot diffuse populations. The original classification, based on the 2D International Sun–Earth Explorer 1 (ISEE) velocity distributions, also included intermediate ions (Thomsen, 1985). Subsequent observations with higher time resolution have, however, showed that intermediate distributions often display signatures of gyrating ions, which can be either isotropic or gyrophase bunched (Fuselier et al., 1986; Meziane et al., 2001). The interaction of suprathermal ions with the solar wind results in instabilities able to generate ULF waves (Gary, 1991).

Little attention has been given to the helium component, which is the most important minor species in the solar wind. Although helium constitutes typically only about 4 %–5 % of the total ion number density (Ipavich et al., 1984; Wurz, 2005; Gedalin, 2017), its contribution to the upstream mass density and dynamic pressure can be as large as 20 % (Gedalin, 2017). Thus, He²⁺ effects in shock dynamics and foreshock physics should not be ignored. Scholer et al. (1981) reported on ISEE observations of proton and alpha particle 30–36 keV∕q beams at the edge of the foreshock that exhibited similar time profiles. Ipavich et al. (1988) studied the content of ∼10 keV nuc⁻¹ H⁺ and He²⁺ in field-aligned beams with the AMPTE CCE spacecraft. They found that the alpha particles in the beams have approximately the same velocity as the H⁺ ions, but that the He²⁺ to H⁺ density ratio is dramatically smaller (2 orders of magnitude) than that measured simultaneously in the solar wind. Based on a study of 14 field-aligned beam events recorded with the ISEE satellite, Fuselier and Thomsen (1992) concluded that the ratio is roughly one-tenth of the solar wind ratio.

Fuselier et al. (1990) showed that two types of suprathermal He²⁺ distributions can be observed upstream of the quasi-parallel shock, namely a diffuse distribution (energetic; from several keV∕e up to the detector maximum) and a non-gyrotropic gyrating distribution. These gyrating He²⁺ distributions are observed near the shock, and their velocity components are consistent with near-specular reflection of a portion of the incident solar wind He²⁺ ions. The helium content in these gyrating populations can be roughly the same as in the pristine solar wind when the Alfvénic Mach number is M_A>7. These authors suggested at that time that the near-specularly reflected He²⁺ ions may be the seed population for the more energetic diffuse helium populations.

Using ISEE data, Fuselier et al. (1995) and Fuselier (1995) studied in more detail the origin of diffuse suprathermal ions with energies from a few up to $\sim \mathrm{100}\mathrm{keV}/e$. They found that the ratio of He²⁺ to H densities with suprathermal energies (normalized to solar wind abundances) is dependent on the location within the foreshock. High-energy field-aligned beams ($>\mathrm{10}\mathrm{keV}/e$) near the foreshock edges show significant ${\mathrm{He}}^{\mathrm{2}+}/\mathrm{H}$ ratios (near solar wind quantities), whereas lower energy beams ($\sim \mathrm{1}\mathrm{keV}/e$) deeper within the foreshock exhibit intermediate proton distributions and lower ${\mathrm{He}}^{\mathrm{2}+}/\mathrm{H}$ ratios. This difference in helium fraction was then assumed to be indicative of their origin. Low-energy beam production was explained in terms of magnetospheric leakage or shock reflection, whereas high-energy beams were attributed to shock drift acceleration, which is efficient for both protons and He²⁺. Additionally, Fuselier et al. (1995) found that, still deeper within the quasi-parallel foreshock, He²⁺ distributions are non-gyrotropic partial rings, whereas H distributions are ring beams and density ratios return to solar wind levels. Both of these distribution types are consistent with specular reflection of a portion of the incident solar wind (Fuselier et al., 1990).

Diffuse ion distributions are found throughout the deep foreshock, far from the foreshock edge. For them, the ratio of suprathermal He²⁺ to suprathermal H ion densities is similar to that of the solar wind composition, usually $n_{\mathrm{p}}/n_{\mathit{\alpha }}\sim \mathrm{4}\mathit{\%}$ Thus, Fuselier et al. (1995) suggested that the lower energy field-aligned beams with almost no helium content cannot be the seed population for diffuse ions. They proposed that the very energetic field-aligned beams at the edge of the foreshock propagate upstream much faster than the solar wind flow, and are confined to the edge, and therefore cannot contribute to the diffuse population observed further downstream. These results changed the original paradigm in which the origin of diffuse ions was explained in terms of field-aligned beams evolving into intermediate ions and then diffuse distributions (Thomsen, 1985). The fact that high-concentration He²⁺ gyrating ions are observed in the quasi-parallel foreshock, as are diffuse ions, suggests that gyrating distributions can be the seed population for the diffuse He²⁺ distributions. Similarly, gyrating H⁺ distributions are probably the source of the energetic diffuse H⁺ as well.

Using 1D hybrid simulations, Trattner and Scholer (1991) and Trattner and Scholer (1994) studied the acceleration of protons and He²⁺ ions at the quasi-parallel shock. They found that the concentration of helium in the diffuse population depends on the solar wind Mach number, plasma β, and the shock θ_Bn. In another numerical work, Trattner and Scholer (1993) investigated the thermalization of He²⁺ through the quasi-parallel shock and showed that, even if initially the heavier ions are less decelerated by the cross-shock potential, this difference disappears within a few gyroperiods downstream of the shock. The simulation results of Trattner and Scholer (1994) show a non-gyrotropic He²⁺ distribution re-entering the upstream region from the magnetosheath due to its large gyroradius, consistent with the shapes of the distributions observed by ISEE. These authors showed that He²⁺ ions can alter the shock structure and that the occurrence of He²⁺ ion clouds upstream of the shock is dependent on the Mach number.

In the recent years, the Hot Plasma Composition Analyzer (HPCA) instrument (Young et al., 2016) onboard the Magnetosphere Multiscale (MMS) spacecraft (Burch et al., 2016) has allowed new investigations of helium ions near the Earth's bow shock, providing, in particular, He²⁺ velocity distributions. Broll et al. (2018) investigated the reflection of He²⁺ at the quasi-perpendicular bow shock and showed that He²⁺ ions can undergo a similar specular-reflection process as the protons. The study of interstellar He⁺ pick-up ions at the quasi-perpendicular bow shock conducted by Starkey et al. (2019) revealed that single reflection at the shock plays a significant role in accelerating these ions.

In this work, we analyse He²⁺ properties in the foreshock using a global hybrid-Vlasov simulation of near-Earth space performed with the Vlasiator model (Palmroth et al., 2018). Vlasiator ion distribution functions have been compared to spacecraft observations of the foreshock in the past in Kempf et al. (2015). We investigate both the local and global properties of suprathermal He²⁺ ions and their possible influence on wave activity. We compare our results with MMS measurements in the Earth's foreshock and propose some new interpretations of some ISEE observations in the global context provided by our numerical simulation.

2 Methods and data

2.1 Simulations

We investigate the Earth's foreshock region using Vlasiator (Palmroth et al., 2018), a hybrid-Vlasov simulation capable of describing ion kinetics whilst encompassing global scales. Vlasiator solves the Vlasov equation for grid-discretized particle distribution functions, with closure provided by Ohm's law augmented by the Hall term. Electrons are considered a charge-neutralizing massless fluid with no electron pressure gradient term. We investigate the foreshock using a 2D–3V simulation, with 3D moments and velocity distribution functions but a 2D spatial domain. The geocentric solar ecliptic (GSE) simulation domain is $X\in (-\mathrm{48.66}{R}_{\mathrm{E}};\mathrm{64.35}{R}_{\mathrm{E}})$ and $Z\in (-\mathrm{59.65}{R}_{\mathrm{E}};\mathrm{39.24}{R}_{\mathrm{E}})$ and a single cell width in the Y direction. The spatial resolution is 300 km (1.3 times the solar wind ion inertial length) and the velocity–space resolution for both ion species (protons and alpha particles) is 30 km s⁻¹.

Our simulation is set to have solar wind values of β=0.7, M_ms=5.6, and M_A=6.9. We initialize the simulation with a solar wind of $n_{\mathrm{p}}=\mathrm{1}{\mathrm{cm}}^{-\mathrm{3}}$ and $n_{\mathit{\alpha }}={\mathrm{10}}^{-\mathrm{2}}{\mathrm{cm}}^{-\mathrm{3}}$. Due to the mass ratio, we set T_p=0.5 MK, and T_α=1.0 MK. The solar wind speed is set to $u_{\text{sw}}=\mathrm{750}{\mathrm{km}\mathrm{s}}^{-\mathrm{1}}$ in the $-{\widehat{e}}_x$ direction, simulating fast solar wind conditions and ensuring efficient simulation initialization. Despite the use of fast solar wind conditions, the Alfénic Mach number (7) and plasma beta (0.7) are typical for a variety of solar wind conditions, ensuring the validity of the initial conditions.

We set an interplanetary magnetic field (IMF) of B_x=3.54 nT and $B_z=-\mathrm{3.54}\mathrm{nT}$, resulting in a 45^∘ cone angle. Although the simulation plane is meridional, the foreshock dynamics are comparable with an ecliptical Parker spiral setup. The Earth's magnetic dipole is a ${\widehat{e}}_z$-aligned line dipole, resulting in a realistic magnetopause standoff distance (Daldorff et al., 2014). The simulation has an inner boundary at $\mathrm{3}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{km}\approx \mathrm{4.7}{R}_{\mathrm{E}}$, modelled as a perfectly conducting ionosphere. Thus, the run is nearly identical to the simulation presented in Blanco-Cano et al. (2018), with the addition of alpha particles as an independent, self-consistent species. Our alpha particle density is set to only 1 % of the solar wind, so mass loading and effects on ULF wave properties are expected to be small. In order to constrain memory usage, we set the minimum stored phase–space densities (as explained in von Alfthan et al., 2014) to $f_{\text{min,p}}={\mathrm{10}}^{-\mathrm{15}}{\mathrm{s}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{6}}$ and $f_{\text{min,}\mathit{\alpha }}={\mathrm{10}}^{-\mathrm{17}}{\mathrm{s}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{6}}$.

Additionally, in subsection 3.1 we compare our main simulation to an equatorial Vlasiator simulation with a spatial resolution of 228 km and solar wind values of β=2.3, M_ms=5.9, M_A=10, $n_{\mathrm{p}}=\mathrm{3.3}{\mathrm{cm}}^{-\mathrm{3}}$, $n_{\mathit{\alpha }}=\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{2}}{\mathrm{cm}}^{-\mathrm{3}}$, and $u_{\text{sw,}x}=-\mathrm{600}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$. The IMF is set to 5 nT, with a 5^∘ cone angle.

2.2 Observations

In this study, we also analyse observations from the MMS mission (Burch et al., 2016) in the Earth's foreshock. We use data from three different instruments, namely the fluxgate magnetometer (FGM; Russell et al., 2016), the Hot Plasma Composition Analyzer (HPCA; Young et al., 2016), and the dual ion spectrometers (DIS; Pollock et al., 2016), which is part of the fast plasma investigation (FPI) instrument. For all instruments, we use survey-mode data only. The FGM produces magnetic field measurements with 16 s⁻¹ time resolution in fast survey mode.

Energy–time spectrograms from HPCA measurements are used to identify whether the spacecraft are located in the foreshock or in the solar wind. We also calculate partial densities for the different ion species, following the procedure described in the HPCA science algorithms and user manual available on the MMS Science Data Center web page (https://lasp.colorado.edu/mms/sdc/public/datasets/hpca/, last access: 1 October 2020). HPCA data are made available in survey or burst modes. Data are collected in 0.5 spin (10 s) time resolution at 64 energies, 16 azimuth angles, and 16 elevation angles for five different species (H⁺, He⁺, He²⁺, O⁺, and O²⁺). In the survey-mode data used here, the data are reduced in energy and angles to 16 energies, eight azimuths, and eight elevation angles. The total energy range of the data is between $\sim \mathrm{1}\mathrm{eV}/e$ and 40 keV∕e, with a mass resolution of $M/\mathrm{\Delta }M\sim \mathrm{8}$.

Finally, we calculate the ion partial density (without species separation), using measurements from a FPI–DIS, as a means of comparison to the partial densities obtained from HPCA. The FPI produces burst sky maps which consist of ion count arrays of 32 energies × 32 azimuth angles × 16 polar angles that are accumulated every 150 ms. Then, 30 consecutive DIS burst sky maps are summed in order to produce the survey-mode sky maps with a time resolution of 4.5 s. The energy range of the DIS is between 10 eV∕e and 30 keV∕e.

3 Results

In Fig. 1, we show an overview of the foreshock region in the simulation. Figure 1a shows out-of-plane magnetic field B_y fluctuations at time t=1100 s, showcasing the ULF wave fronts seen throughout the foreshock, displayed on a symmetric logarithmic colour scale. The black curves indicate magnetic field lines. Figure 1b also displays the ULF foreshock extent as black contours drawn for $B_y=\pm \mathrm{0.1}\mathrm{nT}$ at time t=1100 s, with the diverging colour map indicating the relative abundances of foreshock suprathermal ions, scaled to the incoming solar wind number density ratio. The solar wind and/or foreshock thermal ion distribution was accumulated from the velocity space constrained within a sphere of 500 km s⁻¹, centred at the solar wind speed $u_{\text{sw},x}=-\mathrm{750}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$, and all ions outside this sphere were considered part of the suprathermal distribution. Suprathermal particle measurements were averaged over 4 min (between 1080 and 1220 s), providing an overview of a steady state foreshock, smoothing over effects due to the gyration of particle populations at the foreshock edge. An animated version of Fig. 1b, showing instantaneous density ratios instead of the 4 min average, is provided as Supplementary video A.

Figure 1 Overview of the foreshock region of the simulation. (a) Out-of-plane magnetic field fluctuations B_y at 1100 s, using a symmetric logarithmic colour scale, indicating the extent of the ULF foreshock region. Black curves indicate magnetic field lines. (b) Ratio of suprathermal densities for helium over proton, normalized to the solar wind ratio of 1 %, averaged over a period of 4 min. The ratio is not shown in the pristine solar wind. Black contours are drawn for $B_y=\pm \mathrm{0.1}\mathrm{nT}$ at time t=1100 s.

3.1 Foreshock edge

As shown in Fig. 1, the ion and ULF foreshocks are not identical in extent. The ULF foreshock edge visible in both panels connects to the bow shock at $Z\sim -\mathrm{5}{R}_{\mathrm{E}}$, whereas the ion foreshock edge intersects the bow shock already at $Z\sim +\mathrm{5}{R}_{\mathrm{E}}$ (Fig. 1b). At the bottom of the figure, at $Z=-\mathrm{50}{R}_{\mathrm{E}}$, we see the ion foreshock extending up to X=40 R_E, whereas the ULF foreshock only extends to a position ∼10 R_E further downstream. We also see that, in Fig. 1b, the ratio of suprathermal alphas to protons in the foreshock shows significant deviations from the incoming solar wind ratio of 1 ∕ 100. Throughout most of the deep foreshock, the $n_{\mathit{\alpha }\text{,st}}{n}_{\text{p,st}}^{-\mathrm{1}}\times {\mathrm{10}}^{\mathrm{2}}$ ratio, shown in Fig. 1b, tends to values ≳2, whereas at the foreshock edge, it falls below 0.2. This abundance of He²⁺ within the deep foreshock is likely a computational artefact, with H ions being efficiently scattered into the diffuse population, which is not tracked, whereas He²⁺ remains in non-gyrotropic partial rings and clumps. Right at the edge of the foreshock we see spatially periodic structures with, again, large alpha-to-proton ratios, likely caused by bursty reflection at the quasi-perpendicular shock and the alpha particles having larger gyroradii, thus gyrating further into the upstream. The ULF and ion dynamics of the foreshock edge are further examined in Fig. 2.

Figure 2 (a) Magnetic field out-of-plane component B_y indicating the ULF foreshock extent. Contours indicate proton (black) and helium (green) suprathermal densities at 0.1 % (dotted contours) and 0.5 % (solid contours) of solar wind values. Four thick lines indicate the cut-throughs across the foreshock boundary. (b–e) Profiles across the foreshock at the positions shown in (a). Profiles are shown for suprathermal proton density (black solid line), suprathermal helium density (green solid line), and magnetic field magnitude (blue solid line). The horizontal dashed lines indicate 0.5 % (dashed line) of the solar wind density.

Figure 2a shows an excerpt from Fig. 1a, featuring a portion of the foreshock edge. The colour map is again the out-of-plane magnetic field component, overlaid with contours of proton (black) and helium (green) suprathermal densities. We focus on the contours upstream of the ULF foreshock. The solid contours are drawn where the suprathermal density is 0.5 % of the species' solar wind density and the dotted contours at 0.1 %, respectively. The thick grey lines indicate four cut-throughs across the foreshock edge. The magnetic field strength and the suprathermal proton and helium densities along these cuts are plotted in Fig. 2b–e. As shown in these plots, there is variation in the profiles of ions across the foreshock edge, moving into the foreshock from right to left. Close to the bow shock (Fig. 2b–c), there is a somewhat rapid increase in ion densities near 10 R_E, followed by a gradual increase over the following 1–2 R_E and finally plateauing values. Further out (Fig. 2d and e), we see a more gradual increase in suprathermal ion densities over several R_E and a less clear plateau. The suprathermal ion density threshold at 0.5 % of the species solar wind density is shown as a horizontal dashed line.

We note that non-thermal particles are found several R_E upstream of the foreshock waves, consistent with previous works showing that no ULF wave activity is observed in conjunction with the region closest to the foreshock edge where field-aligned beams are expected (and found). More importantly, we find significant amounts of suprathermal helium throughout the field-aligned beam region and into the ULF foreshock. As shown by the black and green contours in Fig. 2a, the helium foreshock edge is located equal to or slightly downstream of the proton foreshock edge. This shift is more pronounced in the dotted contours and increases when moving further away from the bow shock. At about 30 R_E from the shock, the difference is of the order of 1 R_E. The shift of the foreshock edge is also noticeable in the line profiles (Fig. 2b–e). This suggests that the proton foreshock is slightly more extended than the helium foreshock, and measurements made at the very edge of the foreshock can suggest significantly lower helium fractions despite the helium abundances rising to proton-comparable levels a few R_E deeper within the foreshock. We also point out that the plateau of alpha particles inside the foreshock edge has more fluctuations to it than the plateau of protons.

We also note that the suprathermal ion density contours in Fig. 2 and the time-averaged ratios seen in Fig. 1 are not smooth, having a wavy shape instead. Supplementary Video A shows the wavy shape is due to intermittent reflection of ions at the bow shock, with bursts of ions propagating away from the shock along the field lines. This periodic and intermittent enhanced reflection of particles may be due to mesoscale reformation of the bow shock (Battarbee et al., 2020e) and can cause further discrepancies and variation in field-aligned beam densities.

Figure 3 MMS3 crossing from the solar wind into the foreshock region on 30 December 2018. (a) The MMS position relative to the bow shock just before the foreshock crossing in a plane that contains B. The solid blue line shows the bow shock position in the plane containing MMS and the dashed line in the plane containing Earth. The grey lines are the magnetic field lines. (b) The survey-mode measurements of the magnetic field magnitude and components, showing the spacecraft entering the ULF foreshock. (c) The number densities for H⁺ and He²⁺ suprathermal ions measured by HPCA with spacecraft frame energy cutoffs indicated in the legend, with the helium density multiplied by 10 to ease density evaluation. The FPI ion number density is shown for comparison. (d, e) The HPCA ion energy spectrograms for protons and helium, with the dashed line indicating the respective suprathermal ion cutoff energies.

Figure 4 MMS3 crossing from the solar wind into the foreshock on 18 November 2018. The format is the same as in Fig. 3 but with He²⁺ multiplied by 25.

We now compare our numerical results with observations from the MMS spacecraft at the foreshock edge. Figures 3 and 4 show two time intervals during which the MMS3 spacecraft crossed from the solar wind into the foreshock region. Figures 3a and 4a show the spacecraft position relative to the bow shock model, which is scaled with the solar wind dynamic pressure (Farris et al., 1991; Schwartz, 1998). The plane is chosen so that B is in the $X_{\text{GSE}}-Y^{\prime }$ plane and that Y^′ is in the positive Y_GSE direction. Figures 3b–e and 4b–e, from top to bottom, show the IMF magnitude and components, the suprathermal densities of H⁺ (black; HPCA), He²⁺ (green; HPCA), and all ions (red; FPI) and proton and He²⁺ energy–time spectrograms from the HPCA instrument. The He²⁺ suprathermal densities have been multiplied by 10 or 25 to ease comparison of suprathermal ion density gradients. The lowest energy used in the calculation of the suprathermal (partial) densities for each species is also stated on the panel. This energy is about 4 times higher for alphas compared to the protons, due to larger mass of He²⁺. We note that the energy–time spectrograms show energy per charge $E_i{Z}_i^{-\mathrm{1}}$, with the lowest energy included in the suprathermal population shown with the dashed line set to the values indicated in Figs. 3c and 4c. Also, only those suprathermal ions measured by HPCA (up to 40 keV∕e) are included in the densities, but we estimate the lost contribution due to higher energy ions to be small due to very small phase–space densities.

During both time intervals the transition between the foreshock and the undisturbed solar wind was not due to any IMF rotational discontinuity, as can be seen from the MMS magnetic field components. We also checked the solar wind measurements propagated to the bow shock nose from the OMNI database (King and Papitashvili, 2005) and found no sharp IMF change during those intervals. This means that the foreshock was undergoing gradual motion due to slow IMF rotations, so the spacecraft did not observe travelling foreshocks (Kajdič et al., 2017) or foreshock cavities (Schwartz et al., 2006; Billingham et al., 2008, 2011). This gradual motion of the foreshock region over the spacecraft location allows for magnetic field and suprathermal density profiles to be compared with those obtained from our simulations (Fig. 2).

Figures 3c and 4c show that the proton and He²⁺ suprathermal densities do not always behave in the same manner. In Fig. 3c, we can clearly see that the proton and He²⁺ suprathermal densities are well correlated across the foreshock edge, although H⁺ has a slightly stronger initial beam before 22:06 Coordinated Universal Time (UTC). Both increase with a similar relative gradient as the spacecraft crosses from the solar wind into the foreshock. In Fig. 4c, the suprathermal proton density starts to increase ∼3 min before the suprathermal He²⁺ density and remains low well into the foreshock. This suggests that the proton foreshock extends further outward than the He²⁺ foreshock. In other words, the proton foreshock extends to field lines connected to larger θ_Bn values than the He²⁺ foreshock.

In order to better understand the different foreshock edge crossing behaviour seen in Figs. 3 and 4, we compare the main Vlasiator simulation used in this study to an equatorial plane simulation with a quasi-radial IMF, as detailed in Sect. 2.1. In Fig. 5, we show an enlarged view of the foreshock edge of both runs, with colour maps indicating the out-of-plane component of the magnetic field and black and green contours indicating proton and helium suprathermal densities at 0.5 % (solid) and 0.1 % (dotted) inflow densities, averaged over 2 min. In comparing the Vlasiator plots with MMS observations, it is important to keep in mind that the Vlasiator definition for suprathermals includes reflected particles with simulation frame energies comparable with the solar wind bulk, whereas the MMS suprathermal density is defined with a markedly higher minimum energy threshold. Comparison of Fig. 5a–b indicates that simulations can reproduce the two different foreshock edge behaviours, with Fig. 5a (our main run) presenting a relatively rapid drop-off and Fig. 5b (the quasi-radial IMF comparison run) showing a much more gradual fall-off of suprathermal ion densities and a greater difference between the two ion species.

The IMF prior to the 30 December 2018 foreshock edge crossing had the IMF pointing in the dawn, anti-sunward direction with its cone angle $\sim \mathrm{55}{}^{\circ }$. Consequently, this was also the orientation of the foreshock. At the time, MMS3 was located near the nose of the Sun–Earth line, at approximately $(\mathrm{12.5},\mathrm{0.8},\mathrm{5.1})R_{\mathrm{E}}$ in GSE coordinates. This location, together with large IMF cone angle, means that the foreshock edge was crossed in a way that is qualitatively similar to that shown in Fig. 5a.

In the case of the 18 November 2018 event, the IMF cone angle was $\sim \mathrm{35}{}^{\circ }$ and pointing in the northern, anti-sunward direction prior to the foreshock encounter, and consequently, this was also the foreshock orientation. The GSE coordinates of MMS3 were $(\mathrm{9.8},\mathrm{11.3},\mathrm{5.7})R_{\mathrm{E}}$, i.e. far from the Sun–Earth line. Thus, this foreshock edge crossing is comparable with the one shown in Fig. 5b.

Figure 5 An enlarged view of the foreshock edge, with the magnetic field out-of-plane component B_oop indicating the ULF foreshock extent. Contours indicate suprathermal ion densities, as in Fig. 2a, but averaged over 2 min. (a) The meridional plane 45^∘ IMF simulation used in the majority of this paper. (b) An equatorial plane 5^∘ IMF run for comparison. The IMF orientation appears to affect the suprathermal ion density gradient at the foreshock edge.

3.2 Velocity distribution functions and their properties

We now examine the properties of ion velocity distribution functions (VDFs) at the edge of and within the foreshock. Figure 6 shows 2D projections of proton and helium velocity distribution functions in the solar wind frame at three positions close to the foreshock edge, extracted from our Vlasiator simulation at time 1000 s. In Fig. 6 (panels 1–2 and 4, subpanels $(\parallel ⟂)$), we show projections that have been generated by averaging the instantaneous VDFs over the v_B×V direction, whereas subpanels labelled $(⟂⟂)$ have been averaged over the v_B direction. We use two different colour scales to differentiate the phase–space densities of the two ion species, with ranges selected to account for the input solar wind abundance ratio. Figure 6 (panels 1–2 and 4) shows VDFs from virtual spacecraft locations (A, B, and C, respectively). The Fig. 6 $(\parallel ⟂)$ subpanels also show an ellipsoid located at the position in velocity space where particles would end up if they were specularly reflected from the solar wind population at the closest point of the bow shock. The shock location was determined according to a plasma compression criterion of n_p>2n_p,sw, and the shock normal direction was estimated to be equal to a vector pointing from the closest shock location to the virtual spacecraft location. We emphasize that this estimate is a rough one, to be improved upon in future studies, and does not account for ion propagation times, drifts, or the existence of an electron pressure gradient cross-shock potential. Figure 6 (panel 3) shows an overview of the foreshock region, indicating the locations of the spacecraft on top of a colour map depicting the temperature anisotropy calculated from the whole proton VDF. The dark orange region at the upstream edge of the foreshock, with parallel temperatures in excess of perpendicular temperatures, is indicative of the FAB region of the proton foreshock. Magnetic fields lines are depicted with black curves.

In panel 1 of Fig. 6, at the position of virtual spacecraft A located very close to the foreshock edge, we see very clear proton and helium FABs, though the helium beam appears to have more structure. Parallel velocities of both beams are between 1500 and 2000 km s⁻¹ in the solar wind frame or $\mathit{\gtrsim }{u}_{\text{sw}}=\mathrm{750}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$ in the simulation frame, which translates to roughly 3–5 keV nuc⁻¹. We note that this is larger than the energy at which specularly reflected particles would be found, which suggests that these particles have experienced shock drift acceleration (SDA) at the quasi-perpendicular shock front. This indicates that the source of the FABs in panel 3 of Fig. 6 is located at roughly $X=\mathrm{16}{R}_{\mathrm{E}};Z=\mathrm{0}{R}_{\mathrm{E}}$. Panel 2 of Fig. 6 depicts VDFs at position B, at the boundary between the FAB region and the ULF foreshock. At this location, the ion beam seems to be transitioning into a more gyrating distribution, with a gyrophase-bunched signature (visible in the $(⟂⟂)$ panel at $v_{B\times V}>\mathrm{500}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$) for both species but particularly for helium at $v_{B\times V}>\mathrm{1000}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$. Parallel velocities have begun to decrease, extending down to 1000 km s⁻¹ in the solar wind frame of reference. A similar upstream ion velocity decrease was also seen in Battarbee et al. (2020e). Still, at this position, both helium and protons show a qualitatively similar shape to their distribution functions. Conversely, in panel 4 of Fig. 6, at position C, the proton VDF resembles a low-energy field-aligned beam, extending from 1000 to 1800 km s⁻¹ in v_B, whereas the helium distribution has very little trace of a beam, instead consisting of a broken up gyrating ion population. An animated version of Fig. 6 can be found as Supplementary Video B. We also note that, early in this video, at virtual spacecraft A, the helium VDF in particular shows what appears to be a gyrophase-bunched population but which remains stationary in velocity space. We deduct that it is in fact spatial sampling of the helium foreshock edge, visible due to the large gyroradius of alpha particles. We also point out that the proton beam energy found, in particular, in panel 1 of Fig. 6 does not extend to the tens of keV∕e found in some spacecraft observations, possibly a result of the sparse velocity space implementation in Vlasiator.

Figure 6 Velocity distribution functions (VDFs) for protons and alpha particles and their locations in the foreshock. (3) Map of the foreshock, with three VDF positions indicated with capital letters. The colour indicates the proton temperature anisotropy, with magnetic field lines in black. Panels (1), (2), and (4) show sets of four projections of ion VDFs at virtual spacecraft locations in the solar wind frame. Coordinates are chosen to represent two positions at the foreshock boundary (A and B) and one just within the ULF foreshock (C). In each set, the subpanels are labelled as v_B versus $v_{B\times (B\times V)}$ ($\parallel ⟂$; left columns) or v_B×V versus $v_{B\times (B\times V)}$ ($⟂⟂$; right columns) and protons (top rows) or helium (bottom rows). The black and white circles estimate the areas of specular reflection.

Figure 7 is similar to Fig. 6 but with virtual spacecraft locations chosen to represent regions further within the foreshock and also closer to the quasi-parallel bow shock. Panels 1 and 2 of Fig. 7 show spacecraft D and E, which are located close to the bow shock and depict mostly gyrating and intermediate ion populations. The helium populations of gyrating ions appear to have more structure to them. It is also noteworthy that both proton and helium $(\parallel ⟂)$ subpanels of panel 1 in Fig. 7 have what looks like a FAB propagating towards the shock at a parallel velocity greater than the solar wind speed, with velocities v_B<0. Although the estimate of the velocity space location of specularly reflected particles does not account for particle travel times or shock shape evolution, we find that the near-circular ellipsoids indicating potential specular reflection are usually found in VDF regions where there are enhancements, suggesting that specular reflection indeed plays a role in the generation of these populations. As the bow shock reforms as a mesoscale process with distinctly non-planar features, the direction of specular reflection varies, leading to the intermittent and gyrating partial rings seen in these VDFs. This mapping of specular reflection can also be seen in Supplementary Video C, which is an animated version of Fig. 7. Panel 4 of Fig. 7 shows proton and helium VDFs further within the deep foreshock and away from the bow shock. The proton population appears to be a gyrating ion population, and the helium population is a low-energy beam population, but viewing the VDFs at different times (see Supplementary Video C) shows that both ions usually resemble gyrating ion populations.

Figure 7 Velocity distribution functions for protons and alpha particles and their locations in the foreshock. (3) Map of the foreshock with three VDF positions indicated with capital letters. The colour indicates proton temperature anisotropy, with magnetic field lines in black. Panels (1), (2), and (4) show sets of four projections of ion VDFs at virtual spacecraft locations in the solar wind frame. Coordinates are chosen to represent two positions close to the quasi-parallel bow shock (D and E) and one deep within the outer foreshock (F). In each set, subpanels are labelled as v_B versus $v_{B\times (B\times V)}$ ($\parallel ⟂$; left columns) or v_B×V versus $v_{B\times (B\times V)}$ ($⟂⟂$; right columns) and protons (top rows) or helium (bottom rows). The black and white circles estimate the areas of specular reflection.

In light of the complex simulated VDF shapes, and in order to obtain a better understanding of global foreshock ion characteristics, we show, in Fig. 8, plots of global per-species temperature anisotropies and a measure of per-species non-gyrotropy. We apply the temperature anisotropy to the whole VDF, allowing us to identify regions where VDFs have FAB-like features showing up as anisotropy values much smaller than 1. The agyrotropy measure (Swisdak, 2016) is as follows:

$$\begin{array}{}\text{(1)}& Q_{\text{ag}}={\displaystyle \frac{P_{\mathrm{12}}^{\mathrm{2}}+P_{\mathrm{13}}^{\mathrm{2}}+P_{\mathrm{23}}^{\mathrm{2}}}{P_{⟂}^{\mathrm{2}}+\mathrm{2}{P}_{⟂}{P}_{\parallel }}},\end{array}$$

where $P_{\parallel }=P_{\mathrm{11}}$ is the parallel pressure, $P_{⟂}=\mathrm{0.5}(P_{\mathrm{22}}+P_{\mathrm{33})}$ is the perpendicular pressure, and P₁₂, P₁₃, and P₂₃ are off-diagonal pressure tensor components that can be used to evaluate the complexity of the VDF and the role of gyrating ions. For a completely gyrotropic distribution, Q_ag=0 and Q_ag=1 would signify a maximal deviation from gyrotropy. Figure 8a–b show agyrotropies for protons and helium, respectively, and Fig. 8c–d show temperature anisotropies. The panels in Fig. 8 also show out-of-plane magnetic field contours at a level of ±0.2 nT, indicating the extent and structure of the ULF foreshock. We find that the FAB region of protons is clearly visible in Fig. 8c at the edge of the ion foreshock as a dark orange band ($T_{⟂,\mathrm{p}}{T}_{\parallel ,\mathrm{p}}^{-\mathrm{1}}\approx \mathrm{0.2}$). In Fig. 8a we see the same structure at the foreshock edge. It is visible as a band of medium green blobs (Q_ag,p∼0.01) close to the shock nose, transitioning to paler green thin streaks (Q_ag,p∼0.001) away from the shock.

Figure 8 Foreshock characteristics for proton and helium at 1100 s. (a, c) Protons and (b, d) helium are shown. (a, b) Agyrotropy Q_ag (see Eq. 1), where a value of 0 corresponds with perfect gyrotropy. (c, d) Temperature anisotropy $T_{⟂}{T}_{\parallel }^{-\mathrm{1}}$ calculated for the total distribution function, including both solar wind and suprathermal parts. Black contours show the out-of-plane magnetic field B_y fluctuations at a level of ±0.2 nT, indicating the extent of the ULF foreshock region.

For helium, in Fig. 8d, this FAB region at the edge of the foreshock is also clearly visible, with very low anisotropy values. Interestingly, helium also shows a parallel pressure signature (dark orange low anisotropy values) deeper in the foreshock. This region coincides with a weakening of the ULF foreshock, visible in the destructuring of wave fronts (black contours) in the vicinity of ($X=\mathrm{10}{R}_{\mathrm{E}};Z=-\mathrm{40}{R}_{\mathrm{E}}$). Referring to Fig. 1b, we see that this band of low-energy, FAB-type alpha particles coincides with an increase in the measured helium fraction, which continues downstream from that point.

Figure 8a–b show, in particular for protons, an enhancement in agyrotropy at the boundary between the FAB region and the ULF foreshock. This dark green region is a signature of gyrating and gyrophase-bunched ions at this boundary, in agreement with Meziane et al. (2004), Mazelle et al. (2007), and Andrés et al. (2015). This effect can also be seen at virtual spacecraft location B in Fig. 6 (panel 2; ($p⟂⟂$) and ($\mathit{\alpha }⟂⟂$) subpanels), as a gyrophase-bunched extension of the ion distribution. We note that this increase in agyrotropy and, thus, gyrating ions matches the ULF foreshock and previous studies well, down to about $Z=-\mathrm{25}{R}_{\mathrm{E}}$, but the boundary becomes less well defined further out, away from the shock. We also note that Fig. 8b shows darkened bands of increased agyrotropy on both the upstream and downstream edges of the inner-foreshock-heightened parallel pressure alpha particle band, which suggests there might be similar gyrophase sampling taking place as for the foreshock FAB beam proper.

For both protons and helium, we see striped enhancements of agyrotropy right at the outer ion foreshock boundary, indicative of spatial sampling of gyrating ion beams right at the outermost foreshock edge. The gyration of these ion beams, accelerated at the quasi-perpendicular shock front and made non-uniform by the rippling of the shock front, is particularly visible in Supplementary Video D, which is an animated version of Fig. 8. Finally, we note that large regions of the foreshock close to the quasi-parallel bow shock show signatures of temperature anisotropies ≳1 and enhanced agyrotropies, which are likely a signature of specular reflection of ions at the quasi-parallel bow shock.

3.3 Foreshock waves

Figure 9 Time series and wavelets at virtual spacecraft positions C (a–h) and D (i–p) at the edge of the ULF foreshock. Panels (a–d) and (i–l) show proton (helium) quantities in black (green). (a, i) Total number densities. (b, j) Suprathermal number densities on a logarithmic scale. (a, b, i, and j) The helium number densities have been multiplied by 100. (c, k) Temperature anisotropies. (d, l) Agyrotropy Q_ag measured as a 5 s running average. (e, m) Magnetic field magnitudes $\left|B\right|$. (f, n) Magnetic field out-of-plane components B_y. (g, o) Wavelet power spectra of $\left|B\right|$ fluctuations. (h, p) Wavelet power spectra of B_y fluctuations. The black contours show the 95 % confidence level, and the cross-hatched regions denote the cone of influence. The white dashed lines show the expected frequency of B_y fluctuations (see text).

Figure 10 Time series and wavelets at virtual spacecraft positions E (a–h) and F (i–p) in the ULF foreshock. Panels (a–d) and (i–l) show proton (helium) quantities in black (green). (a, i) Total number densities. (b, j) Suprathermal number densities on a logarithmic scale. (a, b, i, and j) The helium number densities have been multiplied by 100. (c, k) Temperature anisotropies. (d, l) Agyrotropy Q_ag measured as a 5 s running average. (e, m) Magnetic field magnitudes $\left|B\right|$. (f, n) Magnetic field out-of-plane components B_y. (g, o) Wavelet power spectra of $\left|B\right|$ fluctuations. (h, p) Wavelet power spectra of B_y fluctuations. The black contours show the 95 % confidence level, and the cross-hatched regions denote the cone of influence. The white dashed lines show the expected frequency of B_y fluctuations (see text).

Figures 9 and 10 show measurements from virtual spacecraft placed at locations C, D, E, and F in the foreshock (see Figs. 6 and 7). The top row in each plot displays the total number density n_tot for protons and helium, and the second row displays their suprathermal number densities n_st. In both rows, the helium number densities have been multiplied by a factor of 100 to ease the comparison with proton number densities. Note that the suprathermal number densities are shown on a logarithmic scale, since the values vary significantly. In the third row, temperature anisotropy $T_{\perp }{T}_{\parallel }^{-\mathrm{1}}$ is displayed, and the fourth row displays the agyrotropy measure Q_ag. The time series of Q_ag have been smoothed out with a 5 s running average to help readability due to a large number of high-frequency fluctuations. The fifth and sixth rows display the total magnetic field $\left|B\right|$ and its out-of-plane component B_y, respectively.

The fluctuations of H⁺ and He²⁺ total densities follow each other quite closely at all locations, though the amplitude of the He²⁺ density oscillations is larger than that of protons at point D and E. These larger He²⁺ density variations seem to be well correlated with the fluctuations of the He²⁺ suprathermal density at position D (Fig. 9j) but not so much at point E (Fig. 10b). At all locations analysed with these virtual spacecrafts, the suprathermal ion ratio $n_{\mathit{\alpha }\text{,st}}{n}_{\text{p,st}}^{-\mathrm{1}}$ is larger than the solar wind ion ratio $n_{\mathit{\alpha }\text{,sw}}{n}_{\text{p,sw}}^{-\mathrm{1}}$, as shown already by Fig. 1b. The agyrotropy is also more pronounced for He²⁺ than for H⁺, as illustrated in Fig. 8a–b, and with the VDF discussed in the previous section.

The two bottom rows of Figs. 9 and 10 (panels g, h, o and p) display the wavelet power spectra of $\left|B\right|$ and B_y, in order to investigate foreshock wave activity. The wavelet transform (Torrence and Compo, 1998) is calculated using the Morlet wavelet. The black contours in the power spectra show the 95 % confidence level, and the cross-hatched regions bordering the spectra depict the cone of influence, where edge effects arising due to the time series' endpoints become important. The horizontal dashed line on the B_y wavelet power spectra is drawn at 50 s, which is close to the expected spacecraft frame period of foreshock fast magnetosonic waves for these upstream conditions, according to empirical models (Le and Russell, 1996; Takahashi et al., 1984). At all four positions, an enhancement in the B_y component power spectra can be seen in the vicinity of this period, which is in agreement with previous works showing that fast magnetosonic waves permeate the foreshock in Vlasiator simulations as in spacecraft observations (Palmroth et al., 2015; Turc et al., 2018). We note that none of the virtual spacecraft selected here display the typical quasi-monochromatic foreshock waves shown in these previous studies, due to their relative proximity to the bow shock (≲7 R_E). In the vicinity of the shock, the wave activity is more complex due to nuanced interactions between the waves and the diffuse ion population (Greenstadt et al., 1995; Turc et al., 2018).

At point F we find weaker wave activity, despite its location deep within the foreshock. This is most likely due to the low density of suprathermal gyrating or beam particles in this part of the foreshock (see Fig. 10j). This results in a lower wave growth rate as this parameter depends on the beam density for the beam–beam instabilities at play in the foreshock (Gary, 1993). This weaker wave activity is accompanied by a lower temperature anisotropy for the He²⁺ ions, due to the suprathermal He²⁺ population being in the form of low-energy, field-aligned beams in this region (see Fig. 7, panel 4). Comparing the positions of our virtual spacecraft with Fig. 8 and the contours indicating well-structured or more broken up ULF wave fronts, we see that the weakest wave power (point F) is seen at the most broken up position, and the strongest power (point D) is found at a position of well-structured wave fronts.

According to previous works, the plasma rest frame frequency of foreshock fast magnetosonic waves generated by proton beams is of the order of 10 % of the proton gyrofrequency, probably due to the cyclotron resonance which gives rise to the waves (Hoppe and Russell, 1982; Eastwood et al., 2005; Wilson, 2016). He²⁺ ions have a gyrofrequency that is half that of protons. He²⁺ ion beams could, thus, generate fast magnetosonic waves at a period of ∼100 s with the upstream parameters used in our simulation. The wavelet power spectra in Figs. 9 and 10 show enhanced wave power around ∼100 s, but this part of the spectra is mostly inside the cone of influence of the wavelet transform, making it difficult to draw firm conclusions. Moreover, the analysis of the foreshock wave properties in an identical Vlasiator run without a helium population (not shown) reveals similar enhancements of the wave power at ∼100 s. It is therefore unlikely that these fluctuations are due to helium beam instabilities. Thus, despite a 1 % solar wind helium content providing non-negligible mass loading, we find the helium component to not have a significant impact on foreshock wave populations.

4 Discussion and conclusions

In this study, we present how sometimes the proton foreshock extends further out upstream than the helium foreshock. Our model of the edge of the foreshock shows that, within a few R_E of the foreshock edge, the ratio of suprathermal alphas to protons (as shown in Fig. 1b) often tends towards a value of about 10 times less than the solar wind. This seems to be in excellent agreement with Fig. 1b of Fuselier and Thomsen (1992). They restricted the analysis of beam ions to high energies (1.6–5.7 times the solar wind energy), which will likely require measurements very close to the foreshock edge. At any given distance from the shock, slower particles will have travelled for a greater period of time and will thus have drifted deeper into the foreshock due to, for example, E×B drift. In the central section of the FAB region, but still outside the ULF foreshock, we report a suprathermal alpha fraction of the order of the solar wind value. At the very sunward edge of the foreshock in Fig. 1b, some regions of relatively more abundant suprathermal alphas can be seen, and we suggest that this is an effect stemming from the larger gyroradii of alpha particles, as the feature is clearly spatially periodic.

In our simulation, we find that, well inside of the upstream edge of the foreshock, there is a region with strong, well-structured ULF wavefronts and a relative decrease in suprathermal alpha particles, as shown in Fig. 1b. We suggest that this could be due to these very strong proton-driven ULF waves being strong enough to scatter even alpha particles which are not in resonance with the wave oscillation, due to a larger mass-to-charge ratio than that of protons. Further inside the foreshock, the ULF wave front breaks up into less uniform waves. This break-up is inherent to proton dynamics and not caused by the presence of alphas. Previously, the growth of waves with distance from the foreshock edge has been reported in Le and Russell (1992), but their study was performed close to the nose of the shock, whereas the destructuring we report takes place far along the flank of the bow shock. In this region we report a helium fraction higher than that of the solar wind. One possible explanation for this is that the destructured ULF waves are still able to scatter protons, but alpha particles propagate in a less disturbed fashion, more akin to a low-energy FAB. This is seen as a second region of parallel pressure enhancement for alphas in Fig. 8d.

Based on the time–energy spectrogram results in panels d and e of Figs. 3 and 4, it appears that the foreshock edge transitions for H⁺ and He²⁺ were in more agreement for the December 2018 event than for the November 2018 event. During the event depicted in Fig. 3, the IMF upstream of the foreshock was dominated by the B_y component, and MMS was located at approximately $X=\mathrm{13},Y=\mathrm{1},$ and Z=5R_E (roughly at the nose of the shock). If we assume a parabolic bow shock shape, the field lines at the foreshock edge were thus connected to a region of the shock where ${\mathit{\theta }}_{B\mathrm{n}}\sim \mathrm{45}{}^{\circ }$. For each field line, we can estimate the derivative of θ_Bn as follows:

$${\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_{B\mathrm{n}}}{\mathrm{d}{\stackrel{\mathrm{̃}}{r}}_{⟂\text{FSedge}}}},$$

where ${\stackrel{\mathrm{̃}}{r}}_{⟂\text{FSedge}}$ is a normalized spatial distance vector perpendicular to the foreshock edge. At the location of MMS with the listed IMF conditions, the derivative of θ_Bn respective to the distance to the foreshock edge is large. Conversely, during the event depicted in Fig. 4, MMS was located at approximately $X=\mathrm{10},Y=\mathrm{11},$ and Z=6R_E, i.e. somewhat at the flank, and the IMF had a strong −B_x component. Thus, the field lines at the foreshock edge can be assumed to connect to a region where the corresponding derivative of θ_Bn for each field line is small. When comparing two Vlasiator simulations with different IMF directions in Fig. 5 and corresponding qualitatively with the MMS observation situations, we see qualitatively similar behaviour. We note that the contours in Fig. 5 were averaged over time, smoothing out variations due to ion gyroradii, ion gyrotimes, and reflection variations due to mesoscale bow shock reformation. Variations such as these are detectable in virtual and real spacecraft measurements. Based on this comparison, we suggest that the different profiles of the foreshock edge transition, seen particularly well in the suprathermal ion densities, may be related to the derivative of θ_Bn for the connecting field line, as we are able to replicate the MMS observation variation in the gradient of suprathermal ion profiles using two simulation runs with different IMF cone angles. In a previous study, Sibeck et al. (2008) investigated foreshock edge gradients in radial IMF hybrid simulations, finding a correlation with model dipole tilt. They stated that their bow shock edges were still propagating outwards, and the foreshock in their Fig. 1 appears tilted to the south, possibly as a result of non-radial magnetic field components and, perhaps, a similar underlying cause as in our explanation. This can be investigated further in both observational and simulational studies.

Figure 8 shows that helium exhibits a greater agyrotropy throughout the ion foreshock. Interestingly, enhanced helium agyrotropy appears to coincide with well-structured ULF regions. This suggests that the mostly proton-induced ULF waves might be efficient drivers of agyrotropy for gyrating He²⁺ ions which have half the gyrofrequency of protons, leading them to scatter into the diffuse distribution. Proton-induced waves have been shown to heat the solar wind helium populations (Hollweg and Turner, 1978; Dusenbery and Hollweg, 1981). We propose that this additional heating by ULF waves is one reason for the helium ions exhibiting greater agyrotropy values than protons throughout the foreshock. Another source of agyrotropy and VDF break-up in helium may be that He²⁺ ions have greater gyroradii than H⁺ ions, and thus, each ion scans a greater extent of foreshock waves and structures.

We note the absence of a diffuse population of ions in our proton and helium VDFs in Figs. 6 and 7 but deduce that it is due to the low phase–space density of the diffuse ion population extending below our velocity space sparsity threshold, leading to those ions being discarded. We do not expect the diffuse ion population to strongly affect the wave dynamics, maintaining reasonable validity of the rest of our results. We do admit that some of our deductions about proton and alpha particle scatterings depend on this assumption that the suprathermal density decrease can be interpreted as the scattering of ions from a gyrating population to a diffuse population. At least some of the reported relative alpha particle abundance in Fig. 1 can be attributed to this, which agrees with the proton-induced ULF waves being more efficient at scattering protons, leading to a greater portion of the proton suprathermal population scattering into the diffuse part of velocity space, whereas alphas remain as gyrating ions tracked by the simulation (and resulting in high agyrotropies as well). We also note that our study is limited to lower energy FABs with energies of the order of 1–1.5 times the solar wind energy. The lack of a deka–keV FAB in our simulation may result from our choice of cold electrons, i.e. no electron pressure gradient term at the shock front. Since protons with a mass-to-charge ratio of 1 would feel the effect of a cross-shock potential stronger than alpha particles would, this can also partially explain the predominantly large helium fraction in our simulated foreshock.

We evaluate virtual spacecraft wavelets, looking for signatures of waves driven by helium beam instabilities, which could develop together with the proton beam instabilities in the foreshock, but we do not find convincing evidence for such waves. This result is as expected because of the low abundance of He²⁺ ions in our simulation, i.e. only 1 % of the solar wind proton number density. This number density is still enough to cause mass loading, pushing the bow shock 0.5–1.0 R_E further towards the Earth than in a comparative run without helium (comparison not shown here; for the other simulation see Blanco-Cano et al., 2018). Future simulations with a higher helium-to-proton ratio will allow further investigation of helium-driven waves in the foreshock. We note that, due to the ratio of gyrofrequencies, simulations must be run for an extended period of time in order to accurately capture potential helium-induced waves.

Providing a point of comparison with a similar model, Jarvinen et al. (2019) simulated the magnetosphere and the foreshock of Mercury with a global hybrid cloud-in-cell model, including 4 % of He²⁺ ions in solar wind. Jarvinen et al. (2019) presented a fast-mode ULF wave field at wave periods of approximately 5 s, and populations of both H⁺ and He²⁺ were present in the quasi-parallel bow shock region, with He²⁺ backscattering found to be more efficient than for H⁺. Similar to our model, they did not include an electron pressure gradient term. We note that, despite how the Hermean magnetosphere is much smaller and the impinging solar wind is quite different from that at the Earth, some properties of ULF waves at different planetary foreshocks appear to scale with respect to the interplanetary magnetic field magnitude, as illustrated by Hoppe and Russell (1982). Additionally, Fig. 1b shows proportionally enhanced suprathermal He²⁺ densities compared to H⁺, which is similar to the point result in Jarvinen et al. (2019), but we also show spatial structure in the ratio of the suprathermal populations.

In summary, we show how, for the simulated solar wind conditions, He²⁺ is a significant foreshock species which has dynamics similar to those of protons but with distinct properties, such as preferential heating by proton-induced ULF waves. Both protons and helium are found in the FAB at the foreshock edge, with beam energies decreasing from ≳ 5 to $\sim \mathrm{3}\mathrm{keV}{\mathrm{nuc}}^{-\mathrm{1}}$ when going from the foreshock edge to the inner edge of the FAB region. Helium ions are found to exhibit more agyrotropy than protons, probably due to their interaction with the proton-dominated ULF waves.

The profiles of suprathermal abundances of proton and helium at the foreshock edge show variability, both in Vlasiator simulations and in MMS data, and we show how this may be explained with the prevailing IMF orientation affecting the derivative of θ_Bn at the field-line-connected position at the bow shock.