Vlasov simulation of electrons in the context of hybrid global models: An eVlasiator approach

. Modern investigations of dynamical space plasma systems such as magnetically complicated topologies within the Earth’s magnetosphere make great use of supercomputer models as well as spacecraft observations. Space plasma simulations can be used to investigate energy transfer, acceleration, and plasma ﬂows on both global and local scales. Simulation of global magnetospheric dynamics requires spatial and temporal scales achievable currently through magnetohydrodynamics or hybrid-kinetic simulations, which approximate electron dynamics as a charge-neutralizing ﬂuid. 5 We introduce a novel method for Vlasov-simulating electrons in the context of a hybrid-kinetic framework in order to examine the energization processes of magnetospheric electrons. Our extension of the Vlasiator hybrid-Vlasov code utilizes the global simulation dynamics of the hybrid method whilst modelling snapshots of electron dynamics on global spatial scales and temporal scales suitable for electron physics. Our eVlasiator model is shown to be stable both for single-cell and small-scale domains, and the solver successfully models Langmuir waves and Bernstein modes. We simulate a small test-case section 10 of the near-Earth magnetotail plasma sheet region, reproducing a number of electron distribution function features found in spacecraft measurements. and code description. and TB devised the solver method. MA assisted with data anal- ysis, model development, and comparisons with observations. performed the dispersion tests. MP is the PI of Vlasiator and leads the investigation. YPK, MG, KP, AJ, LT, MD, and MP participated in discussion and ﬁnalization of the manuscript.


Introduction
Physical processes in near-Earth space are dominated by plasma effects such as non-thermal particle distributions, instabilities, plasma waves, shocks, and reconnection. Modern research into space phenomena utilizes both spacecraft measurements and 15 supercomputer simulations, investigating how ions, electrons, and electric and magnetic fields interact in the vicinity of plasma structures. Spacecraft provide point-like observations, limited in their ability to investigate spatial structures, although modern constellation missions can have multiple satellites close by allowing for multipoint analysis to decipher, e.g., current sheet directions (Escoubet et al., 2001;Burch et al., 2016a). Computer simulations on the other hand are limited by spatial resolution, time stepping, and the large difference between ion and electron temporal and spatial scales (see, e.g., Tóth et al., 2017). 20 Simulations capable of modelling the whole near-Earth geospace have historically used magnetohydrodynamics, neglecting kinetic effects and implementing electrons only as e.g. the Hall term correction to Ohm's law. These models can be run for extended periods of time, but as they model plasma motion as a fluid, they use coarse grids, e.g. 0.25 R E (Janhunen et al., 2012) or 0.1 R E  (where R E is the Earth radius) and cannot model kinetic effects but are sufficient for some global dynamics. Recent advances have allowed global investigations into hybrid-kinetic models, where ions are treated 25 as a kinetic self-consistent species and electrons are a charge-neutralizing fluid. Successful approaches include hybrid-Vlasov models (Palmroth et al., 2018) and hybrid-PIC (particle-in-cell) codes (e.g. Lin and Wang, 2005;Sibeck et al., 2008;Omidi et al., 2009;Karimabadi et al., 2014). Kinetic investigation run times rarely exceed one hour or hundreds to a few thousand ion gyroperiods. The simulation spatial resolution is chosen to be relevant to the scales of investigation, with the most usual metric being the ion inertial length d i . The simulation domain must encompass the necessary global dynamics with sufficient space 30 to manage boundary effects.
In order to understand electron physics, kinetic modelling of electrons has been investigated by a number of methods such as full-PIC (ions and electrons as interacting particles, e.g., Hesse et al. 2005), full-Vlasov (ions and electrons as interacting distribution functions, e.g., Umeda et al. 2009;Schmitz and Grauer 2006;Pezzi et al. 2019), hybrid-PIC-electrons (dynamic electron particles, ions as a static background, e.g., Lapenta et al. 2007) and hybrid-Vlasov-electrons (dynamic electron dis-35 tribution function with ions as a static background, e.g., Nunn 2005). In fully kinetic numerical investigations, the standard approach is to alter the ion-to-electron mass ratio of ∼ 1836 to, e.g., 50 (Hesse et al., 2005) or 25 (Wilson et al., 2016) in order to achieve interesting dynamics with available computational resources. Using explicit solvers, resolving waves and kinetic electron instabilities to prevent simulation self-heating requires the spatial resolution to encompass the Debye length λ D (Birdsall and Langdon, 2005) and the time stepping must resolve the electron plasma oscillation ω pe . This can, however, be bypassed 40 via semi-implicit or implicit solver methods, though the resolution decrease incurs the loss of some electron physics. Effects such as the Dungey cycle (Dungey, 1961), involving the whole magnetosphere, are unachievable with traditional kinetic electron approaches. Full-PIC approaches have, however, been applied to investigations of, e.g., reconnection in a Harris current sheet ( Harris 1962, investigated in, for example, Lapenta et al. 2015;Daughton et al. 2011) or asymmetric reconnection . Pritchett (2000) presents a historical review of magnetospheric PIC simulations and anticipates the development model. They note that they do not capture reconnection to an acceptable accuracy and have yet to publish global simulation results. Global ten-moment results for the Hermean magnetosphere have been published in Dong et al. (2019). Furthermore, approaches which focus on electron effects at lower frequencies (negating effects at plasma oscillation timescales) have been investigated in, for example, Lin and Chen (2001) and Tronci and Camporeale (2015).
Several processes that occur in the magnetosphere that depend on electron behavior are still poorly understood. Recently, 70 missions such as Magnetospheric MultiScale (MMS; Burch et al., 2016a) have enabled plasma measurements that are able to better resolve electron-scale physical processes. MMS in particular has provided data to many publications on magnetic reconnection (e.g., Burch and Phan, 2016;Phan et al., 2018;Huang et al., 2018;Hoilijoki et al., 2019a;Fargette et al., 2020), the most popular topic of electron physics investigations. Reconnection-driven jets and dipolarization fronts cause magnetic flux pileup and excitation of waves such as whistlers, affecting energy conversion and dissipation (Khotyaintsev et al., 2011;75 Breuillard et al., 2016;Zhang et al., 2018). Bulk flows along the tail lead to electrons heating up as they approach the Earth (Runov et al., 2015;Artemyev et al., 2013) with the electron-to-proton temperature ratio approaching even 1 (Wang et al., 2012). These flows interact with strong currents found in the plasma sheet (Nakamura et al., 2008;Artemyev et al., 2017).
The dynamics of electrons near the current sheet include strong Hall fields and current sheet thinning (Lu et al., 2019(Lu et al., , 2016. Electron anisotropies can excite electron-driven waves and time-domain structures, such as have been observed recently in 80 different regions of the magnetosphere (e.g., Cattell et al., 2005;Mozer et al., 2015;Ergun et al., 2016). They have been characterized as whistler mode waves, electrostatic solitary waves and lower hybrid waves among other types. These waves interact strongly with electrons, causing effects such as heating, changes to temperature anisotropy, and particle energization.
These energized electrons can then add to energetic particle precipitation, leading to the generation of auroras (Ni et al., 2016). This paper introduces an alternative, novel method for simulating electron distribution function physics in the context of 85 global ion-determined fields. The aim is to investigate how much of the global electron physics and distribution functions can be understood by utilising ion-generated field as modelled by hybrid-kinetic codes, as opposed to a numerically unfeasible global full-kinetic approach. The paper is organised as follows. In Section 2, we introduce the ion-kinetic hybrid-Vlasov code Vlasiator and how the Vlasov equation is solved. In section 3 we introduce the eVlasiator modifications implemented for the analysis of electron distribution functions. Section 3.1 describes how our electron simulation is set up from fields and moments 90 modelled by an ion-kinetic simulation and Section 3.2 details the field solver changes implemented. Section 3.3 describes the sample test simulation used in this study. In Section 4 we perform rigorous validation and stability tests for our electron solver, and in Section 5 we present some electron distribution functions achieved by running our solver on a test dataset, comparing them with existing literature. Finally, Section 6 draws conclusions of our analysis and lays out a plan for future research approaches.

95
2 The Vlasiator ion-kinetic hybrid-Vlasov code Vlasiator (von Alfthan et al., 2014;Palmroth et al., 2018) is, at the present time, the only hybrid-Vlasov code capable of simulating the global magnetospheric system of the Earth, accounting for ion-kinetic effects on spatial and temporal scales which model both magnetopause and magnetotail dynamics. Vlasiator solves the Vlasov equation for particle distribution functions discretized on cartesian grids, with closure provided by Ohm's law augmented by the Hall term. Each particle 100 population is described using a uniform cartesian three-dimensional velocity space grid (3V) with a resolution chosen to accurately model the solar wind inflow Maxwellian distribution and with extents chosen to encompass heated ion populations associated with the magnetosheath and flux transfer events. A standard Vlasiator global run proton velocity-space grid has a resolution of 30 km s −1 , extending between ±4020 km s −1 . To constrain computational cost and memory usage, those parts of the velocity distribution function which have a phase-space density below a sparsity threshold are discarded, except for buffer 105 regions which allow the correct growth of the VDF in these parts (von Alfthan et al., 2014). The proton sparsity threshold is usually set to a value between 10 −17 and 10 −15 m −6 s 3 .
In the spatial domain, Vlasiator can be run in 1D, 2D, or 3D, with 2D the most usual choice in order to evaluate global dynamics. Simulations have used spatial resolutions of, e.g., ∆x = 228 km or ∆x = 300 km, enough to accurately model ion cyclotron waves though not resolving the ion inertial length in all regions of the simulation domain. Large-scale global 3D 110 runs will be made possible in the near future by adaptive mesh refinement (AMR), using non-uniform cell sizes in the spatial domain, thus cutting down on the computational cost of the simulation.
Vlasiator models standard collisionless space plasmas dominated by protons but can also model other particle species in the same self-consistent simulation. However, until now, the electron population has been treated in the usual way of implementing it as a massless charge-neutralizing fluid. This method does not track the evolution of electrons beyond assuming charge 115 neutrality, and therefore, these standard Vlasiator simulations cannot be used to infer electron dynamics. This paper presents a novel approach for investigating how a global plasma current structure can influence electrons with limited self-consistency enforced through plasma oscillation and current densities.

Solving the Vlasov equation
Vlasiator uses the hybrid-Vlasov ion approach to model the near-Earth space plasma environment. The full six-dimensional 120 (6D) phase space density f s (x, v, t), with x the ordinary space variable, v the velocity space variable, and t the time variable, for each ion species s of charge q s and mass m s is evolved in time using the Vlasov equation (1). The electric and magnetic fields, denoted E and B respectively, are evolved using Maxwell's equations: Faraday's Law (2), Gauss's Law (3) and Ampère's Law (4), in which µ 0 and ε 0 are the vacuum permeability and permittivity, respectively, and j is the total current density.
The solenoid condition in Gauss's Law (3) is ensured via divergence-free magnetic field reconstruction (Balsara, 2009). In 125 the hybrid approach, electrons are assumed to maintain plasma neutrality, resulting in the charge density ρ q in Gauss's law vanishing. In the Darwin approximation, also used in many hybrid codes, propagation of light waves is neglected by removing the displacement current term ε 0 ∂E ∂t in Ampère's law (4). The Vlasiator field solver follows the staggered-grid approach of Londrillo and Del Zanna (2004), and is described in detail in Palmroth et al. (2018).
The generalized Ohm's Law providing closure for the Vlasov system is

135
where V is the plasma bulk velocity, σ is the conductivity, e is the elementary charge, n e is the electron number density, and P e is the electron pressure tensor. In hybrid approaches of collisionless plasma, we can assume high conductivity, neglecting the first term on the right-hand side. In the limit of slow temporal variations, the electron inertia term (the last term on the right-hand side) also vanishes. The remaining two terms on the right-hand side of the equation are the Hall term, J × B/(n e e), and the electron pressure gradient term, ∇ · P e /(n e e). In hybrid models, a true description of electron pressure is unavailable 140 so it must be described via some approximation such as adiabatic, isothermal or polytropic electrons or a fixed ion-to-electron temperature ratio, or by neglecting the small electron pressure gradient term altogether. The standard ion-hybrid Vlasiator code supports isothermal fluid electrons but existing simulations have always set this temperature to zero. This along with assuming charge-neutrality (proton number density n p = n e ) results in the ion-hybrid Vlasiator using the simplified MHD version of Ohm's Law with the Hall term included: As Vlasov methods do not propagate particles but rather evolve distribution functions, we now briefly explain the semi-Lagrangian method employed by Vlasiator (Palmroth et al., 2018). Vlasiator propagates distribution functions of particles following the SLICE-3D method (Zerroukat and Allen, 2012) and utilizing Strang splitting with advection (the second term of Vlasov's equation 1) and acceleration (the third term of Vlasov's equation 1) calculated one after the other with a Leapfrog 150 offset of 1 2 ∆t. In this manuscript, ∆ denotes steps on the full simulation grid and associated time step and δ is used to indicate calculations performed as sub-stepping. For each time step, a Vlasov acceleration is evaluated with time step length ∆t which is, amongst other things, limited to a maximal Larmor orbit gyromotion rotation of 22°. For each acceleration step, a transformation matrix is initialized as an identity matrix. The transformation matrix decomposes rotation and field-parallel acceleration into three shear transformations. The matrix is updated with substepping of δ ts where each δ ts corresponds to 155 a 0.1°Larmor gyration. Instead of applying linear acceleration by electric fields, a method similar to the Boris-push method (Boris, 1970) is applied, where first a transformation is performed to move to a frame in which electric fields vanish, then the rotation is applied, and then a frame transformation back to the original frame is added. In the standard hybrid formalism, the frame without electric fields is found via the MHD Ohm's law with the Hall term included (6). This Hall frame estimates the frame of reference of electrons, assuming electrons generate a current density which corresponds to the local magnetic field 160 structure, in accordance with Ampère's law. After substepping is evaluated, the transformation matrix is applied to the gridded velocity distribution function by the SLICE-3D algorithm.

Simulation initialization
Modelling the evolution of electron distribution functions in response to global magnetic field structures requires input from the large-scale fields and moments of a Vlasiator simulation of near-Earth space. In the eVlasiator approach, we read magnetic field vectors and proton plasma moments for the chosen simulation domain and apply user-defined temperature scaling to generate initial Maxwellian electron velocity distribution functions. 175 We do not model electrons throughout the whole global domain, choosing instead a region of interest to reduce the computational cost, though our method is designed to work with any subset of and up to the whole global domain. For the selected domain, we read in the Vlasiator ion-hybrid simulation proton moments, cell-face-average magnetic field components and cell-edge-average electric field components (the latter being used by the staggered-grid field solving algorithm from Lon-

The eVlasiator field solver
For non-relativistic electrons, the advection step of solving the Vlasov equation requires no adjustments, but the acceleration step must be adjusted to account for electron oscillation and the relevant electric fields.
In our field solver we maintain static magnetic fields as read from the input Vlasiator simulation. We model the electric field by including additional terms in Ohm's law (5), allowing for interaction of distribution functions with electron-oscillation 190 electric fields. Whistler mode propagation is not included in this study. We do not include any electric field due to Gauss' law.
We will consider each term of the eVlasiator field solver separately: -As we keep magnetic fields static, we do not implement Faraday's law (2).
-Collisionless plasma physics assumes that electrons are fast enough to balance out any charge imbalance, and in hybridkinetic simulations this holds true. We do not implement Gauss' law (3) in order to quantify to what extent charge 195 neutrality holds in the eVlasiator context.
-The last term in Ampère's law (4) is the displacement current, which is neglected in the Darwin approximation. However, electron motion can be very rapid and thus we now include this term in our model, though still maintaining static magnetic fields. This approach thus constrains electrons to the defined static magnetic fields and does not introduce light waves.

200
-As our plasma remains collisionless, we maintain our assumption of infinite conductivity, and thus the J/σ term in the generalized Ohm's law (5) remains zero.
-The Hall term, J × B/(n e e), is used to estimate the electron reference frame, and is discussed further below.
-As eVlasiator models electrons with full distribution functions, we include the full electron pressure tensor P e and implement the electron pressure gradient term using spatial gradients calculated for electron pressure.

205
-The final term of the general Ohm's law is the electron inertia term. Much like with our choice of including the displacement current, we now include the electron inertia term in our solver.
For electron dynamics to be modelled, electron gyration and plasma oscillation must both be considered. We limit the acceleration time step to a maximum of 22°of Larmor rotation or 22/360 of a single plasma oscillation. The value of 22°is used to ensure our VDF remapping algorithm SLICE-3D remains stable and the value 22/360 was chosen for equal resolution of 210 both oscillations. Due to the computational cost of SLICE-3D remapping, a substepping approach is used in order to accurately model the electron gyromotion and plasma oscillation. The electron gyroperiod is τ ce = 2πω −1 ce and the plasma oscillation time is τ pe = 2πω −1 pe , where the electron plasma frequency is In transformation matrix generation, substepping is constrained to a maximum of δ ts ≤ min(τ pe , τ ce )/3600. For each substep, a procedure similar to the hybrid method is applied, with the improvement that instead of performing gyration in the Hall frame (estimating the electron frame within the hybrid context) the gyration is performed in the actual substep-updated electron frame.
Electron oscillation is handled in parallel with gyration by tracking an additional cell-volume-averaged electric field compo-220 nent E Je which is itself derived from the small-scale electron oscillation. For each substep, we perform two parallel 4th order Runge-Kutta propagations. The first one is tracking electron bulk velocity response δV e to the E Je field. This simple acceleration term is in fact equal to evaluating current changes via the electron inertia term in Ohm's law with the E Je field included in the left-hand-side electric field. The second 225 Runge-Kutta propagation tracks the evolution of the E Je field due to changing current density, according to the displacement current on the right-hand side of Ampère's law (4), whilst maintaining static magnetic fields. The ∇ × B term in Ampère's law is fixed to the static input magetic fields. Thus, for each Runge-Kutta step, the electric field E Je is updated with where c is the speed of light, and B, n p and the proton bulk velocity V p are assumed constant throughout the substep. Each of the four δV e Runge-Kutta coefficients are updated with the latest estimate for δE Je , and vice versa. Values for E Je are stored between acceleration steps to ensure continuity of the oscillation. The change δV e calculated via each Runge-Kutta step is then applied to the transformation matrix, allowing the solver to proceed to perform gyration in the electron frame of reference. The 235 substepping procedure is visualized in Figure 1. Further details of the solver and advection methods in Vlasiator can be found in Palmroth et al. (2018).
After substepping is completed, the transformation matrix describing Vlasov acceleration is passed to the SLICE-3D algorithm, which decomposes the transformation into three cartesian shears and updates the velocity distribution function for the particle species.

Transformation matrix generation substepping
Dual RK4 for each substep line dipole scaled to result in a realistic magnetopause standoff distance (Daldorff et al., 2014). The simulation has an inner boundary at 3 · 10 6 m ≈ 4.7 Earth radii, modelled as a perfectly conducting sphere.

250
For this eVlasiator sample run, we choose a region from the magnetotail with 70 × 1 × 40 simulation cells in the X, Y, and Z directions, respectively. The subregion extent is from X − = −75.6 · 10 6 m to X + = −54.6 · 10 6 m, from Y − = −0.15 · 10 6 m to Y + = +0.15 · 10 6 m, and from Z − = −6 · 10 6 m to Z + = +6 · 10 6 m. Within this domain, visualized with a small rectangle in Figure 2a, the electron plasma period τ pe ranges from ∼ 0.7 ms in the magnetotail plasma sheet up to ∼ 2.5 ms in the near-sheet X-line site. The distributions are discretized onto eVlasiator velocity meshes, with the electron velocity mesh consisting of 400 3 cells, extending from −4.2·10 7 to +4.2·10 7 m s −1 in each direction, resulting in an electron velocity space resolution of 210 km s −1 .
The electron VDF sparsity threshold was set to 10 −21 m −6 s 3 , ensuring good representation of the main structure of the VDF.
Discretizing a hot and dense electron distribution onto a cartesian grid is numerically challenging without using vast amounts of 260 memory. As portions of our simulation domain have proton temperature up to 10 8 K, we use an empirical estimate of T i /T e ∼ 4 as magnetosheath temperature ratios are usually around 4 to 12 (Wang et al., 2012). Paterson and Frank (1994), Hoshino et al. (2001), Artemyev et al. (2011), and Grigorenko et al. (2016) show similar proton-electron temperature ratios in the magnetotail.
In order to constrain the extent of our velocity space and numerical requirements of our solver, we implement our electrons with a mass of 10 times the true electron mass, resulting in an ion-to-electron mass ratio of m i /m e = 183.6. As mentioned 265 above, we calculate the required electron bulk velocity for each cell using the local volumetric (cell-average) derivatives so that the ion and electron fluxes in each cell correspond with the current density J required for fulfilling Ampère's law (4) (with the displacement current neglected at initialization). This is equal to performing a transformation to the Hall frame of reference.
Proton densities, magnetic field lines, proton temperatures, proton bulk velocities and electron bulk velocities calculated for simulation initialization are shown in Figure 2 along with an overview of the input Vlasiator simulation and the selected 270 electron sub-domain.

Single-cell stability of electron oscillation
To validate the performance of our electron solver, we performed single-cell tests, with resultant electron bulk velocities V e and plasma oscillation electric fields E Je shown in Figure 3. These single-cell tests did not have magnetic field curvature or an 275 ion population present, resulting in the electron motion oscillating around a stability point at V e = 0 and E Je = 0. We set the electron number density to n e = 0.1 cm −3 and the magnetic field to B x = 20 nT (panels a through d) or B x = 200 nT (panels e and f). We set an initial velocity perturbation of V e,0 = (−100, −150, 200) km s −1 , close to but below our electron velocity resolution of ∆v = 210 km s −1 . As can be seen from Figure 3, the electron oscillatory motion is well resolved and remains stable over an extended period. In panels e) and f) where the magnetic field strength was artificially increased in order to set 280 the plasma and gyroperiods to values closer to each other (1.11 ms and 1.79 ms, respectively), we see a gradual evolution of oscillation amplitude and, thus, E Je field magnitude as the two types of electron motion interact. Over longer periods of time this growth becomes unstable, but it can be counteracted by using a smaller substep time step. Also, this instability occurs only when τ ce ≈ τ pe which does not occur in our full simulation domain.

Dispersion relation analysis 285
Although our method is geared towards solving electron motion at coarse spatial resolutions, to further validate the solver, a wave dispersion test was run (Kilian et al., 2017;Kempf et al., 2013). As waves are a collective, emergent phenomenon of an electron number density of n e = 0.4 · 10 6 m −3 , an electron temperature of T e = 2.5 MK, and a magnetic field magnitude of 50 nT. In one simulation, the magnetic field direction was chosen to coincide with the extended simulation direction (resulting in parallel plasma wave modes to be resolved), in the other one, the magnetic field was set up perpendicular to the long dimension, resulting in perpendicular mode resolution. The plasma had zero bulk velocity in the simulation frame, with an added white noise velocity fluctuation ofṽ = 1000 m/s. The simulation was run for 0.037 seconds (433 ω −1 pe ).
295 Figure 4 shows the dispersion data resulting from spatial and temporal Fourier transform (using a von Hann window).
Overlaid are analytic dispersion curves for the Langmuir wave (black dashed curve) and electron Bernstein modes (black solid curves). The wave behaviour in the simulation shows good agreement in both parallel and perpendicular directions. One noteworthy additional feature visible in the parallel direction (Figure 4a) is the presence of an entropy wave feature at low wave number k and angular frequency ω that shows a quantization consistent with the electron velocity space resolution.

Stability within larger simulation domain
We also evaluate the stability of our solver over the larger simulated domain described in Section 3.3, with initialization values derived from the Vlasiator hybrid-Vlasov simulation. These graphs are shown in Figure 5. Panels a through e show the evolution of electron temperature values over a simulation of 1.0 s, covering hundreds of electron plasma periods and, for the most part, tens of gyroperiods. We evaluate minimum, maximum, mean, and median values for total, B-parallel, and B-perpendicular 305 electron temperatures. The system is seen to relax somewhat towards a final state, though some evolution is still apparent at the end of the simulation, possibly due to boundary effects. The maximum temperature plot in panel b is of particular interest as the hottest plasma cells appear to diffuse into their surroundings until t ∼ 0.4 s when dynamic gyration processes overtake this temperature diffusion with perpendicular heating.
Panel f shows the agyrotropy measure (Swisdak, 2016) calculated from the electron pressure tensor, indicating that in 310 the majority of the simulation domain electrons remain gyrotropic and even peak values do not grow past 10 −3 . Panel g shows statistics for the electron number density deviation from the initialisation value, indicating loss of plasma neutrality due to the motion of electrons. The minimum value oscillating between approximately 10 −9 and 10 −6 cm −3 indicates the level of numerical fluctuations, and the maximum, mean and median values show how charge imbalance does grow initially but stabilises within about 0.1 s and does not grow beyond 10 −1 cm −3 .

315
In panels h through k of Figure 5 we show how the instantaneous plasma oscillation electric fields E Je are well-behaved throughout the simulation box, converging towards stable values. We note that as E Je fields oscillate around zero, the averages are indeed zero throughout (not shown) and the values used for inferring minimum, maximum, mean and median values are instantaneous values from a arbitrary phase of the oscillation. In panel l we show the normalized current density J departure from the balance current J B = ∇×B µ0 which would be required to maintain the magnetic field structure according to Ampère's 320 law (4). This metric is seen to also stabilize, mostly at values well below unity. We expect the maximum value outliers to be due to locally small values of J B . Panels m and n show statistics for the parallel and perpendicular components of the electric field caused by electron pressure gradients, that is, the − ∇·Pe nee term. As expected due to the ability of electrons to propagate along field lines, perpendicular components are much larger than parallel components. All components remain stable at roughly their initial values. A minimum value is not shown as the use of a numerical slope limiter in the calculation of pressure gradients 325 gives identically zero field components at local extrema of pressure.
As part of our evaluation of solver stability, we performed a comparison run where our electron solver performed the rotation transformation corresponding with gyromotion in the Hall frame instead of in the substep-associated electron bulk frame. This transformation choice resulted in unstable growth of, in particular, E Je , as could be expected (not shown).

330
Results from the electron simulation after 1.0 s of evolution are presented in Figure 6. Figures 6a, b show parallel and perpendicular acceleration or deceleration of electrons as the ratio of end-of-simulation temperature to initial temperatures. Heating is found in particular near the X-line configuration and where the PSBL meets the magnetosphere, with parallel heating more localized than perpendicular heating. Figure 6c shows the agyrotropy measure (Swisdak, 2016) calculated for electrons, indicating where the electron distribution 335 has become non-gyrotropic. In most of the simulation domain, the value is very small, but enhanced agyrotropy (still relatively small values below 10 −3 ) are found in the PSBL regions and at the magnetic field X-line. Some of this agyrotropy may be due to spatial sampling of electron gyromotion with a magnetic field gradient leading to larger gyroradii further away from the plasma sheet. Figure 6d shows the electric field due to ∇P e , with the field strongest where the PSBL meets the magnetosphere. The field 340 direction is pointed towards the tail sheet or the magnetosphere, as expected. Magnitudes remain of the order of a few millivolts per metre.
Figures 6e,f quantify the charge imbalance resulting from electrons evolving due to static magnetic fields and the electric field resulting from the Ohm's law terms presented in this paper. Figure 6e, shows the level of charge imbalance as change in electron number density, and Figure 6f as the change scaled by the original electron number density. In the majority of Figure 5. Evolution of electron and solver parameters over the whole simulation domain. a-d: Minimum, maximum, mean, and median values for electron temperature Te and its components parallel and perpendicular to the local magnetic field. e: Minimum, maximum, mean, and median values for electron temperature anisotropy. f: Minimum, maximum, mean, and median values for electron agyrotropy QAg,e. g: Minimum, maximum, mean, and median values electron density deviation from initial state, indicating charge imbalance. h-k: Minimum, maximum, mean, and median values for the plasma oscillation electric field EJ e and its components parallel and perpendicular to the local magnetic field. l: Minimum, maximum, mean, and median normalized values for current density J deviation from the value J B = ∇×B µ 0 required to fulfill Ampère's law for the local magnetic field. m,n: Maximum, mean, and median values for parallel and perpendicular components of the electric field due to electron pressure gradients. the simulation domain, imbalance remains below 10 −2 cm −3 . The electric field response is unable to maintain full plasma neutrality with some regions near the magnetosphere showing greater deviation from the initial state. Some stronger imbalance at the domain edges is likely a boundary effect which shall resolve itself with a larger simulation domain.
In Figure 7 we display electron velocity distribution functions after 1.0 s of simulation. Figure 7a shows the evolved electron temperature anisotropy T ⊥,e T −1 ,e , and Fig. 7b displays the maximum of instantaneous values of E Je , taken over 10 measure-   of the magnetosphere. As we have bulk flows of both ions and electrons towards the tail current sheet, some small part of this heating can be attributed to betatron acceleration as electrons convect towards stronger magnetic fields just adjacent to the actual high-beta plasma sheet. Other effects causing anisotropies may arise from spatial leakage of electrons undergoing plasma oscillation, with gyromotion binding perpendiculary heated electrons to the oscillation region and parallel accelerated electrons propagating along field lines to the near-magnetosphere PSBL regions.
The maximum of instantaneous values of E Je , shown in Figure 7b, indicate that the strongest electron oscillations on our 360 simulated scales are found in or near the PSBL, which would be consistent with observations of electron-driven waves in the PSBL (Onsager et al., 1993). Some increase in E Je is seen also at the X-line location, but not in other parts of the current sheet.
We note that the X-line included in the Vlasiator simulation snapshot was not actively reconnecting. Comparison with Figure 7a and virtual spacecraft measurements indicate that parallel features, akin to electron beams, are indeed found in regions with enhanced E Je . of hot plasma, but rather shows that the acceleration may not be strong enough to be discerned over the main hot VDF.
Parallel heating near the magnetotail plasma sheet has been reported to coincide with bi-directional electron distributions (Hada et al., 1981) with temperature ratios going up to 2-3, as in our simulation. Our VSC [2] and [5] show clear bi-directional distributions. Due to our static background magnetic field, our parallel heating cannot be due to conventional Fermi acceleration. However, Hada et al. (1981) propose that adiabatic plasma processes where curvature drifts dominate over gradient drifts 375 (Yamamoto and Tamao, 1978) can lead to significant parallel heating. Our VSC [1] is from close to the X-line and shows parallel elongation of the central part of the distribution, reminiscent of the football or shifted-football distributions of Figure 2 of Hoshino et al. (2001). Asano et al. (2006) describe streaming 500 eV electrons at the PSBL, associated with a substorm event and variation of B y , especially at small scales. Scaling with our electron mass, this corresponds to approximately 4000 km s −1 electron velocities, 380 which is reasonably within the range of our VDFs in Figure 7. We note that our simulation produces a background B y profile with ∇B y in agreement with Figure  We also note our VSC [3] displaying a disjoint parallel beam, matching the ISEE-2 observations in Figure 5 of Onsager et al. (1991).
Observations of perpendicular crescents are shown in MMS data in Burch et al. (e.g. 2016bBurch et al. (e.g. , 2019 at electron diffusion regions (EDRs), in conjunction with dayside magnetopause reconnection sites. These observed structures are produced at very 390 small spatial scales, not captured by our current model. We do, however, observe similar agyrotropic crescents in our results further out (in particular in Figure 7j), suggesting successful capture of a level of electron dynamics. These perpendicular crescents are found at very low phase-space density values, as could be expected by the low agyrotropy values seen in Figure 6e.
Something akin to a parallel electron crescent (Burch et al., 2016b) can be seen in Figure 7c, and bi-directional distributions as reported in Figures 6 and 7 of Burch and Phan (2016) are qualitatively similar to our Figures 7k and m.

Conclusions
In this method paper we have presented a novel approach to investigating electron distribution function dynamics in the context of global ion-hybrid field structures. Our method exploits global dynamics provided by hybrid-Vlasov simulations in order to evaluate the response of gyrating and plasma oscillating electrons to global magnetic field structures.
We have shown our solver to behave in a stable manner, resolving electron inertia and updating a responsive electric field 400 E Je derived from the displacement current. If run at much finer spatial resolutions, our model replicates Langmuir waves and electron Bernstein modes. Electron temperatures evolve in response to the field structure but do not experience uncontrolled growth. Our sample simulation produces multiple features associated with spacecraft observations of VDFs, such as parallel acceleration, bi-directional distributions, and perpendicular crescents.
Our model has several built-in limitations as it does not treat electrons as a fully self-consistent species. Magnetic fields 405 gathered from the Vlasiator simulation are kept constant and thus force electron bulk motion to adhere to the required current density structure. As the initialization information is gathered from a hybrid-Vlasov simulation, it has a spatial resolution far below that required for resolving electron-scale waves such as whistlers, Bernstein waves and chorus waves. Scattering of electrons via these missing waves is somewhat accounted for by initializing every simulation from a Maxwellian isotropic distribution. These features together limit the applicability of the model to short periods of time. On the other hand, our model 410 is efficient, and much larger spatial domains of investigation are easily achievable. Also, multiple eVlasiator runs can be performed from a single Vlasiator magnetosphere run to evaluate different driving conditions such as temperature ratios and anisotropies. The method builds on the efficiently parallelized Vlasiator codebase and will benefit from future numerical and computational improvements to Vlasiator solvers.
Our model can be applied to investigate electron dynamics on global spatial scales, with the current version applicable to 2D 415 investigations, e.g., in the noon-midnight meridional plane. Electron velocity distribution functions generated by the model can be used to investigate, e.g., energetic electron precipitation into the Earth's auroral regions. The generated electron anisotropies can be used to infer regions where, for example, whistler waves can be expected to grow. The model can be run for several different initialization time steps to evaluate long-term evolution of precipitating electron distributions. This could be used to, for example, evaluate electron distribution changes as bulk flows and dipolarization fronts in the Earth's magnetotail propagate 420 earthward. Li et al. (2020) observe electron Bernstein modes driven by perpendicular crescent distributions. As we have shown in Figures 4 and 7, with sufficient resolution we can reproduce electron Bernstein waves and agyrotropic electron distributions.
Thus, we are in position to investigate this connection further in eVlasiator.
Future improvements to our model will allow simulation initialization from non-uniform 3D-3V Vlasiator meshes, allowing investigation of spatially three-dimensional topologies including tail plasma sheet clock angle tilt. A possible path of future 425 investigation would be to upsample the initialization fields and moments in order to achieve better resolution, but we emphasize that the model does not attempt to solve electrons in a fully self-consistent manner as magnetic fields are still kept constant. Increasing resolution by interpolating the input moments to a finer grid might not significantly improve plasma sheet density and temperature profiles. Increasing spatial resolution introduces numerous caveats including increased computational cost and possible charge imbalance resulting from spatially resolved electron oscillations, though our dispersion tests did not