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 ﬁeld (IMF). It is permeated by shock-reﬂected ions and electrons, low-frequency waves, and various plasma tran-sients. We investigate the extent of the He 2 + foreshock using Vlasiator, a global hybrid-Vlasov simulation. We perform the ﬁrst numerical global survey of the helium foreshock and interpret some historical foreshock observations in a global context.Theforeshock edge is populated by both proton and helium ﬁeld-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 2 + 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 reﬂection of ions at the bow shock. We show that specular reﬂection 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 2 + to be scattered less than protons.


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;, waves (Hoppe et al., 1981;Blanco-Cano et al., 2009;, and non-Published by Copernicus Publications on behalf of the European Geosciences Union. linear 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 2+ 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 −1 H + and He 2+ 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 2+ 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 2+ 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 nongyrotropic gyrating distribution. These gyrating He 2+ 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 2+ 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 nearspecularly reflected He 2+ 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 ∼ 100 keV/e. They found that the ratio of He 2+ to H densities with suprathermal energies (normalized to solar wind abundances) is dependent on the location within the foreshock. High-energy field-aligned beams (> 10 keV/e) near the foreshock edges show significant He 2+ / H ratios (near solar wind quantities), whereas lower energy beams (∼ 1 keV/e) deeper within the foreshock exhibit intermediate proton distributions and lower He 2+ / 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 2+ . Additionally, Fuselier et al. (1995) found that, still deeper within the quasi-parallel foreshock, He 2+ 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 2+ to suprathermal H ion densities is similar to that of the solar wind composition, usually n p /n α ∼ 4 % 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 highconcentration He 2+ gyrating ions are observed in the quasiparallel foreshock, as are diffuse ions, suggests that gyrating distributions can be the seed population for the diffuse He 2+ 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 2+ 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 2+ 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 2+ distribution reentering 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 2+ ions can alter the shock structure and that the occurrence of He 2+ 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  has allowed new investigations of helium ions near the Earth's bow shock, providing, in particular, He 2+ velocity distributions. Broll et al. (2018) investigated the reflection of He 2+ at the quasi-perpendicular bow shock and showed that He 2+ 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 2+ properties in the foreshock using a global hybrid-Vlasov simulation of near-Earth space performed with the Vlasiator model . 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 2+ 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.

Simulations
We investigate the Earth's foreshock region using Vlasiator , 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 ∈ (−48.66 R E ; 64.35 R E ) and Z ∈ (−59.65 R E ; 39.24 R 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 −1 .
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 p = 1 cm −3 and n α = 10 −2 cm −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 sw = 750 km s −1 in the −ê 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 = −3.54 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ê z -aligned line dipole, resulting in a realistic magnetopause standoff distance (Daldorff et al., 2014). The simulation has an inner boundary at 3 × 10 4 km ≈ 4.7 R 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 min,p = 10 −15 s 3 m −6 and f min,α = 10 −17 s 3 m −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 p = 3.3 cm −3 , n α = 3.3 × 10 −2 cm −3 , and u sw,x = −600 km s −1 . The IMF is set to 5 nT, with a 5 • cone angle.

Observations
In this study, we also analyse observations from the MMS mission  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 −1 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 2+ , O + , and O 2+ ). 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 ∼ 1 eV/e and 40 keV/e, with a mass resolution of M/ M ∼ 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.

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 = ±0.1 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 −1 , centred at the solar wind speed u sw,x = −750 km s −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.

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 ∼ −5 R E , whereas the ion foreshock edge intersects the bow shock already at Z ∼ +5 R E (Fig. 1b). At the bottom of the figure, at Z = −50 R 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 α,st n −1 p,st × 10 2 ratio, shown in Fig. 1b, tends to values 2, whereas at the foreshock edge, it falls below 0.2. This abundance of He 2+ within the deep foreshock is likely a computational artefact, with H ions being efficiently scattered into 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 = ±0.1 nT at time t = 1100 s. the diffuse population, which is not tracked, whereas He 2+ 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 2a shows an excerpt from Fig. 1a, featuring a portion of the foreshock edge. The colour map is again the outof-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 . 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.
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 GSE −Y 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 2+ (green; HPCA), and all ions (red; FPI) and proton and He 2+ energytime spectrograms from the HPCA instrument. The He 2+ 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 2+ . We note that the energy-time spectrograms show energy per charge E i Z −1 i , 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., 2008Billingham et al., , 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 2+ suprathermal densities do not always behave in the same manner. In Fig. 3c, we can clearly see that the proton and He 2+ 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 2+ density and remains low well into the foreshock. This suggests that the proton foreshock extends further outward than the He 2+ foreshock. In other words, the proton foreshock extends to field lines connected to larger θ Bn values than the He 2+ 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 ∼ 55 • . 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 (12.5, 0.8, 5.1) R 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 ∼ 35 • 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 (9.8, 11.3, 5.7) R E , i.e. far from the Sun-Earth line. Thus, this foreshock edge crossing is comparable with the one shown in Fig. 5b.

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 ( ⊥)), 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 ( ⊥) 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 −1 in the solar wind frame or u sw = 750 km s −1 in the simulation frame, which translates to roughly 3-5 keV nuc −1 . 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 = 16 R E ; Z = 0 R 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×V > 500 km s −1 ) for both species but particularly for helium at v B×V > 1000 km s −1 . Parallel velocities have begun to decrease, extending down to 1000 km s −1 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 −1 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 gyrophasebunched 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 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 ( ⊥) 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.
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 nongyrotropy. We apply the temperature anisotropy to the whole VDF, allowing us to identify regions where VDFs have Q ag = P 2 12 + P 2 13 + P 2 where P = P 11 is the parallel pressure, P ⊥ = 0.5(P 22 +P 33) is the perpendicular pressure, and P 12 , P 13 , and P 23 are offdiagonal 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 ⊥,p T −1 ,p ≈ 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.
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 = 10 R E ; Z = −40 R E ). Referring to Fig. 1b, we see that this band of lowenergy, 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 (α ⊥⊥) 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 = −25 R 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.

Foreshock waves
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 ⊥ T −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 |B| and its out-of-plane component B y , respectively.
The fluctuations of H + and He 2+ total densities follow each other quite closely at all locations, though the amplitude of the He 2+ density oscillations is larger than that of protons at point D and E. These larger He 2+ density variations seem to be well correlated with the fluctuations of the He 2+ 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 α,st n −1 p,st is larger than the solar wind ion ratio n α,sw n −1 p,sw , as shown already by Fig. 1b. The agyrotropy is also more pronounced for He 2+ 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 |B| 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 mod- els (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 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 2+ ions, due to the suprathermal He 2+ 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;. He 2+ ions have a gyrofrequency that is half that of protons. He 2+ 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.

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 2+ 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 = 13, Y = 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 θ Bn ∼ 45 • . For each field line, we can estimate the derivative of θ Bn as follows: wherer ⊥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 = 10, Y = 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 2+ ions which have half the gyrofrequency of protons, leading them to scatter into the diffuse distribution. Protoninduced 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 2+ 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 crossshock 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 2+ 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 2+ 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 2+ were present in the quasi-parallel bow shock region, with He 2+ 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 2+ 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 2+ 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 ∼ 3 keV nuc −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.
Code and data availability. Vlasiator (https://www.helsinki.fi/en/ researchgroups/vlasiator; Palmroth, 2020) is distributed under the GPL-2 open source license at https://github.com/fmihpc/vlasiator/, https://doi.org/10.5281/zenodo.3640593 (Palmroth and the Vlasiator team, 2020). Vlasiator uses a data structure developed inhouse at the University of Helsinki, and it is available at (https: //github.com/fmihpc/vlsv/, Sandroos, 2019). The Analysator software is available at (https://github.com/fmihpc/analysator/; Hannuksela and the Vlasiator team, 2020) and was used to produce the presented figures. The run described here takes several terabytes of disk space and is kept in storage maintained within the CSC -IT Center for Science. Data presented in this paper can be accessed by following the data policy on the Vlasiator website.
Video supplement. The supplementary videos A, B, C, and D provide movie extensions of Figs. 1b, 6, 7, and 8, showcasing the temporal evolution of foreshock features and particle population shapes and properties. Movie A (Battarbee, 2020a) is a filmed extension of Fig. 1b. The animation of the foreshock region of the simulation includes over 250 s of simulation. The ratio of suprathermal densities for helium over protons is normalized to the solar wind ratio of 1 %. The ratio is not shown in the pristine solar wind. Black contours are drawn for B y = ±0.1 nT. Movie B (Battarbee, 2020b) is a filmed extension of Fig. 6. It shows the animation of velocity distribution functions for protons and alpha particles and their locations in the foreshock. Panel 3 shows a 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 (namely A and B) and one just within the ULF foreshock (C). In each set, the subpanels are labelled as v B versus v B×(B×V ) ( ⊥; left columns) or v B×V versus v B×(B×V ) (⊥⊥; right columns) and protons (top rows) or helium (bottom rows). Movie C (Battarbee, 2020c) is a filmed extension of Fig. 7. It shows an animation of the velocity distribution functions for protons and alpha particles and their locations in the foreshock. Panel 3 shows a 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 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×(B×V ) ( ⊥; left columns) or v B×V versus v B×(B×V ) (⊥⊥; right columns) and protons (top rows) or helium (bottom rows). Movie D (Battarbee, 2020d) is a filmed extension of Fig. 8. It shows an animation of the foreshock characteristics for proton and helium over 363 s of simulation. The left column (panels a and c) shows protons, and the right column (panels b and d) shows helium. The top row shows agyrotropy Q ag (see Eq. 1), where a value of 0 corresponds with perfect gyrotropy. The bottom row shows temperature anisotropy T ⊥ T −1 , calculated for the total distribution function including both solar wind and suprathermal parts. Black contours show the out-ofplane magnetic field B y fluctuations at a level of ±0.2 nT, indicating the extent of the ULF foreshock region.
Author contributions. This paper was outlined and drafted at the Third International Vlasiator Science Hackathon held in Helsinki, Finland, on 19-23 August 2019. MB, XBC, LT, and PK performed the preliminary investigation leading to the team effort. MB led the investigation and the writing process. AJ provided the MMS plots and SF and KT verified them. VT performed the wavelet analysis with the help of LT. MP led the Vlasiator team and organized the hackathon. UG and YPK were instrumental in adding helium support to Vlasiator. All coauthors (including MA, TB, MAT, TK, SR, MD, and JS) provided scientific input during the Hackathon event, during the pre-submission review process, or both.