Hybrid-Vlasov modelling of nightside auroral proton precipitation during southward interplanetary magnetic field conditions

Particle precipitation plays a key role in the coupling of the terrestrial magnetosphere and ionosphere by modifying the upper atmospheric conductivity and chemistry, driving field-aligned currents, and producing aurora. Yet, quantitative observations of precipitating fluxes are limited, since ground-based instruments can only provide indirect measurements of precipitation while particle telescopes onboard spacecraft merely enable point-like in-situ observations with inherently coarse time resolution above a given location. Further, orbit time scales generally prevent the analysis of whole events. On the other 5 hand, global magnetospheric simulations can provide estimations of particle precipitation with a global view and higher time resolution. We present the first results of auroral (∼1–30 keV) proton precipitation estimation using the Vlasiator global hybridVlasov model in a noon-midnight meridional plane simulation driven by steady solar wind with southward interplanetary magnetic field. We first calculate the bounce loss cone angle value at selected locations in the simulated nightside magnetosphere. Then, using the velocity distribution function representation of the proton population at those selected points, we study the 10 population inside the loss cone. This enables the estimation of differential precipitating number fluxes as would be measured by a particle detector onboard a low-Earth-orbiting spacecraft. The obtained differential flux values are in agreement with a well-established empirical model in the midnight sector, as are also the integral energy flux and mean precipitating energy. We discuss the time evolution of the precipitation parameters derived in this manner in the global context of nightside magnetospheric activity in this simulation, and we find in particular that precipitation bursts of< 1min duration can be self-consistently 15 and unambiguously associated with dipolarising flux bundles generated by tail reconnection. We also find that the transition region seems to partly regulate the transmission of precipitating protons to the inner magnetosphere, suggesting it has an active role in regulating ionospheric precipitation.


Introduction
The terrestrial atmosphere and ionosphere are known to be affected by the precipitation of particles coming from the magnetosphere and the solar wind.Precipitating protons in the keV energy range produce diffuse auroral emission, principally in the based optical instruments can be used in conjunction with satellite observations (e.g., Donovan et al., 2008).Besides, since energetic proton precipitation produces ionisation enhancements in the E region of the ionosphere, ionospheric radars such as incoherent scatter radars can provide indirect observations of proton precipitation (e.g., Lyons et al., 2010).
In numerical studies, two aspects of proton precipitation can be distinguished.First, empirical models describing precipitation fluxes and locations have been developed.One example of these is the Hardy model (Hardy et al., 1989(Hardy et al., , 1991) ) which provides the average precipitation patterns as a function of the Kp index obtained by compiling two years of DMSP spacecraft data.This statistical model describes proton precipitation in a spatial grid containing 30 bins in corrected geomagnetic latitude above 50 • and 48 bins in MLT, at seven levels of geomagnetic activity from Kp = 0 to Kp ≥ 6−.The model gives the proton precipitation differential number flux in the 20 DMSP energy channels (from 30 eV to 30 keV with regular spacing in logarithmic scale) as well as the integral number flux, the integral energy flux and the mean precipitating energy.The Hardy model has been used as an input in a number of ionospheric models, e.g., Transcar and IPIM (Blelly et al., 2005;Marchaudon and Blelly, 2015).More recently, the OVATION Prime model (Newell et al., 2014), derived from a combination of DMSP particle data and global ultraviolet imager (GUVI) data from the TIMED spacecraft, has been developed and includes proton precipitation.However, due to the lack of DMSP satellites orbiting in the postmidnight sector, the magnetic local time (MLT) sectors comprised between 00:15 and 03:30 MLT are described using a linear interpolation.
Second, models based on first principles focus essentially on ionospheric effects of precipitating ions, such as the linear transport model described in Basu et al. (2001).In such approaches, precipitating fluxes are taken as model inputs, and the transport code provides ion production rates and auroral emission intensity profiles in the ionosphere.From a magnetospheric point of view, while simulations of electron precipitation have been achieved with various models (e.g., Palmroth et al., 2006;Raeder et al., 2008), attempts to model proton precipitation are scarce.A few magnetohydrodynamics (MHD) models have addressed the issue of relating particle precipitation characteristics to the magnetotail dynamics (Gilson et al., 2012;Ge et al., 2012, see below); however, currently no global kinetic simulations of the near-Earth environment have undertaken this task.
One obvious aspect through which proton precipitation and magnetotail dynamics can be tied is the relationship between auroral streamers and bursty bulk flows (BBFs).Bursty bulk flows are fast plasma flows propagating in the magnetotail along the Sun-Earth direction, either earthwards or away from the Earth (Angelopoulos et al., 1992(Angelopoulos et al., , 1994)).They are also characterised by an enhanced B z (northward component of magnetic field) and reduced plasma pressure (e.g., Ohtani et al., 2004).A given BBF may embed several dipolarising flux bundles (DFBs), which are themselves coherent magnetic structures exhibiting a B z increase by up to a few tens of nanoteslas (Liu et al., 2013) and may play a role in particle energisation in the magnetotail (Runov et al., 2017).Auroral streamers are thin auroral structures oriented in the north-south direction (e.g., Nishimura et al., 2011).They have been associated with BBFs based on observational evidence combining spacecraft data with ground-based auroral imager, magnetometer and coherent scatter radar data (e.g., Fairfield et al., 1999;Nakamura et al., 2001a, b;Amm and Kauristie, 2002;Sergeev et al., 2004;Gallardo-Lacourt et al., 2014).On modelling aspects, Ge et al. (2012) reproduced proton precipitation enhancements associated with approaching BBFs with a global MHD model and found good agreement with spacecraft and ground-based auroral observations.In this paper, we present an overview of nightside proton precipitation in a global magnetospheric hybrid-Vlasov simulation under southward interplanetary magnetic field (IMF) conditions, using the Vlasiator model (von Alfthan et al., 2014;Palmroth et al., 2018).Contrary to previous work in which the main focus was on the ionosphere, we here look at precipitation from the perspective of the magnetosphere and unambiguously tie the precipitation characteristics to the self-consistent simulation results.In particular, following the study by Juusola et al. (2018a) investigating flow bursts inside a BBF in the same Vlasiator run, we examine how proton precipitation can be related to dipolarising flux bundles.Section 2 describes the Vlasiator model as well as the methods developed to estimate proton precipitation from its outputs.Section 3 presents the results of the study, i.e., the features of nightside proton precipitation in the simulation under southward IMF.Section 4 discusses those results, and sect.5 summarises the main conclusions.

Vlasiator
This study relies on numerical simulations using Vlasiator, a global hybrid-Vlasov model of the near-Earth plasma environment (von Alfthan et al., 2014;Palmroth et al., 2018).In the hybrid-Vlasov approach, protons are described as velocity distribution functions (VDFs) in a grid of the phase space (ordinary space and velocity space), whereas electrons are treated as a massless, charge-neutralising fluid.The time evolution of the VDFs is obtained by solving the Vlasov equation, and closure of the system is achieved with Ohm's law including the Hall term (Palmroth et al., 2018).
The simulation run used in this study is a 2D-3V (two-dimensional in ordinary space, three-dimensional in velocity space) simulation in the noon-midnight meridional plane of the magnetosphere, i.e., XZ in the geocentric solar magnetospheric (GSM) coordinate system.The origin of the simulation domain is set in the Earth's centre.The simulation box in ordinary space spans from X = −94R E (Earth radii, 1R E = 6371 km) on the nightside to X = +47R E on the dayside, and is comprised between boundaries at Z = ±57R E in the north-south direction, with a grid resolution of 300 km.The inner boundary of the magnetospheric domain lies at ∼ 4.7R E from the center of the Earth.In each cell of the ordinary space, a velocity grid extends between ±4020 km s −1 with a resolution of 30 km s −1 in V X , V Y , and V Z .A sparsity threshold is implemented in Vlasiator, below which phase-space density values are discarded to limit the computational load of the simulation.In the run used in this study, the sparsity threshold is set to 10 −15 m −6 s 3 .The near-Earth environment is driven by incoming solar wind from the +X wall of the simulation domain.The solar wind population consists of protons characterised by Maxwellian VDFs with a number density of 1 cm −3 , a temperature of 500 kK, a bulk velocity of −750 km s −1 along the X axis, a dynamic pressure of about 0.9 nPa, and carrying a purely southward interplanetary magnetic field (IMF) with a magnitude of 5 nT.Other boundary conditions are as follows: von Neumann condition for the −X and ±Z walls, periodic conditions in the out-of-plane direction (±Y ), whereas the inner boundary is a perfectly conducting cylinder with a static Maxwellian VDF for protons.Finally, the geomagnetic field is given by a 2D line dipole along the Z axis, centred at the origin, and scaled to obtain a realistic magnetopause standoff distance (Daldorff et al., 2014).Figure 1 gives an overview of the simulated area at a timestep (t = 1800 s) when nightside reconnection is taking place.The plotted parameter is the plasma temperature, and the black lines correspond to magnetic field lines.The incoming solar wind at 500 kK with southward IMF can be seen on the right-hand side of the figure, and the major geospace regions (bow shock, magnetosheath, dayside magnetopause, polar cusps, magnetotail lobes, plasma sheet) can be identified.A plasmoid is forming in the nightside magnetosphere (Palmroth et al., 2017), whose signature is visible between X = −30R E and X = −20R E .This simulation run was previously used in several studies focusing on various phenomena and regions in the near-Earth space, such as dayside magnetopause reconnection (Hoilijoki et al., 2017), ion acceleration in the magnetosheath (Jarvinen et al., 2018), magnetotail reconnection (Palmroth et al., 2017;Juusola et al., 2018a), and magnetotail current sheet flapping (Juusola et al., 2018b).In each of these studies, the presented results showed good correspondence with earlier findings.
Here, this run will be used to characterise nightside auroral proton precipitation in the simulation by making use of the VDF

Precipitation spectra calculation
Figure 2 illustrates the definitions of angles and vectors used in this section.The configuration is shown for particles precipitating in the southern hemisphere to be consistent with the results shown in the rest of the paper.For particles precipitating in the northern hemisphere, the configuration is symmetrical.
At a given three-dimensional location r of the ordinary space, the differential intensity J (more generally called directional differential flux) of protons at energy E for the direction along unit vector u is related to the velocity distribution function f by (see 6.4.1. in Baumjohann and Treumann, 1997) where m p is the proton mass and v = 2E/m p is the velocity magnitude.J , when expressed in SI units, is given in particles Particle detectors onboard spacecraft estimate J by taking the average of the incoming particle flux in a given energy channel over the detector's surface S det and viewing angle Ω det , which can be written where u is the unit vector of a given velocity direction inside Ω det , n det is the unit vector normal to the detector's surface, and r sc is the location of the spacecraft, for instance near the top of the ionosphere.Strictly speaking, Ĵ is a differential intensity measured within a fixed solid angle and detector area.In practice, telescopes onboard spacecraft generally collect particles within a cone with an opening angle θ det of the order of 15 • (e.g., Rodger et al., 2010).Assuming that the directional or, by defining µ = cos θ, with µ det = cos θ det .In the symmetrical configuration where the spacecraft observes precipitating particles in the northern hemisphere, n det is aligned with −b sc , but angle definitions remain the same.
As a first-order approximation, protons remain attached to a given magnetic flux tube; one can hence trace particles backwards in time out to the magnetosphere, following magnetic field lines.Assuming smooth quasilinear propagation, for a given particle of pitch angle θ, the quantity sin 2 θ/B is conserved along the trajectory.Let r 0 be a location in the magnetosphere mapping to the spacecraft in terms of magnetic field topology.For particles travelling between r 0 and r sc , we have where B 0 and θ 0 correspond to the magnetic field magnitude and the particle pitch angle at r 0 , respectively, and B sc and θ correspond to the same parameters at r sc .This can also be written with µ 0 = cos θ 0 .By differentiating eq. ( 6), one gets According to Liouville's theorem, f (r, v, µ, ϕ) is conserved along the trajectories of the system, i.e., in the absence of potential fields Using eq. ( 8), and then equations ( 6) and ( 7), yields from eq. ( 4) where b 0 is the unit vector with the direction of the local magnetic field at r 0 , and In other words, the directional differential particle flux observed by the telescope onboard the spacecraft can be estimated by averaging the relevant subset of the particle VDF at a chosen conjugate location in the magnetosphere.The angular boundaries of the averaged phase space domain are obtained by scaling the viewing angle of the detector to θ det,0 = arccos µ det,0 based on the ratio of magnetic field magnitudes at the spacecraft and at the conjugate location in the magnetosphere where the VDF is analysed (cf.eq. ( 10)).
For particles located at r 0 to reach a spacecraft at the top of the ionosphere at ∼800 km altitude, their pitch angles must be within the bounce loss cone, whose opening angle is determined by In the nightside equatorial plane, θ lc takes values of a few degrees at most, and is smaller than 1 • beyond 10 R E (see, e.g., Sergeev and Tsyganenko, 1982).In Vlasiator, knowing the VDF at a given location in the magnetosphere, one can calculate the bounce loss cone angle value and, at a given velocity magnitude v (i.e., at a given energy), average the phase-space density inside the loss cone to evaluate Ĵ (E, b 0 , r 0 ).It should be noted that, with this approach, the obtained directional differential precipitating flux corresponds to the one which would be measured by a telescope onboard a spacecraft at the topside ionosphere with a viewing angle Ω det = 2π sr (cf.eq. ( 2)).
In practice, the methodology to estimate differential precipitating proton fluxes at r 0 is as follows.First, the loss cone angle θ lc is calculated.This is achieved by following the magnetic field line from r 0 to the inner boundary of the simulation domain, at about 4.7 R E , and then extrapolate the magnetic field to the top of the ionosphere using the line dipole approximation.The ratio of the magnetic field magnitudes at r 0 and at the mapped region in the topside ionosphere enables the calculation of the bounce loss cone angle using eq.( 11).Second, we estimate the differential precipitating flux by calculating where f (r 0 , v, θ, ϕ) θ<θ0 is the average value of the phase-space density at speed v inside the bounce loss cone.This approximation of eq. ( 9) is reasonable since θ 0 is very small and hence we have µ 0 = cos θ lc ≈ 1 inside the integral.

Examples of differential precipitating proton fluxes with Vlasiator
In the simulation used in this study, full velocity distributions of protons are saved every 50 simulation cells in the X and Z direction due to limitations in the file sizes.Therefore, the estimation of the differential precipitating proton flux based on eq. ( 12) can be applied only in a few selected locations in the nightside magnetosphere.(approximately aligned with the Y direction, i.e., into the simulation plane).The grid spacing is 1000 km/s.The velocity distribution at S 1 consists of a core population centred near the origin and a crescent-shaped beam at about 1000 km/s.It corresponds to the superposition of a cold plasma with a low bulk velocity in the magnetic field direction with a more energetic population travelling earthwards.The velocity distribution at S 2 is nearly Maxwellian with a broad loss cone in the magnetic field direction (this loss cone will be discussed below).The width of the distribution, centred near the origin, indicates that it corresponds to a hotter plasma (in comparison with the one at S 1 ) nearly at rest in the local magnetic frame, with protons having velocities up to about 2000 km/s.This contrast between the velocity distributions at S 1 and S 2 is reflected in the plasma temperature, which is 15 MK at S 1 and 32 MK at S 2 .In panels (b) and (c), the bounce loss cone is indicated with pink lines.
As can be seen in the two examples, the loss cone is not empty, and hence a directional differential flux of precipitating protons can be derived using eq.( 12).We note that the empty part of the phase space in panel (c) is due to the presence of the inner boundary.Indeed, the inner boundary absorbs incoming protons, and particles bouncing back from a mirror point earthwards from the inner boundary in the southern hemisphere are therefore absent at virtual spacecraft S 2 in the simulation.This feature does not have consequences in our study, since the relevant part of the velocity distribution function is the one located inside the bounce loss cone.

Nightside proton precipitation
Figure 4 shows four examples of velocity distributions at virtual spacecraft S 1 and S 2 during the simulation run (panels a-d) and the associated precipitating proton differential fluxes obtained using the method described in section 2.2 (panels e-h).
The velocity distributions are slices in the plane defined by the magnetic field direction, v B , and the electric field direction, eV −1 ).Observations from the DMSP 12 spacecraft reported in Lummerzheim et al. (2001) indicate mean precipitating energy of the order of 10 keV and flux values around 10 3 cm −2 s −1 sr −1 eV −1 .These values were however obtained in the evening MLT sector, where proton precipitation exhibits statistically harder energy spectra than in the midnight sector (Galand et al., 2001).
An overview of the dynamics of the nightside magnetosphere from t = 1000 s until the end of the simulation at t = 2150 s is provided with supplementary animation S1. Figure 5 shows the time evolution of the proton precipitation at virtual spacecraft S 1 and S 2 during this time interval.Panels (a) and (b) show the precipitating proton differential flux (colour scale) as well as the mean precipitating energy (black line) as a function of time, while panels (c) and (d) show the integral energy flux.
It can be seen that, at virtual spacecraft S 1 , broad-energy proton precipitation above 1 keV starts around t = 1360 s, as the magnetic field line observed by S 1 is becoming slightly stretched.The energetic tail of the proton precipitation spectrum reaches up to 20 keV and the mean precipitating energy fluctuates between 2 and 5 keV until t ≈ 1750 s, i.e., about 90 s after the global magnetotail reconfiguration is initiated by the dominant X-line near X = −13R E .Around t = 1750 s, the precipitation becomes narrower in energy and, after t = 1850 s, the mean precipitating energy increases from 4 keV to 20 keV, while fluxes decrease by almost two orders of magnitude.This corresponds to the times during which the nightside reconnection speeds up and S 1 observes magnetic field lines becoming more dipolar.Between t = 1900 s and t = 1925 s, S 1 does not observe any precipitation, as the virtual spacecraft is momentarily observing lobe-type plasma and the local loss cone is empty, but precipitation resumes with high energy (10-20 keV) and relatively low differential flux values (∼ 10 2 cm −2 s −1 sr −1 eV −1 ) after t = 1925 s, as S 1 is magnetically connected to the transition region.A low-energy (< 1 keV) precipitating population appears between t = 1990 s and t = 2040 s, presumably associated with heating of the core plasma.At the end of the simulation, the precipitating mean energy decreases while differential flux values are enhanced to ∼ 10 3 cm −2 s −1 sr −1 eV −1 .The integral energy flux remains within 1-6 × 10 7 keV cm −2 s −1 sr −1 during the broad-energy precipitation phase, but then drops between 10 6 and 10 7 keV cm −2 s −1 sr −1 during the phase when field lines passing by S 1 become dipolar (t =1800-1900 s).When precipitation observations from S 1 resume (t ≈ 1930 s), the integral energy flux remains close to 10 7 keV cm −2 s −1 sr −1 before briefly increasing to reach 10 8 keV cm −2 s −1 sr −1 and abruptly decreasing around 10 6 keV cm −2 s −1 sr −1 until the end of the simulation.Those integral energy flux values are in agreement with the Hardy model, according to which the nightside maximum total energy flux ranges between about 4 × 10 7 and 2 × 10 8 keV cm −2 s −1 sr −1 , depending on the Kp index value (Hardy et al., 1989).
At virtual spacecraft S 2 (panels b and d), proton precipitation above 1 keV appears around t = 1760 s in the simulation, corresponding roughly to the time when the broad-energy precipitation at S 1 becomes narrower in energy.For about one minute, broad-energy precipitation with up to 30 keV protons and with a mean energy of 6-7 keV is observed, leading to integral flux values reaching 2 × 10 8 keV cm −2 s −1 sr −1 .The precipitating energy spectrum then becomes narrower and centred around 2 keV (during t =1850-1900 s) and gradually broadens again until the end of the simulation, when the mean precipitating energy decreases below 1 keV.
While the 2D simulation setup with a scaled line dipole does not enable a direct mapping of the virtual spacecraft locations in terms of geomagnetic latitudes or even L values, it appears clearly that the properties and the dynamics of proton precipitation are very different at S 1 and S 2 .From the geometry, S 1 maps to higher geomagnetic latitudes than S 2 , and the analysis of Fig. 5 suggests that S 1 represents well geomagnetic latitudes usually mapping to the auroral oval, while S 2 would correspond to slightly lower latitudes where the proton auroral arc can be observed around the onset of ionospheric substorms, after drifting equatorwards during the growth phase (Liu et al., 2007).Indeed, precipitation is observed at S 1 before the global magnetotail reconfiguration due to reconnection is initiated (around t = 1660 s; Palmroth et al. ( 2017)), which suggests that S 1 maps to geomagnetic latitudes where aurora is visible during the growth phase of substorms.On the other hand, proton precipitation is only observed at S 2 around the time when the first particles accelerated earthwards following the global magnetotail reconfiguration reach the inner magnetosphere.

Flow bursts and precipitation
In what follows, we will focus on the proton precipitation associated with dipolarising flux bundles by considering virtual spacecraft S 1 between t = 1920 s and t = 2080 s.During this time span, S 1 is on magnetic field lines mapping to the transition region near X = −9R E , as can be seen in supplementary animation S2.Juusola et al. (2018a) showed that, during this time  bulk properties of the plasma between S 1 and the current sheet through the transition region.As can be seen in the figure, they are not located on a same field line, but rather along the region of positive V x , which corresponds to the path followed by precipitating protons as field lines are being convected earthwards.Full velocity distributions are not available at those points due to limitations in the file sizes.The spacecraft indicated with an orange plus sign, which is the furthest in the current sheet, will be called S 0 .
Figure 7 shows the precipitation flux observed at S 1 between t = 1920 s and t = 2080 s, corresponding to the time interval marked with a green line in panels (a) and (c) of Fig. 5. Panel (a) of Fig. 7 shows the differential number flux of precipitating protons (colour scale) and the mean precipitating energy (black line).One can visually identify successive bursts of precipitating protons at energies ranging between about 5 keV and 50 keV.The energy dispersion of precipitating protons is visible, as within a given precipitation burst the highest energies are observed first and the lowest energies are observed last.This is also visible in the time variations of the mean precipitating energy.Between t = 1985 s and t = 2035 s, a low-energy (< 1 keV) population of precipitating protons is observed in addition to the ∼10 keV precipitation, but with fluxes lower than the main precipitating population around 10 keV by almost two orders of magnitude.The differential flux of the main precipitating proton population at 5-50 keV has values around 10 2 cm −2 s −1 sr −1 eV −1 , which is relatively low compared to flux values at other times in the simulation that can reach 10 4 cm −2 s −1 sr −1 eV −1 (cf.Fig. 5a).Panel (b) shows the integral energy flux (blue line) alongside the component of the plasma bulk velocity at S 1 parallel to the magnetic field (red line), V .Signatures of the precipitation bursts seen in panel (a) can be identified in both the integral energy flux and the plasma parallel bulk velocity.
This suggests that the parallel bulk velocity can be used as a proxy for precipitating protons in this context.
In panel (c), the time series of the plasma parallel bulk velocity at the virtual spacecraft shown in Fig. 6 are indicated with a colour code consistent with that of the symbols indicating virtual spacecraft locations.The thick red line corresponds to the parallel velocity at S 1 , i.e., shows the same data as the red line in panel (b).A given fluctuation in the parallel velocity at S 1 can be traced back in time and tailwards by visually identifying its corresponding signature at successive tailward virtual spacecraft.For instance, the parallel velocity enhancement peaking at S 1 around t = 1947 s is consistent with parallel velocity enhancements visible at each of the virtual spacecraft, from the magenta one which is the closest to S 1 till the orange one (S 0 ) which is the deepest in the current sheet, peaking around t = 1927 s.This suggests that some dipolarising flux bundles, which are identified through short-lived enhancements in V x ≈ V in the current sheet, directly lead to proton precipitation into the ionosphere in the regions mapped to the current sheet.However, not all the V enhancements at S 1 can be related to V x enhancements in the current sheet.For instance, the V peak at S 1 around t = 2030 s cannot be traced back to a specific V peak at S 0 .
Table 1 gives the peak times of the integral energy flux, plasma bulk parallel velocity at S 1 and plasma bulk parallel velocity at S 0 for the ten main precipitating flux enhancements between t = 1920 s and t = 2080 s.These times correspond to those of local maxima of each parameter, as indicated with vertical lines in panel (b) of Fig. 7 for the integral energy flux and the parallel velocity at S 1 .For the four precipitation bursts observed at S 1 which could not be traced back to S 0 in the current sheet, the time of S 0 V peak is replaced with an X in the last column.In all cases but one, the peak times for the integral energy flux  Table 1.Times of peak integral energy flux and plasma bulk parallel velocities at S1 and S0 for the main bursts of proton precipitation.In the last column, X indicates whenever a given precipitation burst observed at S1 could not be traced back to S0.In the four cases where a precipitation burst cannot be directly linked to a DFB passing by S 0 , the signatures are lost between the cyan and blue spacecraft (burst number 3), between the blue and indigo spacecraft (burst number 6), between the indigo and purple spacecraft (burst number 8), and between the light and dark green spacecraft (burst number 9).These virtual spacecraft are located in the transition region, which suggests that the transition region can act like a buffer for DFBs, either directly transmitting them and leading to bursts of precipitating protons, or releasing precipitating particles with no direct correlation with an incoming DFB from the tail.

Discussion
This paper presents the first unambiguous investigations of the global terrestrial magnetic field dynamics and its effects on proton precipitation in the midnight sector using a global hybrid-Vlasov model of the near-Earth environment, Vlasiator.Mende et al. (2002) listed four possible causes for particles to precipitate: (i) injection of particles on closed geomagnetic field lines, (ii) interaction of particles with electric fields, in particular waves, (iii) compression of the magnetic flux tube, and (iv) scattering on stretched magnetic field lines.All these causes can be studied using Vlasiator, except for some of the wave-particle interactions, as waves may not all be resolved by the model.It is known that proton precipitation can result from interactions with electromagnetic ion cyclotron (EMIC) waves (Erlandson and Ukhorskiy, 2001;Sakaguchi et al., 2008) as well as magnetosonic waves (Xiao et al., 2014).While magnetosonic waves are well resolved in Vlasiator, it is unclear whether this is also the case for EMIC waves, as their expected spatial scale is close to the ordinary-space grid resolution in the current runs.Wave-particle interactions are expected to be particularly important in the case of loss-cone scattering of trapped particles, i.e., ring current protons in the case of auroral proton precipitation.In this study, however, the analysis of Fig. 7 suggests that the precipitating protons essentially come from the magnetotail and are hence precipitating due to pitch-angle scattering resulting from the small curvature radius of the magnetic field near the neutral sheet (Sergeev and Tsyganenko, 1982).While it would prove interesting to estimate whether there is a non-negligible contribution to this precipitation from wave-particle interactions, this has to be left for a future study.
One limitation coming from the simulation set-up is related to the fact that the simulation run used in this study is 2D in ordinary space, which besides preventing a 3D description of the tail plasma requires using a line dipole instead of a point dipole for the geomagnetic field.The line dipole is scaled to reproduce a realistic dayside magnetopause stand-off distance, but it implies that geomagnetic field line equations have the form r(θ) = C cos θ in polar coordinates, with C a constant and θ the latitude, instead of r(θ) = C cos 2 θ in the case of a centre dipole.This means that it is not possible to map L values with geomagnetic latitudes at ionospheric altitudes in this run; hence, we cannot assess how virtual spacecraft S 1 and S 2 are mapped to the ionosphere in terms of geomagnetic latitude.Such a discussion will be possible once 3D-3V Vlasiator runs are available.
However, this 2D set-up still allows us to investigate the connection between tail dynamics, such as dipolarisation fronts, and precipitation.
Despite those limitations, the precipitating proton fluxes (differential number flux and integral energy flux) observed at the two selected locations (virtual spacecraft S 1 and S 2 ) are in agreement with predictions from the Hardy model (Hardy et al., 1989(Hardy et al., , 1991) ) in the midnight MLT sector and for active geomagnetic conditions (Kp > 3).While direct comparison with particle data from spacecraft would be desirable, the fact that DMSP satellites, which are currently the main source of direct precipitating proton observations, are on Sun-synchronous orbits essentially near the dawn-dusk plane does not allow midnight-sector observations.While velocity distributions are not saved in every simulation cell in Vlasiator runs, results suggest that the parallel component of the plasma bulk velocity can be used to at least qualitatively describe the integral energy flux of precipitating protons (cf.Fig. 7b).The idea that a beam superimposed to the core plasma population could lead to observable effects in the VDF first and second moments (bulk velocity and temperature) was brought forward by Parks et al. (2013), who showed that the apparent slowing down and temperature increase of the solar wind associated with nonlinear structures are due to the presence of a beam of particles propagating in opposite direction and with greater energy than the core solar wind population.In our case, precipitating proton beams along the magnetic field direction affect the parallel component of the plasma bulk velocity in an analogous manner.Future work could therefore include deriving a proxy for proton precipitation relying on plasma bulk parameters, which are saved in every simulation cells, to quantitatively estimate the precipitation parameters such as the integral energy flux or mean energy.
In their study of BBF-associated proton precipitation with a MHD model, Ge et al. (2012) found that dipolarisation fronts led to precipitating proton enhancements into the ionosphere with integral energy fluxes of the order of 0.1 µW m −2 .This corresponds to 6 × 10 4 keV cm −2 s −1 .Assuming a uniform flux along all downwards directions as a rough approximation (2π solid angle), this gives a directional integral energy flux of the order of 10 4 keV cm −2 s −1 sr −1 , which is roughly three orders of magnitude below values obtained with Vlasiator when the global magnetotail reconfiguration is taking place (cf.Fig. 7b).
A likely explanation is that the MHD simulation has low ion temperature T i values compared to our kinetic approach, leading to lower precipitation energy fluxes as these are proportional to T i √ T i in their approach (cf.eq. ( 1) in Ge et al. ( 2012)).We note that the Vlasiator integral flux values are on the other hand in agreement with statistical patterns shown in Galand et al.
(2001) (∼0.1 mW m −2 , equivalent to 10 7 keV cm −2 s −1 sr −1 with the same reasoning as above).One strength of evaluating the proton precipitation parameters using a kinetic model is that not only integral energy fluxes can be calculated, but also differential fluxes, which may enable more detailed future studies of proton precipitation and its link to global magnetospheric dynamics.
The examination of precipitation bursts passing by virtual spacecraft S 1 between t = 1920 s and t = 2080 s suggests that these bursts are associated with dipolarising flux bundles originating from the vicinity of the stable X-line in the current sheet.
Indeed, dipolarisation front signatures in the B z component of the current-sheet plasma can be identified in Fig. 6 of Juusola et al. (2018a) during the studied time interval (written as t = 32:00-34:40 in their figure).While there are not enough such precipitation bursts in the studied time period to carry out a statistical analysis of their properties (integral flux enhancement and parallel velocity enhancement at S 1 , v x enhancement at S 0 ), it can be noted that there does not seem to be a one-to-one correlation between the magnitude of the v enhancement at S 1 and the possibility to trace it back to S 0 .This suggests that the transition region, corresponding to the location in the magnetotail where tail-like geomagnetic field lines become more dipolar (near X = −10.5 R E in Fig. 3), plays a role in regulating these bursts.It is known that fast flows associated with BBFs can bounce when reaching the inner magnetosphere (e.g., Ohtani et al., 2009;Juusola et al., 2013;Nakamura et al., 2013) and may even exhibit multiple overshoots and oscillations around their equilibrium position (Panov et al., 2010).Our results further indicate that contrary to a frequent assumption, the transition region may be more than just a passive mediator between the plasma sheet and the ionosphere.DFBs could perhaps themselves experience some form of bouncing near the transition region, which could explain the regulation of the studied proton precipitation bursts.
Auroral activations concurrent with BBFs have been widely studied (e.g., Nakamura et al., 2001a;Sergeev et al., 2004;Gallardo-Lacourt et al., 2014).However, the time scales that are involved in such processes are of the order of up to 10 min, in contrast to the short-lived DFBs and associated proton precipitation bursts (< 1 min duration) studied here.In terms of auroral emissions observed from the ground, it must be noted that the precipitation bursts observed at S 1 might actually lead to more weakly modulated emissions than can be inferred from the integral energy flux variations at the virtual spacecraft, as the energy dispersion of precipitating protons tends to smooth the integral energy flux as particles get closer to Earth.This warrants future studies involving ground-based optical observations of proton aurora at high enough cadence (a few seconds at most) to investigate whether DFB-related proton aurora signatures can be seen from the ground.If not, this would imply that global magnetospheric simulations, supplemented by data from spacecraft orbiting in the magnetosphere, might be the only tool to investigate the active role played by the transition region in regulating the precipitation of auroral protons to the nightside ionosphere.The evaluation of the differential number flux of precipitating protons at a given location in the nightside magnetosphere is achieved by averaging the velocity distribution function inside the bounce loss cone within energy bins.The integral energy flux of precipitation as well as the mean precipitating energy are also calculated from the differential number flux.The main results of this study can be summarised as follows.
1. From this first case study, we find that Vlasiator reproduces auroral proton precipitation with realistic differential number fluxes and energies, when compared to the Hardy model.
2. During a selected time interval when a single X-line dominates the tail reconnection in the simulation, proton precipitation observed at a virtual spacecraft in the inner magnetosphere occurs in a bursty manner.In this situation, the integral precipitating energy flux exhibits variations mostly similar to the parallel component of the plasma bulk velocity at the virtual spacecraft.This suggests that the integral energy flux can qualitatively be described with the local parallel velocity at locations for which full velocity distributions are not saved in the Vlasiator run.
3. Finally, it is found that proton precipitation bursts can in some cases be traced back to the current sheet and are associated with dipolarising flux bundles.However, not all precipitation bursts correspond to a definite DFB, which suggests that the transition region plays a role in regulating auroral proton precipitation associated with DFBs during BBFs.

Figure 1 .
Figure 1.Simulated plasma temperature in the entire simulation domain, at t = 1800 s.The white area corresponds to near-Earth space located inside the inner boundary at 4.7 RE.Black lines represent magnetic field lines.

Figure 2 .
Figure 2. Geometry used in the derivation of directional differential flux of precipitating protons.

Figure 3 .
Figure 3. (a) Simulated plasma temperature in the nightside magnetosphere at t = 1800 s.The Sun is located beyond the right boundary of the figure.The two black circles filled with white are virtual spacecraft S1 and S2.The arrows indicate the magnetic field directions at S1 and S2.(b) Velocity distribution function of protons at S1 in the (v ,vB×V ) plane.The grid spacing is 1000 km/s.The magenta lines indicate the boundaries of the bounce loss cone.(c) Same for protons at S2.
Figure 3 shows examples of velocity distributions observed at two virtual spacecraft S 1 and S 2 that are used in this study.The main panel (a) shows the proton temperature in the nightside part of the noon-midnight plane at timestep t = 1800 s, with magnetic field lines drawn in black.The white area corresponds to the region of the geospace located earthwards from the inner boundary at ∼ 4.7R E .Black arrows indicate the magnetic field direction at the virtual spacecraft.Panels (b) and (c) show slices of the velocity distributions at S 1 and S 2 , respectively, in the plane defined by the local magnetic field direction v B and the local electric field direction v B×V

v
B×V .Panels (a) and (b) are the distributions at t = 1800.0s at virtual spacecraft S 1 and S 2 , respectively, shown in Fig. 3.The crescent-shaped beam partly in the bounce loss cone shown in panel (a) is associated with a precipitating differential flux narrow in energies (panel e), peaking around 4-5 keV with values close to 10 4 proton cm −2 s −1 sr −1 eV −1 .On the contrary, the hot, nearly-Maxwellian velocity distribution inside the loss cone at S 2 (panel b) is associated to a broad precipitation spectrum Ann.Geophys.Discuss., https://doi.org/10.5194/angeo-2019-59Manuscript under review for journal Ann.Geophys.Discussion started: 24 April 2019 c Author(s) 2019.CC BY 4.0 License.(panelf), with energies of precipitating protons ranging from below 100 eV to nearly 30 keV.The peak differential flux values are of the order of 4000 proton cm −2 s −1 sr −1 eV −1 , at energies of 2-5 keV.The sharp cutoff taking place at ∼25 keV in panel (f) is due to the sparsity threshold, set to 10 −15 m −6 s 3 in this run, hence corresponding to the lower limit of the colour axis in panels a-d.Panel (c) shows a velocity distribution relatively similar to panel (a) but was obtained at virtual spacecraft S 2 with a delay of 100 s, indicating that the region with narrow-beam precipitation moved earthwards.The corresponding differential flux in panel (g) is centred around 2-3 keV, i.e., slightly lower than in panel (e), while peak flux values are almost the same with nearly 10 4 proton cm −2 s −1 sr −1 eV −1 .The last velocity distribution, in panel (d), was taken at S 1 200 s after the one from panel (a).The core population is broader and is partly inside the bounce loss cone, leading to a low flux of low-energy (0.1-1 keV) proton precipitation as can be seen in panel (h).A beam of low phase-space density with a complex structure partly fills the loss cone at velocities of 1300-2100 km/s, which shows in the precipitating spectrum as two peaks, at about 8 keV and 15 keV, with flux values two orders of magnitude lower than in panel (a).These four selected examples suggest that at virtual spacecraft S 1 and S 2 one can observe a variety of precipitating fluxes during the Vlasiator simulation.If compared to observations, these differential fluxes are within the same order of magnitude as the NOAA 6 proton energy spectra shown inBasu et al. (2001) in terms of peak energy (1-10 keV) and slightly higher in terms of flux values (10 2 -10 3 cm −2 s −1 sr −1

Figure 4 .
Figure 4. (top) Examples of proton velocity distribution slices in the local magnetic frame at virtual spacecraft S1 and S2 at various time steps of the simulation.The bounce loss cone is shown with magenta lines.(bottom) Corresponding directional differential precipitating fluxes obtained with the presented method.

Figure 5 .
Figure 5. (a) Differential number flux of precipitating protons observed at virtual spacecraft S1 as a function of time during the simulation.The black line indicates the mean precipitating energy as a function of time.(b) Same but at virtual spacecraft S2.(c) Integral energy flux associated with proton precipitation as a function of time, at virtual spacecraft S1.(d) Same but at virtual spacecraft S2.The green segments in panels (a) and (c) indicate the time interval with precipitation associated with dipolarising flux bundles presented later.

Figure 6 .
Figure 6.x-component of the plasma bulk velocity in the nightside part of the simulation plane at t = 1945.0s.The red circle indicates the location of virtual spacecraft S1, and the eight coloured plus signs show additional virtual spacecraft for which the parallel velocity components are shown in Fig. 7. Black lines represent magnetic field lines.

Figure 7 .
Figure 7. (a) Differential number flux of precipitating protons observed at virtual spacecraft S1 between t = 1920 s and t = 2080 s.The black line indicates the mean precipitating energy as a function of time.(b) Blue line: Integral energy flux associated with proton precipitation as a function of time at virtual spacecraft S1.Red line: Parallel component of the plasma bulk velocity at S1.The vertical lines indicate the times associated with peak values for these two parameters, and the numbers in grey above the panel identify the precipitation bursts discussed in the text.(c) Parallel component of the plasma bulk velocity at S1 and at the virtual spacecraft indicated with plus signs in Fig. 6, with a consistent colour code.
Burst number Time of integral flux peak [s] Time of S1 V peak [s] Time of S0 V peak V at S 1 are the same within 2 s; the exception is burst number 5 for which the triple peak in flux is associated with a single broad peak in V , with 4 s difference in the times of local maxima.
This paper presents the first evaluations of auroral (∼1-30 keV) proton precipitation from a global hybrid-Vlasov magnetospheric model, Vlasiator.The simulation run considered here corresponds to relatively fast solar wind (750 km s −1 ) with purely southward IMF of moderate magnitude (|B| = |B z | = 5 nT).