Hybrid-Vlasov simulation of auroral proton precipitation in the cusps: Comparison of northward and southward interplanetary magnetic field driving

Particle precipitation is a central aspect of space weather, as it strongly couples the magnetosphere and the ionosphere and can be responsible for radio signal disruption at high latitudes. We present the first hybrid-Vlasov simulations of proton precipitation in the polar cusps. We use two runs from the Vlasiator model to compare cusp proton precipitation fluxes during southward and northward interplanetary magnetic field (IMF) driving. The simulations reproduce well-known features of cusp precipitation, such as a reverse dispersion of precipitating proton energies, with proton energies increasing with increasing geomagnetic latitude under northward IMF driving, and a nonreversed dispersion under southward IMF driving. The cusp is also found more polewards in the northward IMF simulation than in the southward IMF simulation. In addition, we find that the bursty precipitation during southward IMF driving is associated with the transit of flux transfer events in the vicinity of the cusp. In the northward IMF simulation, dual lobe reconnection takes place. As a consequence, in addition to the high-latitude precipitation spot associated with the lobe reconnection from the same hemisphere, we observe lower-latitude precipitating protons which originate from the opposite hemisphere’s lobe reconnection site. The proton velocity distribution functions along the newly closed dayside magnetic field lines exhibit multiple proton beams travelling parallel and antiparallel to the magnetic field direction, which is consistent with previously reported observations with the Cluster spacecraft. In both runs, clear electromagnetic ion cyclotron waves are generated in the cusps and might further increase the calculated precipitating fluxes by scattering protons to the loss cone in the low-altitude cusp. Global kinetic simulations can improve the understanding of space weather by providing a detailed physical description of the entire near-Earth space and its internal couplings.


Introduction
The precipitation of particles from near-Earth space into the upper atmosphere is one of the sources of disruption of telecommunication and navigation signals (Smith et al., 2008;Jin et al., 2017), and carries field-aligned currents which find their closure in the ionosphere and may lead to Joule heating in the lower thermosphere (e.g., Sarris et al., 2020) and to the formation of geomagnetically induced currents in the ground, threatening power grids (e.g., de Villiers et al., 2017). Particle precipitation can also durably affect the chemical composition of the upper atmosphere, leading for instance to ozone destruction and HO x and NO x species formation in the stratosphere and mesosphere (e.g., Andersson et al., 2014;Seppälä et al., 2015). As such, it is a major aspect of the space weather chain, linking in particular the magnetosphere (perturbed by solar driving) to the atmosphere and ionosphere. Yet, it remains particularly challenging to model it using first-principle numerical codes and to predict its energy spectrum and flux under a given set of solar wind parameters. Particle precipitation predominantly takes place in the auroral ovals and hence particularly affects the high latitudes, both on the dayside and on the nightside, for which not only observations but also numerical simulations are needed to improve space weather forecasting (Robinson et al., 2019;Heelis & Maute, 2020).
In the dayside magnetosphere, the polar cusps are the two locations where plasma from the magnetosheath has direct access to the Earth's atmosphere and ionosphere. Their footprints in the ionosphere are therefore regions where particle precipitation is frequent, even when the overall geomagnetic activity is low (Hardy et al., 1985(Hardy et al., , 1989. Precipitation of both protons and electrons at auroral energies (up to a few tens of keV) takes place in the cusps and leads to dayside auroral emissions mostly in the Oxygen 630.0 nm (red) and 557.7 nm (green) lines (e.g., Jacobsen et al., 1990), but also UV emissions (Doppler-shifted Ly-a emission line) associated with proton precipitation (e.g., Frey et al., 2002). Recently, Frey et al. (2019) provided an extensive review of dayside aurora, discussing in particular cusp aurora during various types of solar wind driving, complementing their previous review (Frey, 2007).
It is well-established that the cusp location has a strong dependence on the interplanetary magnetic field (IMF) components in the geocentric solar magnetospheric (GSM) system. B y has an effect on the magnetic local time (MLT) location and B z on the geomagnetic latitude, whereas the intensity of the precipitation is controlled by the solar wind dynamic pressure (e.g., Frey et al., 2002Frey et al., , 2003. The ideal picture is such that, during southward IMF driving, the cusp is typically located at about 58-70°geomagnetic latitude (Frey, 2007), and its equatorward part is magnetically connected to the dayside magnetopause, where reconnection takes place. On the other hand, during northward IMF driving, the cusp is generally observed at much higher geomagnetic latitudes (75-85°; Frey, 2007), as magnetic reconnection takes place at the lobes. While most of the time single lobe reconnection is observed, dual lobe reconnection can occur in case the IMF direction is almost exactly northwards (clock angle within about ±10°, Imber et al., 2006). This leads to the formation of new closed field lines in the dayside magnetosphere. Early studies have also shown that the cusp footpoint can migrate to a lower latitude during geomagnetic storms (Meng, 1982). Observations of particle precipitation have also enabled the estimation of the longitudinal extent of the cusp, which is typically covering about 2.5 h in local time and centred around noon (Newell & Meng, 1992).
Cusp precipitation can be observed directly, with particle detectors on spacecraft flying in the cusp region, and indirectly through the auroral emissions and ionospheric plasma remote sensing. One example of satellite missions comprising particle instruments suitable for auroral precipitation studies is the Defense Meteorological Satellite Program (DMSP) fleet, carrying the Special Sensor J (SSJ) instrument measuring precipitating <30 keV proton and electron fluxes, from which a new database was recently published (Redmon et al., 2017). Other past and present spacecraft measuring particle fluxes include the Fast Auroral Snapshot Explorer (FAST) (Carlson et al., 2001) on low-Earth orbit (LEO), but also spacecraft on wider orbits exploring the near-Earth space, such as those of the Cluster mission (Escoubet et al., 2001), which include the composition distribution function (CODIF) instrument (Rème et al., 2001) and had dedicated cusp observation campaigns. A statistical study of Cluster cusp crossings revealed for instance that it takes over 20 min for the cusps to adjust to changes in the solar wind driving, and that the motion and shape alteration of the cusp is facilitated in the sunlit hemisphere (Pitout et al., 2006). Indirect observations of cusp particle precipitation can be achieved by measuring the auroral emissions. This can be done from above the ionosphere, with imagers aboard spacecraft such as the Far Ultraviolet (FUV) instruments on the Imager for Magnetopause-to-Aurora Global Exploration (IMAGE) satellite (Mende et al., 2000), the Spectral Sensor Ultraviolet Spectrographic Imager (SSUSI) on DMSP, or Global Ultraviolet Imager (GUVI) on Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) (Paxton et al., 2003). The advantage of observing the Doppler-shifted Ly-a emission associated with ion precipitation is that it allows measurements even over sunlit regions, which enables the study of both cusps at the same time (e.g., Østgaard et al., 2005). Cusp aurora can also be observed with suitably located ground-based optical instruments, as well as radars. An example of multi-instrument study of cusp precipitation can be found in Lorentzen et al. (2007), who combined observations from the European Incoherent Scatter (EISCAT) Svalbard Radar, optical instruments from the Kjell Henriksen Observatory, a DMSP satellite and a rocket experiment to explain a pulsating dayside aurora event by pulsed lobe reconnection driven by IMF with positive B z and strongly negative B y . Coherent scatter radars can also provide information on the cusp dynamics by providing information on the plasma horizontal convection; the SuperDARN radar network has proven useful for such studies (e.g., McWilliams et al., 2001).
Cusp precipitation as a function of IMF driving has been the topic of many studies. With southward IMF conditions, precipitation is often observed as bursty (e.g., Mende et al., 1990;Lockwood et al., 1993) and poleward-moving auroral forms can be seen. These suggest that the cusp precipitation during southward IMF driving may be the result of pulsed reconnection at the dayside magnetopause (Øieroset et al., 1996). In addition, precipitating ions exhibit an energy dispersion, with energies decreasing with increasing geomagnetic latitude (e.g., Shelley et al., 1976;Reiff et al., 1977;Smith & Lockwood, 1996;Pitout et al., 2009). This can be explained by a velocity filter effect, as ions with larger parallel velocities reach the ionosphere earlier than those having a lower parallel velocity, which will reach the Earth after the field line on which they are locatedwhere the particle injection has taken placehas convected further polewards (Burch et al., 1982). On the other hand, a typical signature observed during northward IMF driving is a high-latitude aurora spot visible in both the proton and electron aurora emission wavelengths (e.g., Sandholt et al., 1996;Milan et al., 2000). This aurora spot is believed to be the ionospheric footprint of the lobe reconnection site. Spacecraft flying through this spot observe a reverse-dispersion signature in precipitating proton energies, that is, increasing energies with increasing geomagnetic latitude (Øieroset et al., 1997).
Numerical simulations of cusp precipitation using firstprinciple models are scarce. One noteworthy example is the study by Connor et al. (2015), which couples a resistive magnetohydrodynamics (MHD) modelthe Open Global Geospace Circulation Model (OpenGGCM)with the Liouville Theorem Particle Tracer (LTPT). These codes have been used to study the precipitation signatures in the cusp as a function of the IMF clock angle. In particular, they showed that a purely northward IMF results in reverse-dispersion signatures in ion cusp precipitation, associated with antiparallel reconnection at the lobes. On the other hand, for purely southward IMF, their simulations showed a "normal" dispersion signature caused by antiparallel reconnection at the dayside magnetopause.
The coupling of MHD and kinetic codes was also used in studies by Berchem et al. (2014Berchem et al. ( , 2016, who ran global resistive MHD simulations providing the electromagnetic fields for large-scale test-particle simulations. Berchem et al. (2014) revealed dawn-dusk asymmetries in the dayside precipitation, whereas Berchem et al. (2016) analysed the effect of the rotation of the IMF from purely southwards to predominantly northwards on the asymmetries of the dayside proton precipitation between the northern and the southern hemispheres.
The main limitation of MHD simulations associated with test particles is that the fields and the plasma are not coupled self-consistently. Kinetic models are a way to circumvent this issue and take into account the feedback from the plasma to the electromagnetic fields. A study by Omidi & Sibeck (2007) used a hybrid-particle-in-cell (hybrid-PIC) model to study the interactions between flux transfer events (FTEs) and the cusps during southward IMF driving. They found that FTEs can inject downward ions into the cusp and produce signatures that are consistent with poleward-moving auroral forms in the ionosphere. In a three-dimensional simulation study using a hybrid-PIC model, Wang et al. (2009) linked the energetic (10 keV-1 MeV) protons observed in the cusps to the quasiparallel bow shock and foreshock and concluded that such high energies could be reached as the result of Fermi acceleration. Additional self-consistent global kinetic simulations are needed to further study ion precipitation into the cusps under various types of IMF driving, in particular to characterise the precipitating particle differential and integral fluxes and relate them to, e.g., magnetic reconnection, FTEs, or plasma waves. This paper presents a comparison of auroral proton precipitation in the cusps during northward IMF and southward IMF driving in self-consistent kinetic (hybrid-Vlasov) simulations with the Vlasiator model. Precipitating proton differential number fluxes, integral energy fluxes and mean precipitating energy are calculated and analysed in light of the global near-Earth space dynamics context, including the transit of FTEs, magnetic reconnection and wave activity. Section 2 gives an overview of the Vlasiator model and the runs used in this study, and briefly describes the method used to derive proton precipitation fluxes. Section 3 presents the results obtained in this study, which are discussed in Section 4. Finally, a summary of the paper is given in Section 5.

Data and Methods
This study is based on two simulations of near-Earth space using the Vlasiator global hybrid-Vlasov model (von Alfthan et al., 2014(von Alfthan et al., , 2020Palmroth et al., 2018). In Vlasiator, protons are represented as distribution functions that are discretised in a cartesian 2D ordinary-space grid and in a cartesian 3D velocityspace grid, hence giving a 5D representation of the phase-space density. The time evolution of the proton phase-space density follows the Vlasov equation. Electrons are described as a cold, massless, charge-neutralising fluid. The electromagnetic fields are self-consistently propagated using Faraday's law, Gauss's law, and the solenoidal condition. Closure of the system is achieved with the generalised Ohm's law including the Hall term. While there is no general consensus regarding the best way to include the electron pressure gradient term in the generalised Ohm's law, with authors assuming, e.g., an adiabatic electron fluid (e.g., Giacalone, 2004;Omidi et al., 2014;Hao et al., 2016), an isothermal electron fluid (e.g., Ofman & Gedalin, 2013), polytropic electrons (e.g., Caprioli & Spitkovsky, 2014), or a fixed electron-to-ion temperature ratio (e.g., Yang et al., 2011), this term is neglected in Vlasiator simulations as well as in some other hybrid models (e.g., Gargaté et al., 2007;Bai et al., 2015;Sundberg et al., 2016). Hoilijoki et al. (2017) showed that despite this approximation, the reconnection rates obtained with Vlasiator are consistent with theoretical values given by Cassak & Shay (2007), which suggests that neglecting the electron pressure gradient term when focusing on ion-scale phenomena is not unreasonable. Similarly, Yang et al. (2009) evaluated that, at ion scales, the Lorentz and Hall terms are the main contributors to Ohm's law, while the electron pressure gradient term is negligible, even in a region with sharp discontinuities such as the shock ramp.
Both Vlasiator simulations considered in this study are in the noon-midnight meridional plane, i.e., in the XZ plane in the geocentric solar ecliptic (GSE) coordinate system, assuming zero dipole tilt. We use a scaled line dipole field for the geomagnetic dipole Daldorff et al. (2014) due to the 2D geometry. In the first simulation run, coined southward IMF run in the following, the simulation domain extends from X = À96 R E in the nightside to X = +47 R E in the dayside and from Z = À56 R E to Z = +56 R E in the north-south direction, where R E = 6371 km is the Earth radius. The inner boundary of the simulation domain is located at 4.7 R E from the centre of the Earth. The second simulation run, hereinafter coined northward IMF run, has the same domain extents, except that the nightside boundary lies at X = À47 R E . Both runs have the same boundary conditions: the northward, southward and nightside walls have Neumann conditions, the inner boundary is a perfectly conducting sphere, and periodic conditions are implemented in the out-of-plane direction (±Y). The ordinary space resolution is 300 km, and the velocity space resolution is 30 km/s. From the dayside boundary, a homogeneous and stationary solar wind inflow drives the near-Earth plasma dynamics throughout the simulations. In both simulations, the solar wind protons have a Maxwellian velocity distribution with a bulk speed of 750 km/s along the ÀX direction, a density of 1 proton per cubic centimetre, a dynamic pressure of about 0.9 nPa and a temperature of 500 kK. The only difference between the two simulations is the IMF, which is purely southwards in the first run and purely northwards in the second run, hence the chosen terminology. In both cases, the IMF magnitude is 5 nT. These solar wind parameters give an ion inertial length k p = 228 km, a proton thermal gyroradius r L = 190 km, and an Alfvén speed V A = 109 km/s. While the ordinary space resolution of 300 km is greater than both k p and r L , Pfau-Kempf et al. (2018) showed that it is nevertheless sufficient to resolve most of the kinetic physics for the protons in regions such as the shock and the foreshock. In the context of this study, the spatial resolution needs to be sufficient to allow the emergence of the waves that can interact with protons and in particular affect the loss-cone population. As will be seen in Section 3.4, EMIC waves do appear in the simulations, which indicates that despite not fully resolving all ion scales the model reproduces the relevant kinetic physics to study proton precipitation. Figure 1 shows a zoomed-in box of the southward IMF simulation at t = 1600 s, centred on the dayside magnetosphere. The background colour shows the proton temperature; the area in white corresponds to the near-Earth space located Earthwards from the inner boundary of the simulation domain. The black lines indicate magnetic field lines. Magnetic reconnection is taking place near the subsolar point and two magnetic islands are present at the dayside magnetopause. Magnetic islands and cusp plasma are particularly easy to identify with enhanced proton temperature values. We note that the dayside magnetospheric temperature is unphysically low in this run, which is an effect of the 2D simulation setup; this however is not expected to affect the results presented in this paper, since the data analysis is performed outside the magnetosphere. The magnetopause is nevertheless at a realistic location, its subsolar point being at a distance of about 8.5 R E from the centre of the Earth. In each cusp, the location marked with a black circle indicates the simulation cell where precipitating fluxes will be derived in the continuation of the paper, which will also be called "virtual spacecraft" in the following. For each of these cells, the corresponding inset panel shows a cut of the velocity distribution function (VDF) in the (v x , v z ) frame. The loss cone is indicated in magenta for precipitation into the northern cusp and in brown for precipitation into the southern cusp; this convention will be kept consistent throughout the paper. The half loss cone angle values are of the order of 2°. Supplementary Animation S1 consists of frames built in a similar way as Figure 1 from t = 1350 s to the end of the simulation at t = 2150 s every 0.5 s. This southward IMF run has been previously used to study dayside magnetosphere reconnection and flux transfer events (Hoilijoki et al., , 2019, magnetotail reconnection Juusola et al., 2018a), nightside proton precipitation (Grandin et al., 2019), magnetotail current sheet flapping (Juusola et al., 2018b), and magnetosheath ion acceleration (Jarvinen et al., 2018). Figure 2 presents a similar zoomed-in box for the northward IMF simulation at t = 1450 s. Here, magnetic reconnection takes place at the lobes rather than at the dayside magnetopause because of the different IMF orientation; in the southern lobe, a magnetic island can be identified. In a similar manner as in Figure 1, five inset panels show the cuts in the (v x , v z ) frame of the VDFs at five chosen locations indicated with black circles, with the loss cone indicated in magenta and brown for precipitation in the northern and southern cusp, respectively. In the continuation of the paper, the two locations near X = 0 will be referred to as northern and southern cusp spot, whereas the three locations on the dayside will be referred to as dayside northern, subsolar and southern locations. An animated version of Figure 2 is provided in the Supporting Information as Supplementary Animation S2, from t = 1100 s till the end of the simulation at t = 1938 s, with one frame every 0.5 s.
In the following, precipitating proton fluxes will be calculated at the locations shown with black circles in Figures 1 and 2 using the methodology presented in Grandin et al. (2019), where it was applied to the study of nightside auroral proton precipitation and its relation with dipolarising flux bundles. In short, the method consists in averaging the phasespace density f at a given location r within the loss cone and at energy E to compute the differential number flux of precipitating protons as where v ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2E=m p p is the proton speed in the velocity space corresponding to energy E, m p is the proton mass, h is the proton pitch angle, h lc is the loss cone angle, and the brackets denote averaging. From the differential number flux, the integral energy flux and the mean precipitating energy can be calculated by integration and averaging, respectively.  As was the case in Grandin et al. (2019), given that the VDFs are saved only every 50 cells in the X and Z directions to keep the memory space needed to save the output files manageable, only the few locations shown in Figures 1 and 2 are suitable to derive precipitating fluxes in the cusps.

Cusp precipitation under southward IMF
Considering first the southward IMF simulation, we calculate the differential number fluxes of protons precipitating into the northern and southern cusps by analysing the VDFs at the two locations indicated with black circles in Figure 1 between t = 1350 s and the end of the simulation at t = 2150 s. Figure 3 shows the differential number flux (Fig. 3a) and integral energy flux (Fig. 3b) of precipitating protons into the northern cusp as a function of time during the simulation, and the same for the southern cusp precipitation (Figs. 3c and 3d, respectively). The solid black lines in Figures 3a and 3c indicate the mean precipitating energy as a function of time. The time intervals with grey shading correspond to those during which the corresponding virtual spacecraft is not magnetically connected to the inner boundary of the simulation, generally because of the presence of a magnetic island drifting past the virtual spacecraft location (see Supplementary Animation S1). The calculated proton fluxes at the virtual spacecraft during those time intervals therefore do not strictly speaking correspond to precipitating fluxes, due to the magnetic topology. However, similar plasma in nearby cells is still likely to be magnetically connected to the cusp at these times.
The proton precipitation appears to be essentially bursty in both the northern and southern cusps. Most precipitation bursts start with a sudden enhancement in the differential number flux across a wide range of energies (Figs. 3a and 3c), with peak flux values of the order of 10 3 -10 4 protons cm À2 s À1 sr À1 eV À1 . The maximum precipitating energy in a given burst is generally contained within 10-30 keV initially, and then decreases to about 3-5 keV. The mean precipitating energy follows the same trend as the maximum energy, reaching 5-10 keV at the peak and decreasing down to below 1 keV. The integral energy flux also exhibits burst signatures (Figs. 3b and 3d), with maximum flux values shortly after the beginning of the bursts exceeding 10 9 keV cm À2 s À1 sr À1 , which is more than three orders of magnitude above the inter-burst values of~10 6 keV cm À2 s À1 sr À1 . The burst durations, when measured in the integral energy flux time series, are of the order of tens of seconds (20-70 s for the majority of them). Time intervals with less clear burst signatures also exist, for instance during t = 1750-1880 s for the northern cusp and during t = 1980-2150 s for the southern cusp.
In terms of flux values and precipitating energies, the orders of magnitude are similar (as given above) between the northern and southern cusps. However, it is clear that the cusp precipitation is not symmetric, even though the IMF is steady and purely southwards and the geomagnetic dipole axis is not tilted. Looking for instance at the time interval between t = 1350 s and t = 1700 s, one can see that precipitation bursts do not take place at the same time into the northern and southern cusps.
In addition, during this time interval, the southern cusp precipitation bursts have longer duration (~60-70 s) than the northern cusp precipitation bursts (~10-30 s). However, the opposite situation takes place later on during the simulation (between t = 1850 s and t = 2000 s), so this is not a systematic trend but it rather simply illustrates the asymmetry between both cusps. In a study of flux transfer events formed at the dayside magnetopause in this same Vlasiator run, Hoilijoki et al. (2019) also noted an asymmetry between the northern and southern hemisphere, with time intervals favouring more frequent and smaller flux transfer events in the northern hemisphere than in the southern hemisphere, and other time intervals with the opposite situation (see their Fig. 5a).
Focusing on the areas shaded in grey in Figure 3, corresponding to time intervals when the virtual spacecraft are not magnetically connected to the inner boundary of the simulation domain (and hence to the Earth), one can see that these are always coincident with a burst in the differential number flux and integral energy flux. Since those bursts present similar features as those which occur while the virtual spacecraft are indeed magnetically connected to the Earth, it is reasonable to conclude that the precipitation bursts are associated with the passage of flux transfer events in the vicinity of the cusps. This is indeed confirmed by comparing Figure 3 with Supplementary Animation S1, as a magnetic island can be identified near the virtual spacecraft for almost each precipitation burst. In the case of the few exceptions (for instance the burst associated with the northern cusp starting at t = 1442 s), the magnetic field line topology suggests the presence of a very small magnetic island not clearly visible with the chosen number of displayed field lines. This suggests that, in this simulation, cusp precipitation bursts are always associated with FTEs, although in some cases they are too small to be clearly visible in the animation.

Cusp and LLBL precipitation under northward IMF
We now consider dayside precipitation in the northward IMF simulation. The five virtual spacecraft are located in such a way that several aspects of dayside proton precipitation can be studied. We will first study the three dayside virtual spacecraft, and then the two virtual spacecraft located near X = 0, hereinafter coined northern and southern cusp spot virtual spacecraft. As can be seen in Figure 2, the dayside spacecraft are magnetically connected to the low-latitude boundary layer (LLBL) equatorwards of both the northern and southern cusps at t = 1450 s. It is therefore possible to consider proton precipitation in the LLBL equatorwards of both cusps as observed at each of those three locations, and indeed, each of the three associated VDFs does exhibit protons in both loss cones (towards northern cusp and towards southern cusp). At the beginning of the displayed time interval, none of the three dayside virtual spacecraft are magnetically connected to the northern cusp/LLBL. Around t = 1125 s, northern lobe reconnection starts affecting the field lines on which the dayside northern and southern spacecraft are located, implying that after that time, protons appearing inside the northward loss cone (magenta in Supplementary Animation S2) of the VDFs at these spacecraft will precipitate into the LLBL equatorwards of the northern cusp (as indicated by the ending of the grey shading in Figs. 4a and 4c). The same becomes true for the dayside subsolar virtual spacecraft around t = 1190 s (ending of the grey shading in Fig. 4b). Around t = 1150 s, the southern lobe reconnection starts affecting the magnetic field lines on which the three virtual spacecraft are located, hence injecting and accelerating protons northwards. At the dayside southern spacecraft, a field-aligned beam of protons appears at t = 1165.5 s, which is visible both in the corresponding inset panel of  M. Grandin et al.: J. Space Weather Space Clim. 2020, 10, 51 Supplementary Animation S2 and in Figure 4c. It coincides with the arrival of a plasma exhibiting enhanced temperature at the corresponding virtual spacecraft, originating from the southern lobe reconnection site. The mean energy of these precipitating protons is~30 keV. A similar signature appears at the dayside subsolar virtual spacecraft at t = 1240.0 s and at the dayside northern virtual spacecraft at t = 1274.0 s, with the samẽ 30 keV mean precipitating energy, indicating the northward propagation of precipitating protons along the dayside magnetopause.
As the southern lobe reconnection continues to inject protons on the field lines where the dayside virtual spacecraft are located, the northward loss cone remains filled with proton beams. Contrary to the southward IMF case, the precipitating energies are contained within a relatively narrow range of values, with very few particles of energies below 5 keV (except at the dayside southern virtual spacecraft; see Fig. 4c), and the precipitation is less bursty. During the time interval when most of the precipitation is seen in the northern LLBL (roughly t = 1160-1730 s), the precipitating flux characteristics remain relatively constant, except for (i) small fluctuations likely associated with variations in the southern lobe reconnection rate and (ii) a slight overall decreasing trend in mean energy. The latter might be due to the lobe reconnection slowing down as the magnetic tension which had built up prior to its onset is released. The range of precipitating energies is mainly contained within 5-40 keV, the mean precipitating energy remains of the order of 10-20 keV most of the time, and the integral energy flux fluctuates around 10 7 keV cm À2 s À1 sr À1 at the southern location and around 10 8 keV cm À2 s À1 sr À1 at the northern and subsolar locations. Since the three virtual spacecraft, whose locations were determined by the constraint that full VDFs are saved to output files in a limited number of cells (see Sect. 2), are not exactly located on the same field line, part of the differences from one spacecraft to another (besides the delay associated with propagation time from the reconnection site) may be explained by the fact that they do not observe exactly the same precipitating proton beam. However, the precipitating fluxes are very similar to one another, especially at the northern and subsolar virtual spacecraft.
A similar analysis can be carried out for precipitation into the LLBL equatorwards of the southern cusp. Figure 5 shows the proton precipitation fluxes calculated when considering the loss cone associated with the southern LLBL (in brown in Fig. 2 and in Supplementary Animation S2) at the three dayside virtual spacecraft, in the same format as for Figure 4.
Since reconnection is initiated earlier in the northern lobe than in the southern lobe in the simulation, fluxes of protons propagating towards the southern LLBL are detected earlier at the three dayside virtual spacecraft: at t = 1139 s, t = 1214 s and t = 1189 s at the dayside northern, subsolar and southern virtual spacecraft, respectively. The fact that precipitating protons are detected at the southern spacecraft before the subsolar one is likely due to the fact that the subsolar virtual spacecraft lies on magnetic field lines that are sunwards from those on which the other two dayside spacecraft are located, and that are hence affected by reconnection at the northern lobe a bit later on in the simulation.
Precipitating protons are observed at the three dayside virtual spacecraft almost constantly until the end of the simulation. Differential number fluxes have similar values of 10 2 -10 3 protons cm À2 s À1 sr À1 eV À1 and energy ranges (essentially 5-40 keV, except for the dayside northern spacecraft which also exhibits energies lower than 5 keV; see Fig. 5a) as for the northern cusp precipitation. The same is also true for the mean precipitating energy, mostly of the order of 10-20 keV, as well as for the integral energy flux values, which are around 10 7 keV cm À2 s À1 sr À1 at the northern virtual spacecraft and around 10 8 keV cm À2 s À1 sr À1 at the subsolar and southern spacecraft. The energy dispersion of protons is also visible. It should be noted that, although the three virtual spacecraft are not magnetically connected to the southern cusp throughout the end of the simulation, this does not necessarily imply that the observed fluxes are not precipitating into the cusp/LLBL. Indeed, as there is a delay between the time the protons are observed at a given virtual spacecraft and the time when they reach the southern lobe region, the field line on which they are trapped may have time to reconnect at the southern lobe before they reach the southern lobe region (see Supplementary Animation S2 after t % 1750 s).
Overall, the characteristics of proton precipitation fluxes follow the same trend between the northern and southern LLBL. In particular, in both situations, the differential number fluxes include protons with energies below 5 keV only at the virtual spacecraft located closest to the reconnecting site, i.e., at the northern (southern) spacecraft for southern (northern) LLBL precipitation. Moreover, the mean precipitating energy and the integral energy flux are also both lower at that virtual spacecraft, compared to the other two (i.e., subsolar and closest to the LLBL where protons precipitate). This suggests that precipitating protons are energised as they propagate towards the subsolar point. In the absence of additional virtual spacecraft along the precipitating proton path, it is difficult to find the cause for this energisation of the loss-cone population. Finally, it is clear from Figures 4 and 5 as well as from Supplementary Animation S2 that, in the simulation, northern and southern cusp/LLBL precipitation coming from the opposite hemisphere's lobe reconnection site is observed concurrently at the three dayside virtual spacecraft over the period t % 1275-1500 s. The associated velocity distribution functions exhibit during those times several proton beams propagating parallel and antiparallel to the magnetic field, hence bearing clear signatures of dual lobe reconnection.
We now focus on the two virtual spacecraft located near X = 0 in Figure 2, also referred to as cusp spot virtual spacecraft. Figure 6 shows the precipitating fluxes of protons at those cusp spot spacecraft, with similar conventions as in Figures 3-5.  Figures 6a-6b give the differential number flux, mean precipitating energy and integral energy flux of proton precipitation into the northern cusp spot, whereas Figures 6c-6d give these quantities for proton precipitation into the southern cusp spot. One can note that prior to t % 1500 s, the loss cone of the velocity distribution functions at these two virtual spacecraft is empty; hence the displayed time interval in this figure spanning from t = 1450 s until the end of the simulation.
It is obvious, although not unexpected, that high-latitude cusp spot precipitation does not exhibit symmetry between the northern and southern hemispheres; the fact that temporal variations in the reconnection take place independently at the northern and southern lobes is a reasonable explanation for this asymmetry. The precipitating fluxes show a significant variability, with time intervals of low flux (<100 protons cm À2 s À1 sr À1 eV À1 ) of low-energy (<3 keV) precipitation only with a mean energy below 1 keV and other time intervals with very high flux (!10 4 protons cm À2 s À1 sr À1 eV À1 ) of precipitating protons across a wide range of energies (~1-30 keV) with a mean energy of the order of 10 keV. A similar variability can be seen in the integral energy fluxes, which fluctuate from roughly 10 4 to 10 9 keV cm À2 s À1 sr À1 at both high-latitude cusp spots.
Due to the magnetic field topology, protons precipitating into the northern (southern) cusp spot originate from the northern (southern) lobe reconnection, in contrast to the situation for the lower-latitude cusp precipitation discussed above. Besides, although some flux transfer events can be identified at the lobes in this northward IMF simulation (e.g., at the southern lobe around t = 1700 s), these do not seem to be the main drivers of the proton precipitation. Rather, the precipitation seems to be associated with the plasma directly coming from lobe reconnection. This interpretation is supported by Figure 7 which shows the z (north/south) component of the proton bulk velocity in the northward IMF run at t = 1938 s. The northern cusp spot virtual spacecraft is located in the exhaust jet of the northern lobe reconnection site, and the plasma exhibits a southward bulk velocity between the reconnection site and the inner boundary. A symmetric situation can be seen at the southern cusp spot virtual spacecraft situated in the exhaust jet of the southern lobe reconnection site, with plasma showing a northward motion from the reconnection site to the inner boundary.

Cusp extents in the northward and southward IMF simulations
As mentioned in Section 2, Vlasiator output files contain the full velocity distributions only at a few regularly spaced locations in order to limit the file size, which constrains the possible locations for virtual spacecraft in this study. However, in each of the two simulations, there exists a special file, called "restart", containing full distribution functions in each ordinary space cell. In the southward IMF simulation, the restart file is at t = 2119 s, and in the northward IMF simulation it is at t = 1938 s, which is also the last time step of the simulation. This therefore enables us to produce maps of the dayside proton fluxes to obtain a more global view on the cusp precipitation at those specific time steps. Figure 8 presents the integral energy flux (Fig. 8a) and mean energy (Fig. 8b) of protons propagating along the magnetic field direction on the dayside at t = 2119 s in the southward IMF simulation. Pitch angles a < 5°are considered in the northern hemisphere, and a > 175°in the southern hemisphere; both hemispheres are separated with a white dashed line in the figures. Here, a standard "loss-cone" angle of 5°is used everywhere to track protons propagating towards the closest cusp. We note that, when this analysis is performed on magnetic field lines which are not connected to the Earth (such as solar wind field lines beyond the dayside magnetopause, or O-lines inside magnetic islands), the "loss cone" population does not correspond to precipitating protons. For this reason, we refer to such protons as "a < 5°" or "a > 175°" rather than "precipitating protons", although in most of the analysed domain the plasma is on closed or semi-open field lines, and hence the calculated flux does correspond to precipitation. Since the fluxes are obtained by averaging the phase-space density within the loss cone (see Eq. (1)), and since field-aligned proton beams are typically significantly wider than 5°, the effect of choosing a fixed loss cone angle value everywhere rather than calculating its exact value (typically of the order of 2°, as mentioned in Sect. 2) can be neglected. Tick marks with 5°spacing have been added at the inner boundary to facilitate the interpretation of the figure;  we will refer to those angular values as "inner-boundary angles", noted Â IB , in the following.
The first feature that can be noted from Figure 8a is that the equatorward boundary of the cusps under southward IMF driving is sharp, with integral energy fluxes abruptly decreasing from above 10 7 keV cm À2 s À1 sr À1 to less than 5 Â 10 4 keV cm À2 s À1 sr À1 , at Â IB % +73°and Â IB % À71°in the northern and southern hemisphere, respectively. On the other hand, the poleward extent of the cusps is more difficult to determine, as integral energy flux values decrease more smoothly with increasing Â IB . High-flux (>10 8 keV cm À2 s À1 sr À1 ) regions have an extent in Â IB of 5°at both cusps, but moderate integral energy fluxes can be seen up to Â IB % ±85°.
Looking at the mean energy (Fig. 8b), one can see that the precipitating flux is harder near the equatorward edge of the cusps than near their poleward edge, with highest values of about 8 keV and 15 keV in the northern and southern hemispheres, respectively. This is consistent with the energy dispersion of precipitating protons observed in the cusp during southward IMF (e.g., Pitout et al., 2009). A possible explanation for the southern-hemisphere precipitation exhibiting a larger mean energy than in the northern hemisphere is that protons precipitating into the southern cusp originate from a large magnetic island moving towards the nightside near (X = 7 R E , Z = À6 R E ) (see also Supplementary Animation S1). The motion of the magnetic island towards the southern cusp might cause energisation of the loss-cone protons.
Away from the vicinity of the inner boundary, field-aligned proton integral energy fluxes have elevated values near the dayside magnetopause only, along the line where magnetic islands convect from the main reconnection site near the subsolar point towards the nightside. This is a further indication that cusp proton precipitation during southward IMF is associated with flux transfer events.
The analysis of the northward IMF restart file is presented in Figure 9, whose format differs from Figure 8 in that the fluxes and energies are shown separately for a < 5°protons propagating towards the northern cusp (Figs. 9a-9b) and for a > 175°protons propagating towards the southern cusp (Figs. 9c-9d). This enables visualising not only the high-latitude cusp spot precipitation originating from the same hemisphere's lobe reconnection, but also the interhemispheric coupling, with protons originating from the southern (northern) lobe reconnection site precipitating into the northern (southern) cusp.
The poleward extent of the cusp, measured along the inner boundary, is easily determined thanks to the abrupt drop in integral energy fluxes at Â IB % 84°and Â IB % À78°in the northern and southern hemisphere, respectively. This poleward boundary seems to be constrained by the lobe reconnection site location. On the other hand, the equatorward boundary of the proton precipitation is much less clear, as integral energy fluxes gradually decrease with decreasing Â IB values. Taking 10 7 keV cm À2 s À1 sr À1 as the reference value, integral energy fluxes remain elevated down to Â IB % 55°and Â IB % À60°i n the northern and southern hemisphere, respectively. In the northern hemisphere, the cusp appears split in two parts, due to the magnetic field topology, with a clear gap in precipitation within Â IB % 70-72°. This gap is not as prominent in the southern cusp.
Compared to the southward IMF situation, integral energy fluxes have maximum values of the same order of magnitude, slightly exceeding 10 8 keV cm À2 s À1 sr À1 in the southern cusp and around 5 Â 10 7 keV cm À2 s À1 sr À1 in the northern cusp. The mean precipitating energies lie within 8-12 keV in the equatorward portion of the cusp (on dayside field lines) and reach 15 keV in its poleward portion (on lobe field lines).

EMIC wave activity in the cusps
The velocity distributions shown in the inset panels of Figures 1, 2 and 7 are intrinsically unstable, which warrants an analysis of the wave activity in the cusps in the simulations. The ion/ion instabilities which may grow in a plasma are classified into four categories: ion/ion right-hand resonant (magnetosonic or fast MHD), ion/ion nonresonant (firehose), ion/ion left-hand resonant (Alfvén), and ion cyclotron anisotropy (electromagnetic ion cyclotron, EMIC) (Gary, 1991). The conditions for those instabilities to grow depend on the beam-core relative drift speed v 0 , the Alfvén speed v A , the beam density n b , the core population density n 0 , as well as the beam parallel and perpendicular temperatures T ||,b and T \,b .
Taking the example of the beam propagating towards the northern cusp in the distribution function seen at the dayside northern virtual spacecraft in Figure 2, we have v 0 = 1586 km/s, n b /n 0 = 0.02, and T \,b /T ||,b =10. At the virtual spacecraft location, v A = 1504 km/s. Under those conditions, the beam is stable for the firehose instability and at the limit of stability for the magnetosonic instability , weakly unstable for the ion/ion left-hand instability (Gary, 1985), and unstable for the ion cyclotron anisotropy instability (Gary & Schriver, 1987). One can therefore expect the ion cyclotron anisotropy instability to be the dominant one and to observe EMIC wave activity in the simulation.
Supplementary Animation S3 in the Supporting Information shows the y (out-of-plane) component of the magnetic field in the studied part of the simulation domain of the northward IMF run between t = 1100 s and t = 1938 s. Wave activity can be clearly seen in the cusp region, especially after t = 1100 s as lobe reconnection is happening in both hemispheres. Out-of-plane magnetic field disturbances originating at both lobe reconnection sites can also be seen propagating along the dayside magnetopause.
The left panels of Figure 10 show representative observations from a virtual spacecraft located in the northern cusp, at X = 2 R E and Z = 5 R E (dark-green square in Fig. 10d), in the northward IMF run. Fluctuations with a period of a few seconds are visible in the magnetic field strength (Fig. 10a) and the magnetic field components (Fig. 10b). Figure 10c displays the wavelet power spectrum of B y (Torrence & Compo, 1998), which confirms that there is significant wave power at a period of about 3 s during most of this interval. To determine the direction of the wave vector and the polarisation of the waves, we apply minimum variance analysis (MVA) (Sonnerup & Scheible, 1998) to the magnetic field components on 10 s sub-intervals. We consider that the determination of the minimum variance direction is reliable when the ratio of the intermediate to minimum eigenvalues is greater than 10 (Eastwood et al., 2005). The minimum variance direction provided by MVA gives us the direction of the wave vector with a 180°a mbiguity. At this virtual spacecraft, the waves propagate Earthwards in the simulation frame (see Supplementary Animation S3) with a velocity much larger than the plasma flow velocity in the cusps (not shown). This suggests that the waves also propagate Earthwards in the plasma frame. Therefore, we impose that the wave vector obtained from MVA should be directed in the ÀX direction. We find wave vectors that are within 30°from the mean magnetic field direction, consistent with the growth rate of EMIC waves maximising along the magnetic field vector (Gary & Schriver, 1987). Figures 10e and 10f show hodograms of the magnetic field components in the maximum (B l ) and intermediate (B m ) variance directions, that is, in the plane perpendicular to the wave vector, during two 10 s intervals marked with the black and purple bars in Figure 10a. The wave vector is in the out-of-plane direction, pointing towards the reader. Both hodograms reveal that the waves are left-hand polarised in the simulation frame. They retain the same polarisation in the plasma frame because the wave vector and the plasma flow are roughly parallel. This polarisation is again consistent with that of EMIC waves (Gary & Schriver, 1987).
A similar analysis has been done with the southward IMF run. Supplementary Animation S4 in the Supporting Information shows the y component of the magnetic field in the southward IMF run between t ¼ 1350 s and t ¼ 2150 s, wherein wave activity can be seen in the cusps throughout the animation. Although the wave characterisation proves more challenging in this run because of the frequent transit of FTEs in the cusp, the minimum variance analysis could be performed in several locations and confirmed the presence of EMIC waves identified thanks to their period being of the order of a few proton gyroperiods and their left-hand polarisation (not shown).

Discussion
The results presented above will be discussed while keeping in mind two limitations inherent to the simulation setup. First, the described precipitation fluxes are calculated at virtual spacecraft locations located relatively far from the ionosphere. This constraint comes from the fact that the inner boundary in the two Vlasiator simulations discussed here is located at 4.7 R E , implying that it is not possible to study the near-Earth plasma at altitudes where spacecraft monitoring particle precipitation (typically on low-Earth orbits under 1000 km altitude) are flown. One therefore has to assume that the proton populations studied at the virtual spacecraft locations would not undergo significant perturbations other than the increasing geomagnetic field magnitude while propagating towards the ionosphere. In the absence of wave-particle interactions and potential fields between the virtual spacecraft and the ionosphere, this assumption holds, following Liouville's theorem. Liang et al. (2013) estimated that typical potential drops in the auroral acceleration region (a few kV) have a negligible effect on precipitating ions whose energies are )1 keV. While this condition is mostly verified for the bulk of the precipitating populations under northward IMF, shown in Figures 4, 5 and 6, the lower-energy component of the precipitating fluxes and the southward-IMF proton precipitation shown in Figure 3 contain protons of only a few kiloelectronvolts, which may be somewhat affected by potential drops above the ionosphere. Since this process would take place beyond the inner boundary of the Vlasiator simulations, this effect cannot be taken into account here, and the lower-energy part of the precipitating spectra should therefore be considered with some caution, as such protons could undergo a deceleration by the upwarddirected potential drops (Coumans et al., 2004). In addition, we cannot rule out that the electromagnetic ion cyclotron (EMIC) waves seen Earthwards from the spacecraft along the same field line in Supplementary Animations S3 and S4 could scatter trapped protons into their loss cone, which has been suggested as one possible cusp precipitation mechanism for protons by Xiao et al. (2013). This means that the fluxes which we calculate at the virtual spacecraft locations are likely to be conservative low estimates of those which a particle detector on a spacecraft on a low-Earth orbit would observe. Nevertheless, as the integral energy fluxes and mean precipitating energies shown in Figures 8 and 9 do not appear to change significantly between the locations of the virtual spacecraft and the inner boundary, we can trust that the precipitating fluxes measured at the virtual spacecraft in the Vlasiator simulation domain are a relatively good representation of those that would be measured just above the ionosphere.
The presence of waves along the dayside magnetopause in the northward IMF simulation could explain the noted increase in precipitating proton energies between the dayside virtual spacecraft farthest from the cusp and the one closest to the cusp where given protons precipitate, as the associated precipitating beam gets slightly isotropised by interacting with those waves, thus adding beam protons into the corresponding loss cone. Testing this hypothesis would however require additional virtual spacecraft along the magnetopause, which is not possible with this run as full velocity distribution functions are not saved everywhere at each time step. This could be the topic of a future study.
The second limitation of these two simulations is related to the fact that they are 2D in ordinary space. Besides preventing the study of the MLT extent of the cusp, this requires using an adapted geomagnetic dipole. Indeed, using a point dipole in a meridional-plane simulation would violate the solenoidal condition; therefore, a line dipole centred at the origin and scaled to reproduce a realistic magnetopause standoff distance (Daldorff et al., 2014) is used in these two simulations. The main downside of this line dipole is that the field lines have a different shape compared to those obtained with a point dipole. While this only marginally affects the magnetic topology far enough from the Earth, this precludes any mapping of field lines to the Earth surface in terms of geomagnetic latitudes. Therefore, we cannot interpret the calculated fluxes and the precipitation maps shown in Figures 8 and 9 in terms of geomagnetic latitudes but rather described them in terms of inner-boundary angles instead. What however can be interpreted is the differences between the northward IMF and southward IMF simulations, even in terms of mapping, since the same line dipole has been used in both cases. For instance, the well-known fact that the cusp is located more polewards under northward IMF driving than under southward IMF driving (e.g., Sandholt et al., 1998;Frey, 2007;Němeček and Šafránková, 2008) is clearly seen in the comparison of our two simulations (see Figs. 8 and 9).
During southward IMF driving, the cusp precipitation can take place in the form of successive pulses (Mende et al., 1990;Lockwood et al., 1993) and often exhibits polewardmoving auroral forms, which are interpreted as being the result of pulsed reconnection at the dayside magnetopause Milan et al., 2000). The global simulation results presented above suggest that, rather than being due to pulsed reconnection, precipitation bursts seen in spacecraft data can be self-consistently tied to flux transfer events transiting in the cusp. Fuselier et al. (2007) reported an event observed with IMAGE FUV driven by predominantly southward IMF, which exhibited clear fluctuations in keV-proton precipitation over time scales of $5 min. These times scales are larger than those of the fluctuations associated with the precipitation bursts in our southward IMF simulation, which are closer to $1 min (Fig. 3). A likely explanation for the shorter time scales seen in Vlasiator is the 2D setup of the simulation, which imposes that each IMF field line must reconnect at the dayside magnetopause. This leads to the formation of more FTEs which in turn translates into more frequent proton precipitation bursts at the cusps.
Recently, Hoilijoki et al. (2019) presented a statistical study of the flux transfer events which are formed in our southward IMF run (Run A in Hoilijoki et al., 2019). In particular, they studied the location of the flux transfer events along the dayside magnetopause as well as their sizes as a function of time (see their Fig. 5a). One can see that most of the major flux transfer events in either the northern or southern hemisphere tracked in their study are associated with a precipitation burst in our study (Figs. 3a and 3c) slightly after they drift out of their tracking area. This gives a further confirmation that the proton precipitation bursts in our southward IMF simulation are driven by incoming magnetic islands formed at the dayside magnetopause by magnetic reconnection.
The formation of these FTEs, their migration towards the northern or southern cusp and the resulting precipitation bursts exhibit hemispheric asymmetries despite the simulation setup being an idealised situation with purely southward IMF and no dipole tilt. Such asymmetries likely arise from numerical effects (such as small rounding errors) introducing very low levels of small-scale fluctuations already early in the simulation. This breaks the ideal symmetric topology and results in large-scale differences between the northern and southern cusps. The fact that small-scale fluctuations can evolve into prominent inter-hemispheric asymmetries reveals how complex the study of those asymmetries is and how difficult it may be to understand what the main sources of differences between the two hemispheres are, besides seasonal dipole tilt effects (e.g., Newell & Meng, 1988;Shi et al., 2012), intrinsic geomagnetic field asymmetries (e.g., Cousins & Shepherd, 2012), and IMF orientation (e.g., Newell et al., 1989). Nevertheless, when considering time scales larger than that of individual FTE transit, the hemispheric asymmetries tend to be averaged out (Hoilijoki et al., 2019).
Once the FTEs reach the poleward region of the cusp, they start reconnecting with the magnetospheric lobe, which magnetically connects the hot plasma they contain to the inner boundary (and ultimately to the ionosphere) while leading to their disintegration. This process, which can be observed throughout Supplementary Animation S1, is consistent with the findings by Omidi & Sibeck (2007) who showed that FTEs get entirely disintegrated "within or tailward of the cusp" and inject their plasma into the cusp, which can be accelerated by reconnection jets. Omidi & Sibeck (2007) found that such plasma injection from FTEs into the cusp can have signatures in ground-based auroral observations consistent with poleward-moving auroral forms, since the reconnection site between a FTE and the lobe moves polewards as the FTE is advected towards the nightside. While we cannot check whether we would obtain similar signatures because of the limited number of possible virtual spacecraft to study precipitation in the southward IMF run, the examination of Supplementary Animation S1 suggests that there could be a qualitative agreement between our results and those from Omidi & Sibeck (2007) in that respect.
On the other hand, in the northward IMF simulation, the high-latitude cusp spots of proton precipitation exhibit much more stable precipitating fluxes compared to the southward IMF case (see Fig. 6). This suggests that lobe reconnection during northward IMF is more continuous than the subsolar magnetopause reconnection during southward IMF which generates more flux transfer events. Besides, at the high-latitude cusp spots in the northward IMF simulation, the mean precipitating energy increases with increasing inner-boundary angles (see Figs. 9b and 9d, inner-boundary angles near ±80-85°). This is consistent with the reverse-dispersion signature reported in observations of cusp precipitation during northward IMF driving (e.g., Øieroset et al., 1997).
The particular case of cusp auroral precipitation during dual lobe reconnection events is less well-known, as dual lobe reconnection requires a nearly purely northward IMF, with a clock angle within ±10° (Imber et al., 2006). While Østgaard et al. (2005) have produced the first simultaneous observations of cusp aurora spots in both hemispheres during northward IMF, it is likely that the spot in a given hemisphere originated from the same hemisphere's lobe reconnection site, hence analogous to the poleward precipitation areas seen in Vlasiator (calculated at the so-called cusp spot virtual spacecraft). A possible ionospheric signature of particle precipitation in the cusp during dual lobe reconnection and originating from the lobe reconnection site in the opposite hemisphere was reported by Imber et al. (2006), who identified a brightening of the northern-hemisphere dayside aurora at a time when dual lobe reconnection was initiated. They interpreted this brightening as a sign that particles coming from the southern lobe reconnection site were precipitating into the northern cusp, which corresponds to what is seen in the northward IMF Vlasiator simulation.
Besides, Pitout et al. (2012) reported Cluster observations of overlapping populations of precipitating ions above the cusp during an event characterised by strongly northward IMF. They suggested that the second ion population, precipitating into the equatorward part of the cusp and with energies of the order of 10 keV, might be originating from the reconnection at the opposite hemisphere lobe. This seems in good agreement with the precipitation observed at the dayside virtual spacecraft in the northward IMF simulation of our study, which indeed map to the equatorward part of the cusp (see Figs. 4,5,and 8). Chang et al. (2004) also reported two types of precipitation in the vicinity of the cusp during northward IMF in IMAGE FUV and DMSP data, and interpreted the equatorward precipitation as originating from the LLBL, on field lines "that had been merged for a very long time or were closed". It is possible that this LLBL population is analogous to that observed at the dayside spacecraft in our northward IMF simulation. Indeed, in a study of ion distribution functions in the dayside magnetosphere during northward IMF, Fuselier et al. (2014) show a three-ion-population distribution in the LLBL which presents very similar features as those seen at the three dayside virtual spacecraft in our northward IMF simulation (see Fig. 2 and Supplementary Animation S2).
Our results are also consistent with magnetohydrodynamics (MHD) coupled with test-particle simulations carried out by Connor et al. (2015), where reversed-dispersion signatures are seen in ion cusp precipitation in the northward IMF simulation, while nonreversed dispersion signatures (i.e., decreasing energies with increasing latitude) are seen in the southward IMF simulation. In their northward IMF simulation, precipitating ion fluxes in the cusp peak at energies of 10-30 keV, whereas in their southward IMF simulation they peak at energies below 10 keV. This is also in good agreement with what is seen in our kinetic simulations. Finally, our results with the southward IMF run can be compared to those obtained by Tan et al. (2012) with an MHD model, wherein a nonreversed dispersion in the cusp ion precipitation is also seen.
In conclusion, the results shown in this paper therefore not only reproduce results obtained using MHD simulations and hybrid-PIC codes such as in Omidi & Sibeck (2007), but they also reproduce for the first time observations of velocity distributions with several field-aligned beams propagating in opposite directions during dual-lobe reconnection, and they give support to the claim made by Xiao et al. (2013) that EMIC waves can play a role in cusp proton precipitation. Further modelling studies of cusp precipitation using self-consistent hybrid kinetic simulations are foreseen once 3D Vlasiator runs are available. In future runs, particular attention will be given to saving velocity distribution functions in a larger number of cells near the inner boundary to enable more detailed studies of the latitudinal and magnetic local time variations of the calculated precipitating fluxes.

Summary
This paper presented the first self-consistent ion kinetic simulation results on cusp proton precipitation during purely southward and purely northward IMF driving, using the Vlasiator code. We used two global 2D hybrid-Vlasov simulations of the near-Earth space and calculated precipitating proton differential number fluxes, integral energy fluxes and mean energy at several locations coined "virtual spacecraft", using the same method as was applied to nightside precipitation in Grandin et al. (2019).
The first simulation was driven by a solar wind with purely southward IMF. The proton precipitation was studied at two virtual spacecraft located above the northern and southern cusp. The precipitating fluxes exhibited bursts of protons with decreasing energy as a function of time, occurring irregularly at a rate of about one per minute on average. Precipitating fluxes at the northern cusp were independent from those at the southern cusp, yet exhibited similar orders of magnitude in terms of mean energy, integral energy and differential number fluxes. We found an unambiguous correlation between the occurrence of precipitation bursts and the transit of flux transfer events in the vicinity of the virtual spacecraft.
The second simulation was driven by the same solar wind as the first simulation, except for the IMF orientation which was purely northward. In this simulation, dual lobe reconnection was taking place. Proton precipitation was studied at three virtual spacecraft in the dayside magnetosphere, mapping to the equatorward part of the cusps. We found signatures of dual lobe reconnection in the velocity distribution functions of these three spacecraft, which exhibited superpositions of a core population and several field-aligned (parallel and antiparallel) beams during part of the simulation. Correspondingly, precipitating fluxes were seen at these three virtual spacecraft, with protons originating from the southern (northern) lobe reconnection site precipitating into the northern (southern) cusp. We suggest that such precipitation originating from the opposite hemisphere lobe reconnection, albeit infrequent, could be observed in the equatorward part of the cusps during dual lobe reconnection.
Two additional virtual spacecraft monitored the highlatitude precipitation spot, where precipitating protons originate from the same hemisphere's lobe reconnection site. We found that the precipitating fluxes and energies at those high-latitude cusp spots are relatively steady and are not associated with flux transfer events.
The comparison of the northward IMF and southward IMF cusps showed that in the former, the precipitation takes place more polewards and exhibits reverse dispersion, whereas in the latter the precipitation takes place more equatorwards and exhibits nonreversed dispersion. At the locations where they reach their maximum values, the integral energy flux and mean precipitating energy have values of the same order of magnitude in both simulations:~10 8 keV cm À2 s À1 sr À1 and 8-15 keV, respectively.
In both runs, unambiguous signatures of EMIC wave activity are present in the cusps, underlining one of the assets of using a kinetic model in such a study. Since EMIC waves can scatter protons into the loss cone, the precipitating fluxes calculated at the virtual spacecraft in the simulations are likely to be conservative low estimates of the fluxes which would reach the ionosphere.
This work will contribute to improving the understanding of auroral precipitation thanks to global kinetic magnetospheric models, whose continued development is a key endeavour for space weather understanding (Robinson et al., 2019;Heelis & Maute, 2020). While the current results are obtained in 2D and at a few locations within the cusp region, future Vlasiator simulations will aim at studying particle precipitation in 3D and within each cell of the simulation domain, thus allowing a more detailed description of the spatial variability and temporal dynamics of the auroral oval during geomagnetic events.
Acknowledgements. We acknowledge the European Research Council for Starting grant 200141-QuESpace, with which Vlasiator  was developed, and Consolidator grant 682068-PRESTISSIMO awarded to further develop Vlasiator and use it for scientific investigations. The Finnish Centre of Excellence in Research of Sustainable Space (FORESAIL), funded through the Academy of Finland grant number 312351, supports Vlasiator development and science as well. The work of L.T. is supported by the Academy of Finland (grant number 322544). PRACE (http://www. prace-ri.eu) is acknowledged for granting us Tier-0 computing time in HLRS Stuttgart, where Vlasiator was run in the HazelHen machine under project number 2014112573. The CSC -IT Center for Science in Finland is acknowledged for the Sisu and Taito supercomputer usage, which led to the results presented here.
The runs described here take several terabytes of disk space and are 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 (https://www.helsinki.fi/en/researchgroups/vlasiator/rules-ofthe-road). The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Supplementary Materials
The supplementary material of this article is available at http://www.swsc-journal.org/10.1051/swsc/2020053/olm The Supporting Information consists of four animations numbered S1-S4, whose descriptions are given below.
Animation S1. Animation of the proton temperature in a zoomed-in box centred on the dayside magnetosphere in the southward IMF simulation. The two locations in the northern and southern cusps where the precipitating fluxes can be derived at each time step are shown with black circles, and the inset axes show cuts of the corresponding VDFs in the (v x , v z ) frame, with the loss cone indicated in magenta (northern cusp VDF) or brown (southern cusp VDF). The black lines represent magnetic field lines. The velocity grid spacing is 1000 km/s. The animation starts at t = 1350 s and lasts until the end of the simulation at t = 2150 s, with one frame every 0.5 s. (IMF: interplanetary magnetic field; VDF: velocity distribution function).
Animation S2. Animation of the proton temperature in a zoomed-in box centred on the dayside magnetosphere in the northward IMF simulation. The five locations shown with black circles indicate the cells where the precipitating fluxes are calculated, and the inset axes show cuts of the corresponding VDFs in the (v x , v z ) frame, with the loss cone indicated in magenta and brown for precipitation in the northern and southern cusps, respectively. The black lines represent magnetic field lines. The velocity grid spacing is 1000 km/s. The animation starts at t = 1100 s and lasts until the end of the simulation at t = 1938 s, with one frame every 0.5 s. (IMF: interplanetary magnetic field; VDF: velocity distribution function).
Animation S3. Animation of the y (dawn/dusk) component of the magnetic field in a zoomed-in box centred on the dayside magnetosphere in the northward IMF simulation. The five locations shown with black circles indicate the cells where the precipitating fluxes are calculated, and the inset axes show cuts of the corresponding VDFs in the (v x , v z ) frame, with the loss cone indicated in magenta and brown for precipitation in the northern and southern cusps, respectively. The black lines represent magnetic field lines. The velocity grid spacing is 1000 km/s. The animation starts at t = 1100 s and lasts until the end of the simulation at t = 1938 s, with one frame every 0.5 s. (IMF: interplanetary magnetic field; VDF: velocity distribution function).
Animation S4. Animation of the y (dawn/dusk) component of the magnetic field in a zoomed-in box centred on the dayside magnetosphere in the southward IMF simulation. The two locations in the northern and southern cusps where the precipitating fluxes can be derived at each time step are shown with black circles, and the inset axes show cuts of the corresponding VDFs in the (v x , v z ) frame, with the loss cone indicated in magenta (northern cusp VDF) or brown (southern cusp VDF). The black lines represent magnetic field lines. The velocity grid spacing is 1000 km/s. The animation starts at t = 1350 s and lasts until the end of the simulation at t = 2150 s, with one frame every 0.5 s. (IMF: interplanetary magnetic field; VDF: velocity distribution function).