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 flows on both global and local scales. Simulation of global magnetospheric dynamics requires spatial and temporal scales currently achievable through magnetohydrodynamics or hybrid-kinetic simulations, which approximate electron dynamics as a charge-neutralizing fluid. 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 of the near-Earth magnetotail plasma sheet region, reproducing a number of electron distribution function features found in spacecraft measurements.

1 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 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 for example Tóth et al., 2017).

Simulations capable of modelling the whole near-Earth geospace have historically used magnetohydrodynamics, neglecting kinetic effects and implementing electrons only as the Hall term correction to Ohm's law for example. 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 (Wang et al., 2020) (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 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 1 h 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 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 distribution 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 for example 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 via semi-implicit or implicit solver methods. If such an approach is used and the resolution is decreased, selecting a very low resolution may result in 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 for example reconnection in a Harris current sheet (Harris, 1962, investigated in, for example, Lapenta et al., 2015; Daughton et al., 2011) or asymmetric reconnection (Hesse et al., 2016). Pritchett (2000) presents a historical review of magnetospheric PIC simulations and anticipates the development of more realistic, global 3D magnetosphere models with increasing computational resources.

More recent simulation studies of electron physics in the magnetosphere such as the PIC simulations by Bessho et al. (2014, 2016) and Hesse et al. (2016) have focused on local regions, modelling for example electron diffusion regions (EDRs) and extracting resultant electron velocity distribution functions (eVDFs). Liu et al. (2013) investigated the small-scale three-dimensional structure of EDRs with a realistic proton–electron mass ratio with a small configuration, and extended to a larger local 3D configuration with a reduced proton–electron mass ratio. These simulations are always local with prescribed driving. A more global approach, MHD-EPIC, is presented by Daldorff et al. (2014), with a two-way coupling of a global, 2D Hall MHD magnetosphere model and local implicit PIC model at regions of interest, where a proton–electron mass ratio of 25 was used. Notably, these PIC regions handled by implicit solvers do not resolve the Debye length. MHD-EPIC has since been employed to study the magnetosphere of Ganymede in 3D with large embedded PIC domains by Tóth et al. (2016) and Zhou et al. (2019).

An example of small-scale global electromagnetic implicit PIC modelling for a weak comet has been performed by Deca et al. (2017, 2019) with a reduced proton–electron mass ratio of 100, and local simulations for a lunar minimagnetosphere (Deca et al., 2015) with a reduced proton–electron mass ratio of 256.

Ricci et al. (2002) discuss the effect of the ratio between the proton mass m_p and the electron mass m_e as a part of the GEM challenge, concluding that reconnection rates are well captured by smaller mass ratios of $m_{\mathrm{p}}/m_{\mathrm{e}}=\mathrm{180}$, although with modified electron kinetics. Lapenta et al. (2010) discusses modifications to electron microphysics at reconnection sites in more detail in relation to proton–electron mass ratios of 64, 256, and 1836 using an implicit PIC model.

Another approach compared to PIC simulations is to represent particle velocity distributions with moments beyond the MHD approach (Wang et al., 2015). For example, Huang et al. (2019) have developed a six-moment multi-fluid full-Maxwell model. They note that they do not capture reconnection to an acceptable accuracy and have yet to publish global simulation results. Global 10-moment results for the Hermean magnetosphere have been published by Dong et al. (2019). Furthermore, approaches which focus on electron effects at lower frequencies (neglecting effects at plasma oscillation timescales) have been investigated by, for example, Lin and Chen (2001) and Tronci and Camporeale (2015).

Several processes that occur in the magnetosphere that depend on electron behaviour are still poorly understood. Recently, 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; 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, 2016). Electron anisotropies can excite electron-driven waves and time-domain structures, such as have been observed recently in 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 global ion-determined fields. The aim is to investigate how much of the global electron physics and distribution functions can be understood by utilizing ion-generated field as modelled by hybrid-kinetic codes, as opposed to a numerically unfeasible global full-kinetic approach. The paper is organized as follows. In Sect. 2, we introduce the ion-kinetic hybrid-Vlasov code Vlasiator and how the Vlasov equation is solved. In Sect. 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 modelled by an ion-kinetic simulation. Section 3.2 describes the time propagation of the eVDF, and Sect. 3.3 details the field solver changes implemented. Section 3.4 describes the sample test simulation used in this study. In Sect. 4 we perform rigorous validation and stability tests for our electron solver, and in Sect. 5 we present some electron distribution functions achieved by running our solver on a test dataset, comparing them with existing literature. Finally, Sect. 6 draws conclusions of our analysis and lays out a plan for future research approaches.

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 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⁻¹, extending between $\pm \mathrm{4020}\mathrm{km}{\mathrm{s}}^{-\mathrm{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 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⁻¹⁷ and ${\mathrm{10}}^{-\mathrm{15}}{\mathrm{m}}^{-\mathrm{6}}{\mathrm{s}}^{\mathrm{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 Δx=228 km or Δx=300 km for example, 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 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. The method does not track the evolution of electrons beyond assuming charge 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.

2.1 Solving the Vlasov equation

Vlasiator uses the hybrid-Vlasov ion approach to model the near-Earth space plasma environment. The full six-dimensional (6D) phase space density $f_{\mathrm{s}}(\mathit{x},\mathit{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 Eq. (1). The electric and magnetic fields, denoted E and B respectively, are evolved using Maxwell's equations: Faraday's law (Eq. 2), Gauss's law (Eq. 3) and Ampère's law (Eq. 4), in which μ₀ and ε₀ are the vacuum permeability and permittivity, respectively, and j is the total current density. The solenoid condition in Gauss's law (Eq. 3) is ensured via divergence-free magnetic field initial-condition reconstruction (Balsara, 2009). In 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 ${\mathit{\epsilon }}_{\mathrm{0}}\frac{\partial \mathit{E}}{\partial t}$ in Ampère's law (Eq. 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).

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial t}}+\mathit{v}\cdot {\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial \mathit{x}}}+{\displaystyle \frac{q_{\mathrm{s}}}{m_{\mathrm{s}}}}\left(\mathit{E}+\mathit{v}\times \mathit{B}\right)\cdot {\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial \mathit{v}}}=\mathrm{0}\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\times \mathit{E}=-{\displaystyle \frac{\partial \mathit{B}}{\partial t}}\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{B}=\mathrm{0}\mathrm{and}\mathrm{\nabla }\cdot \mathit{E}={\displaystyle \frac{{\mathit{\rho }}_q}{{\mathit{\epsilon }}_{\mathrm{0}}}}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\times \mathit{B}={\mathit{\mu }}_{\mathrm{0}}\left(\mathit{J}+{\mathit{\epsilon }}_{\mathrm{0}}{\displaystyle \frac{\partial \mathit{E}}{\partial t}}\right)\end{array}$$

The generalized Ohm's law providing closure for the Vlasov system is

$$\begin{array}{}\text{(5)}& \mathit{E}+\mathit{V}\times \mathit{B}={\displaystyle \frac{\mathit{J}}{\mathit{\sigma }}}+{\displaystyle \frac{\mathit{J}\times \mathit{B}}{n_{\mathrm{e}}e}}-{\displaystyle \frac{\mathrm{\nabla }\cdot {\mathcal{P}}_{\mathrm{e}}}{n_{\mathrm{e}}e}}+{\displaystyle \frac{m_{\mathrm{e}}}{n_{\mathrm{e}}{e}^{\mathrm{2}}}}{\displaystyle \frac{\partial \mathit{J}}{\partial t}},\end{array}$$

where V is the plasma bulk velocity, σ is the conductivity, e is the elementary charge, n_e is the electron number density, and 𝒫_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, $\mathit{J}\times \mathit{B}/\left(n_{\mathrm{e}}e\right)$, and the electron pressure gradient term, $\mathrm{\nabla }\cdot {\mathcal{P}}_{\mathrm{e}}/\left(n_{\mathrm{e}}e\right)$. In hybrid models, a true description of electron pressure is unavailable 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:

$$\begin{array}{}\text{(6)}& \mathit{E}+\mathit{V}\times \mathit{B}={\displaystyle \frac{\mathrm{1}}{e{n}_{\mathrm{p}}{\mathit{\mu }}_{\mathrm{0}}}}(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}.\end{array}$$

As Vlasov methods do not propagate particles but rather evolve distribution functions, we now briefly explain the semi-Lagrangian method employed by Vlasiator (for a full description, see chapter 5.3.1 in 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 (also referred to as translation, the second term of Vlasov's Eq. 1) and acceleration (the third term of Vlasov's Eq. 1) calculated one after the other with a leapfrog offset of $\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }t$. In this paper, Δ denotes steps on the full simulation grid and associated time step and δ is used to indicate calculations performed as substepping. 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 value, which is usually set to 22^∘. This value is constrained by the SLICE-3D shear approach, with values much above 22^∘ resulting in unphysical deformation of the VDF and smaller values increasing the computational cost of the simulation. For each acceleration step of length Δt, a transformation matrix is initialized as an identity matrix. The transformation matrix is first composed to apply the uniform electric field acceleration and the gyromotion due to the magnetic Lorentz force. Then, the transformation matrix is decomposed into three shear transformations. For a detailed explanation of the approach see chapter 3.5.1 of Palmroth et al. (2018). The transformation matrix is incrementally built with substepping of δt where each δt corresponds to a 0.1^∘ Larmor gyration, with the gyration step derived from convergence tests. Instead of applying linear acceleration by the motional electric field, 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 the electric field vanishes, 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 an electric field is found via the MHD Ohm's law with the Hall term included (Eq. 6). This Hall frame estimates the frame of reference of electrons, assuming electrons generate a current density which corresponds to the local magnetic field 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.

3 The eVlasiator global electron solver

In this section we introduce a novel method of simulating electron dynamics within the Earth's magnetic domain by building on the strengths of Vlasiator simulations. The method, called eVlasiator, focuses on the evolution of accurately modelled velocity distribution functions based on global plasma dynamics and structures evolved by the hybrid model. The spatial scales used in Vlasiator are not sufficient to resolve in detail small-scale phenomena such as electron-dominated reconnection, but this balances out with a realistic representation of global structures and asymmetries of the whole magnetosphere. The eVlasiator model solves the Vlasov equation for electron distribution functions using mostly the same methodology as Vlasiator itself but applies a simplified field solver, neglecting magnetic field evolution.

3.1 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. 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 Londrillo and Del Zanna, 2004). Both protons and electrons for the eVlasiator simulation are initialized from the read moments as Maxwellian distribution functions, with electron bulk velocity selected so that Ampère's law (Eq. 4) is fulfilled. Re-mapping input-run Vlasiator proton VDFs as Maxwellians does not affect the simulation results as eVlasiator only considers the proton number density and bulk velocity for current density calculations and does not propagate the proton distribution functions, instead keeping their characteristics completely constant for the duration of the eVlasiator simulation. For each simulation cell, we use the Balsara (2009) approach for calculating cell-average volumetric magnetic fields and respective derivatives. The eVlasiator solver uses volumetric field derivatives for calculating ∇×B.

3.2 The eVlasiator electron solver

eVlasiator solves the evolution of electron eVDFs similar to how Vlasiator simulates proton VDFs (for a detailed explanation, see Palmroth et al., 2018, in particular chapter 5.3.1). Solving the Vlasov Eq. (1) is split into two sections, translation and acceleration, with each of these steps performed in a staggered leapfrog approach. This approach is described in Fig. 1 with the first row indicating the spatial advection of electrons (translation) and the second row describing the effect of the Lorentz force on electrons through electric field acceleration and gyromotion. At time t (or t₀ at the initial state) we have the 5D (2D-3V) or 6D (3D-3V) electron velocity distributions (and, by extension, moments) in the whole simulation domain as well as proton moments and volumetric magnetic fields. Proton and magnetic field data are as read from the Vlasiator simulation and kept constant throughout the eVlasiator simulation. At simulation start, the leapfrog stepping is initialized with a half-length acceleration step (shown in red as step 0. in Fig. 1).

During each translation step, as depicted in Fig. 1 and described by the equation

$$\begin{array}{}\text{(7)}& {\left.{\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial t}}\right|}_{\mathrm{trans}}+\mathit{v}\cdot {\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial \mathit{x}}}=\mathrm{0},\end{array}$$

we perform a semi-Lagrangian spatial advection operation using second-order polynomial remapping, in an identical fashion to in a regular Vlasiator. This is evaluated separately for each cell in the gridded electron distribution functions using the velocity for that cell and is evaluated for one Cartesian direction at a time.

During each acceleration step, as depicted in Fig. 1 and described by the equation

$$\begin{array}{}\text{(8)}& {\left.{\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial t}}\right|}_{\mathrm{acc}}+{\displaystyle \frac{q_{\mathrm{s}}}{m_{\mathrm{s}}}}\left(\mathit{E}+\mathit{v}\times \mathit{B}\right)\cdot {\displaystyle \frac{\partial f_{\mathrm{s}}}{\partial \mathit{v}}}=\mathrm{0},\end{array}$$

we perform a semi-Lagrangian velocity space SLICE-3D update of the whole local distribution function, separately for each spatial cell. This method evaluates the acceleration due to electric fields (uniform movement in velocity space) and the rotation due to the Lorentz force magnetic component (a rotation in velocity space). The uniform movement and the rotation are composed into a transformation matrix. To apply the transformation with the SLICE-3D scheme, the matrix is then decomposed into three shear motions, one along each Cartesian velocity coordinate axis, and performed using semi-Lagrangian fourth-order polynomial remapping, similar to how the regular Vlasiator Vlasov solver works. This approach is applicable as long as velocities are non-relativistic. For a detailed description, see chapter 5.3.1 of Palmroth et al. (2018).

Due to the inherent connection between rapid electron motion and the local electric field response, we update electric fields in tandem with electron acceleration. The approach is detailed in the next subsection.

3.3 The eVlasiator field solver

In the eVlasiator field solver we maintain static magnetic fields as read from the input Vlasiator simulation, only calculating electric field evolution. We model the electric field by including additional terms in Ohm's law (Eq. 5), allowing for the interaction of electron distribution functions with electron–oscillation electric fields. Whistler mode propagation is not included in this study. We do not include any electric field due to Gauss's 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 (Eq. 2).

• Collisionless plasma physics assumes that electrons are fast enough to balance out any charge imbalance, and in hybrid-kinetic simulations this holds true. We do not implement Gauss's law (Eq. 3) in order to quantify to what extent charge neutrality holds in the eVlasiator context.

• The last term in Ampère's law (Eq. 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.

• As our plasma remains collisionless, we maintain our assumption of infinite conductivity, and thus the J∕σ term in the generalized Ohm's law (Eq. 5) remains zero.

• The Hall term, $\mathit{J}\times \mathit{B}/\left(n_{\mathrm{e}}e\right)$, is used to estimate the electron reference frame. This term is no longer required, as the Lorentz gyromotion rotation can be performed in the actual electron bulk motion reference frame.

• As eVlasiator models electrons with full distribution functions, we include the full electron pressure tensor 𝒫_e and implement the electron pressure gradient term using spatial gradients calculated for electron pressure.

• 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 choose to limit the acceleration time step Δt 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 both oscillations as a result of convergence tests. Much larger values will result in neighbouring simulation cells with different plasma characteristics diverging into an unstable state, and much lower values will needlessly cause an increase in computational cost. Due to the computational cost of SLICE-3D remapping, a substepping approach is used in order to more accurately model the electron gyromotion and plasma oscillation. Whilst the 22^∘ step models eVDF evolution to a high accuracy, the accurate and stable simulation of feedback between electron velocity, plasma oscillation, and the electric field due to the electron inertia term in Ohm's law requires substepping and places strict requirements on the length of the substep δt. This substepping is performed in tandem with the SLICE-3D transformation matrix generation. The electron gyroperiod is ${\mathit{\tau }}_{\mathrm{ce}}=\mathrm{2}\mathit{\pi }{\mathit{\omega }}_{\mathrm{ce}}^{-\mathrm{1}}$ and the plasma oscillation time is ${\mathit{\tau }}_{\mathrm{pe}}=\mathrm{2}\mathit{\pi }{\mathit{\omega }}_{\mathrm{pe}}^{-\mathrm{1}}$, where the electron plasma frequency is

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{pe}}=\sqrt{{\displaystyle \frac{n_{\mathrm{e}}{e}^{\mathrm{2}}}{{\mathit{\epsilon }}_{\mathrm{0}}{m}_{\mathrm{e}}}}}\end{array}$$

and the electron gyrofrequency is

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{ce}}={\displaystyle \frac{eB}{m_{\mathrm{e}}}}.\end{array}$$

In transformation matrix generation, substepping is constrained to a maximum of $\mathit{\delta }t\le min({\mathit{\tau }}_{\mathrm{pe}},{\mathit{\tau }}_{\mathrm{ce}})/\mathrm{3600}$. This value was defined as a result of convergence tests, and its dependence on the relationship between τ_pe and τ_ce is discussed more in Sect. 4.

The electron oscillation and electric field feedback loop is handled in parallel with gyration by tracking a cell-volume-averaged electric field ${\mathit{E}}_{J_{\mathrm{e}}}$ which is itself derived from the small-scale electron oscillation. For each acceleration substep, we update electron motion V and the electric field ${\mathit{E}}_{J_{\mathrm{e}}}$ by performing two parallel fourth-order Runge–Kutta propagations. The RK4 algorithm was chosen instead of a Runge–Kutta–Nyström method as it provides a good balance between general applicability, stability, and computational performance. The first one is

$$\begin{array}{}\text{(11)}& \mathit{\delta }{\mathit{V}}_{\mathrm{e}}=\mathit{\delta }t{\displaystyle \frac{e}{m_{\mathrm{e}}}}{\mathit{E}}_{J_{\mathrm{e}}},\end{array}$$

tracking electron bulk velocity response δV_e to the ${\mathit{E}}_{J_{\mathrm{e}}}$ field. This simple acceleration term is in fact equal to evaluating current changes via the electron inertia term in Ohm's law with the ${\mathit{E}}_{J_{\mathrm{e}}}$ field included in the left-hand-side electric field. The second Runge–Kutta propagation tracks the evolution of the ${\mathit{E}}_{J_{\mathrm{e}}}$ field due to changing current density, according to the displacement current on the right-hand side of Ampère's law (Eq. 4) with the ∇×B term in Ampère's law fixed to the static input magnetic fields. Thus, for each Runge–Kutta step, the electric field ${\mathit{E}}_{J_{\mathrm{e}}}$ is updated with

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathit{\delta }{\mathit{E}}_{J_{\mathrm{e}}}& {\displaystyle }=\mathit{\delta }t\left({\displaystyle \frac{\mathrm{\nabla }\times \mathit{B}}{{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}}-{\displaystyle \frac{\mathit{J}}{{\mathit{\epsilon }}_{\mathrm{0}}}}\right)\text{(13)}& {\displaystyle }& {\displaystyle }=\mathit{\delta }t\left({\displaystyle \frac{\mathrm{\nabla }\times \mathit{B}}{{\mathit{\epsilon }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{0}}}}-{\displaystyle \frac{e{\mathit{V}}_{\mathrm{p}}{n}_{\mathrm{p}}-e{\mathit{V}}_{\mathrm{e}}{n}_{\mathrm{e}}}{{\mathit{\epsilon }}_{\mathrm{0}}}}\right)\text{(14)}& {\displaystyle }& {\displaystyle }=\mathit{\delta }t{c}^{\mathrm{2}}\left(\mathrm{\nabla }\times \mathit{B}+{\mathit{\mu }}_{\mathrm{0}}e\left(n_{\mathrm{e}}{\mathit{V}}_{\mathrm{e}}-n_{\mathrm{p}}{\mathit{V}}_{\mathrm{p}}\right)\right),\end{array}$$

where c is the speed of light, and B, n_p and the proton bulk velocity V_p are assumed to be constant throughout the substep. Each of the four δV_e Runge–Kutta coefficients are updated with the latest estimate for $\mathit{\delta }{\mathit{E}}_{J_{\mathrm{e}}}$, and vice versa. Values for ${\mathit{E}}_{J_{\mathrm{e}}}$ 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 substepping procedure is visualized in the third row of Fig. 1. Further details of the solver and advection methods in Vlasiator can be found in Palmroth et al. (2018).

Figure 1 Electron solver procedure including substepping. At simulation start, a half-length acceleration step (0.) is performed. After that, translation (1, 3, … ) and acceleration (2, 4, … ) steps alternate in a leapfrog approach. Each acceleration step applies a transformation matrix which is generated in substeps, each of which updates electron acceleration δV_e and electric field change $\mathit{\delta }{\mathit{E}}_{J_{\mathrm{e}}}$. Each of these updates is performed via a dual Runge–Kutta 4 algorithm over step lengths δt with Runge–Kutta coefficients $k_{\mathrm{1}\mathrm{\dots }\mathrm{4}}^E$ and $k_{\mathrm{1}\mathrm{\dots }\mathrm{4}}^V$.

With each substep, the transformation matrix is evolved by compounding the following transformations:

1. Apply the acceleration δV_e derived from RK4-substepped ${\mathit{E}}_{J_{\mathrm{e}}}$ acceleration.

2. Accelerate electrons by $\frac{\mathrm{1}}{\mathrm{2}}\mathit{\delta }t{\mathit{E}}_{\mathrm{\nabla }{\mathcal{P}}_{\mathrm{e}}}$.

3. Transform to the frame of reference of the electron bulk motion.

4. Rotate the eVDF around the magnetic field direction for δtω_ce.

5. Transform back from the frame of reference of the electron bulk motion to the simulation frame.

6. Accelerate electrons by $\frac{\mathrm{1}}{\mathrm{2}}\mathit{\delta }t{\mathit{E}}_{\mathrm{\nabla }{\mathcal{P}}_{\mathrm{e}}}$.

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 eVDF.

3.4 Sample simulation set-up

In this method introduction, we use a noon–midnight meridional-plane 2D-3V Vlasiator simulation as our test-case input data. This 2D-3V Vlasiator simulation has been used to investigate global and kinetic magnetospheric dynamics in multiple studies such as Palmroth et al. (2017), Hoilijoki et al. (2017), Juusola et al. (2018a), Juusola et al. (2018b), Hoilijoki et al. (2019b), Grandin et al. (2019), and Akhavan-Tafti et al. (2020). It has solar wind values of β=0.7, magnetosonic Mach number M_ms=5.6, Alfvén Mach number M_A=6.9, proton number density $n_{\mathrm{p}}=\mathrm{1}{\mathrm{cm}}^{-\mathrm{3}}$, and solar wind speed u_sw along the ${\widehat{e}}_x$ (Earth–Sun) direction with $u_{\mathrm{sw},x}=-\mathrm{750}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$, simulating fast solar wind conditions and ensuring efficient simulation initialization. The simulation input interplanetary magnetic field is purely southward with $B_z=-\mathrm{5}\mathrm{nT}$ and the Earth's magnetic dipole is a ${\widehat{e}}_z$-aligned line dipole scaled to result in a realistic magnetopause standoff distance (Daldorff et al., 2014). The simulation has an inner boundary at $\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}\approx \mathrm{4.7}$ Earth radii, modelled as a perfectly conducting sphere. The spatial resolution is Δx=300 km.

For this eVlasiator sample run, we choose a region from the magnetotail with $\mathrm{70}\times \mathrm{1}\times \mathrm{40}$ simulation cells in the X, Y, and Z directions, respectively. The subregion extent is from $x_{-}=-\mathrm{75.6}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$ to $x_{+}=-\mathrm{54.6}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$, from $y_{-}=-\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$ to $y_{+}=+\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$, and from $z_{-}=-\mathrm{6}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$ to $z_{+}=+\mathrm{6}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$. Within this domain, visualized with a small rectangle in Fig. 2a, the electron plasma period τ_pe ranges from ∼0.7 ms in the magnetotail plasma sheet up to ∼2.5 ms in the near-plasmasphere lobes. The electron gyroperiod τ_ce ranges from ∼14 ms in most of the lobes up to ∼770 ms at a tail current sheet X-line site.

Figure 2 Simulation box initialization values. (a) Close-up of the central 16 % section of the Vlasiator input simulation with plasma number density overlaid with magnetic field lines. A small rectangle in the magnetotail region indicates the electron simulation domain (b–f). (b) Proton number density overlaid with magnetic field lines. X-line topology is visible at $X\sim -\mathrm{73}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$, $Z\sim -\mathrm{0.5}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$. (c) Proton temperature as a scalar. Electron initialization temperatures are scaled down by a constant factor 4. (d) Ratio of electron plasma and gyrofrequencies. (e, f) Proton and electron bulk velocity magnitudes with in-plane directions indicated with vectors.

The electron distributions are discretized onto eVlasiator velocity meshes, with the electron velocity mesh consisting of 400³ cells, extending from $-\mathrm{4.2}\times {\mathrm{10}}^{\mathrm{7}}$ to $+\mathrm{4.2}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ in each direction, resulting in an electron velocity space resolution of 210 km s⁻¹. The eVDF sparsity threshold was set to ${\mathrm{10}}^{-\mathrm{21}}{\mathrm{m}}^{-\mathrm{6}}{\mathrm{s}}^{\mathrm{3}}$, ensuring good representation of the main structure of the eVDF. Discretizing a hot and dense electron distribution onto a Cartesian grid is numerically challenging without using vast amounts of memory. As portions of our simulation domain have proton temperature up to 10⁸ K, we use an empirical estimate of $T_{\mathrm{i}}/T_{\mathrm{e}}\sim \mathrm{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_{\mathrm{i}}/m_{\mathrm{e}}=\mathrm{183.6}$. As mentioned 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 (Eq. 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 Fig. 2 along with an overview of the input Vlasiator simulation and the selected electron subdomain.

4 Solver performance

4.1 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 ${\mathit{E}}_{J_{\mathrm{e}}}$ shown in Fig. 3. These single-cell tests did not have magnetic field curvature or an ion population present, resulting in the electron motion oscillating around a stability point at V_e=0 and ${\mathit{E}}_{J_{\mathrm{e}}}=\mathrm{0}$. We set the electron number density to $n_{\mathrm{e}}=\mathrm{0.1}{\mathrm{cm}}^{-\mathrm{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 ${\mathit{V}}_{\mathrm{e},\mathrm{0}}=(-\mathrm{100},-\mathrm{150},\mathrm{200})\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$, close to but below our electron velocity resolution of $\mathrm{\Delta }v=\mathrm{210}\mathrm{km}{\mathrm{s}}^{-\mathrm{1}}$. As can be seen from Fig. 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 the plasma and gyroperiods to values closer to each other (1.11 and 1.79 ms, respectively), we see a gradual evolution of oscillation amplitude and, thus, ${\mathit{E}}_{J_{\mathrm{e}}}$ 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. Also, this instability occurs only when τ_ce≈τ_pe, which does not occur in our full simulation domain.

Figure 3 Graphs of solver stability in relation to electron plasma oscillation and gyromotion in a single-cell simulation. Note the different time axes used. (a, c, e) Oscillation electric field ${\mathit{E}}_{J_{\mathrm{e}}}$ components. (b, d, f) Electron bulk velocity V_e components. (a, b) Graph values in relation to the electron plasma oscillation period (indicated with a thick grey bar) and (c, d) in relation to the electron gyroperiod (indicated with a thick black bar), with a background magnetic field of B=20 nT. (e, f) A simulation with a magnetic field of B=200 nT, resulting in the gyromotions and oscillatory motions interacting over multiple periods.

4.2 Dispersion relation analysis

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 the kinetic simulation approach, a correct reproduction of wave dispersion behaviour is a good indicator of correct physical behaviour of the simulation system.

Two 1D-simulation set-ups with a spatial grid resolution of Δx=300 m (= 0.01 d_e) and N_x=1000 cells were initialized with an electron number density of $n_{\mathrm{e}}=\mathrm{0.4}\times {\mathrm{10}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{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 $\stackrel{\mathrm{̃}}{v}=\mathrm{1000}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The simulation was run for 0.037 s (433 ${\mathit{\omega }}_{\mathrm{pe}}^{-\mathrm{1}}$).

Figure 4 shows the dispersion data resulting from a 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 (Fig. 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.

Figure 4 Dispersion analysis of the electron solver in a 1D test case with an axis-parallel (a) and axis-perpendicular (b) magnetic field. The colour map shows the spatio-temporal Fourier transform of $E_{J_{\mathrm{e}},\parallel }$ (a) and $E_{J_{\mathrm{e}},⟂}$ (b) overlaid with analytical solutions for the Langmuir wave (black dashed curve) and Bernstein modes (black solid curves).

4.3 Stability within larger simulation domain

We also evaluate the stability of our solver over the larger simulated domain described in Sect. 3.4, with initialization values derived from the Vlasiator hybrid-Vlasov simulation. These graphs are shown in Fig. 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 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 the majority of the simulation domain electrons remain gyrotropic and even peak values do not grow past 10⁻³. Panel (g) shows statistics for the electron number density deviation from the initialization value, indicating loss of plasma neutrality due to the motion of electrons. The minimum value oscillating between approximately 10⁻⁹ and ${\mathrm{10}}^{-\mathrm{6}}{\mathrm{cm}}^{-\mathrm{3}}$ indicates the level of numerical fluctuations, and the maximum, mean and median values show how charge imbalance does grow initially but stabilizes within about 0.1 s and does not grow beyond ${\mathrm{10}}^{-\mathrm{1}}{\mathrm{cm}}^{-\mathrm{3}}$.

In panels (h) through (k) of Fig. 5 we show how the instantaneous plasma oscillation electric field ${\mathit{E}}_{J_{\mathrm{e}}}$ is well-behaved throughout the simulation box, converging towards stable values. We note that as the ${\mathit{E}}_{J_{\mathrm{e}}}$ field oscillates 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 an arbitrary phase of the oscillation. In panel (l) we show the normalized current density J departure from the balance current ${\mathit{J}}_{\mathbf{B}}=\frac{\mathrm{\nabla }\times \mathit{B}}{{\mathit{\mu }}_{\mathrm{0}}}$ which would be required to maintain the magnetic field structure according to Ampère's law (Eq. 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 $-\frac{\mathrm{\nabla }\cdot {\mathcal{P}}_{\mathrm{e}}}{n_{\mathrm{e}}e}$ 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 gives identically zero field components at local extrema of pressure.

Figure 5 Evolution of electron and solver parameters over the whole simulation domain. (a–d) Minimum, maximum, mean, and median values for electron temperature T_e 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 Q_Ag,e. (g) Minimum, maximum, mean, and median values for electron density deviation from initial state, indicating charge imbalance. (h–k) Minimum, maximum, mean, and median values for the plasma oscillation electric field ${\mathit{E}}_{J_{\mathrm{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 ${\mathit{J}}_{\mathbf{B}}=\frac{\mathrm{\nabla }\times \mathit{B}}{{\mathit{\mu }}_{\mathrm{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.

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, ${\mathit{E}}_{J_{\mathrm{e}}}$, as could be expected (not shown).

5 Results

Results from the electron simulation after 1.0 s of evolution are presented in Fig. 6. Figure 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 plasma sheet boundary layer (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 has become non-gyrotropic. In most of the simulation domain, the value is very small, but enhanced agyrotropy (still relatively small values below 10⁻³) 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 6 Electron distribution properties within the test domain after 1.0 s of simulation. (a) The ratio of parallel electron temperature at 1.0 s to the parallel temperature at the start of the simulation, indicating parallel heating. (b) The same but for perpendicular temperature. (c) The agyrotropy measure for the electron population. (d) The magnitude and direction of the electron pressure gradient term of the electric field. (e, f) The charge imbalance $n_{\mathrm{e}}-n_{\mathrm{e},\mathrm{0}}$ and relative charge imbalance $(n_{\mathrm{e}}-n_{\mathrm{e},\mathrm{0}})n_{\mathrm{e},\mathrm{0}}^{-\mathrm{1}}$ found at the end of the simulation.

Figure 6d shows the electric field due to ∇𝒫_e, with the field strongest where the PSBL meets the magnetosphere. The field 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 Fig. 6f as the change scaled by the original electron number density. In the majority of the simulation domain, imbalance remains below ${\mathrm{10}}^{-\mathrm{2}}{\mathrm{cm}}^{-\mathrm{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 Fig. 7 we display electron velocity distribution functions after 1.0 s of simulation. Figure 7a shows the evolved electron temperature anisotropy $T_{⟂,\mathrm{e}}{T}_{\parallel ,\mathrm{e}}^{-\mathrm{1}}$, and Fig. 7b displays the maximum of instantaneous values of ${\mathit{E}}_{J_{\mathrm{e}}}$, taken over 10 measurements at 0.05 s intervals near the end of the simulation. Panels (c) through (n) of Fig. 7 show parallel and perpendicular projections of electron eVDFs at virtual spacecraft (VSC) [1] through [6], with positions of VSC indicated in panels (a) and (b).

Figure 7 Electron properties and velocity distribution functions after 1.0 s of simulation. (a) Electron temperature anisotropy $T_{⟂,\mathrm{e}}{T}_{\parallel ,\mathrm{e}}^{-\mathrm{1}}$ overlaid with magnetic field lines and six virtual spacecraft locations, labelled [1]–[6]. (b) Maximum value for displacement current ${\mathit{E}}_{J_{\mathrm{e}}}$, taken over 10 measurements at 0.05 s intervals near the end of the simulation. (c–n) Electron velocity distribution function projections into the parallel v_B and v_V×B or perpendicular $v_{\mathbf{B}\times (\mathbf{B}\times \mathbf{V})}$ and v_V×B planes. Each virtual spacecraft is indicated by the number in the parallel eVDF panel with the panel below showing the corresponding perpendicular eVDF for the same virtual spacecraft.

Figure 7a shows how temperature anisotropy $T_{⟂,\mathrm{e}}{T}_{\parallel ,\mathrm{e}}^{-\mathrm{1}}$ indicates parallel energization in the low-density regions adjacent to the PSBL and perpendicular energization adjacent to the X line and within the tailmost region 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 perpendicularly heated electrons to the oscillation region and parallel accelerated electrons propagating along field lines to the near-magnetosphere PSBL regions.

The maxima of instantaneous values of ${\mathit{E}}_{J_{\mathrm{e}}}$, shown in Fig. 7b, indicate that the strongest electron oscillations on our 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 ${\mathit{E}}_{J_{\mathrm{e}}}$ is also seen 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 Fig. 7a and virtual spacecraft measurements indicate that parallel features, akin to electron beams, are indeed found in regions with enhanced ${\mathit{E}}_{J_{\mathrm{e}}}$.

The temperature anisotropies found in the near-Earth tail region of our simulation are mostly in the 0.5 … 1.5 range. Artemyev et al. (2014) reported on Cluster observations of electron temperature anisotropies ranging from 0.8…1.6 and centred around ∼1.1, in agreement with our results, though their observations were gathered between $-\mathrm{20}{R}_{\mathrm{E}}<x<-\mathrm{15}{R}_{\mathrm{E}}$ ($-\mathrm{127}\times {\mathrm{10}}^{\mathrm{6}}<x<\mathrm{96}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$). Regions where parallel temperatures dominate (anisotropy <1) are found in regions of cold plasma, as can be seen by comparing Figs. 2c and 7a. This does not preclude the possibility of parallel acceleration in regions of hot plasma but rather shows that the acceleration may not be strong enough to be discerned over the main hot eVDF.

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 (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 Fig. 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⁻¹ electron velocities, which is reasonably within the range of our eVDFs in Fig. 7. We note that our simulation produces a background B_y profile with ∇B_y in agreement with Fig. 4 of Asano et al. (2006) (not shown), on top of which the streaming electrons are observed. Onsager et al. (1991) describe a simple 2D Liouville model for the PSBL, as well as some ISEE-1 and ISEE-2 observations supporting their model. The formation mechanisms of eVDFs in Onsager et al. (1991) are listed as time-of-flight, energy conservation and magnetic moment conservation, which are included in our model, though we perform a more robust evaluation of plasma oscillation interplay with gyration. The eVDFs shown in their Fig. 4 agree with our VSCs [1], [2], [5], and [6], for example. We also note our VSC [3] displaying a disjoint parallel beam, matching the ISEE-2 observations in Fig. 5 of Onsager et al. (1991).

Observations of perpendicular crescents are shown in MMS data in Burch et al. (e.g. 2016b, 2019) at EDRs, in conjunction with dayside magnetopause reconnection sites. These observed structures are produced at very small spatial scales, not captured by our current model. We do, however, observe similar agyrotropic crescents in our results further out (in particular in Fig. 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 Fig. 6e. Something akin to a parallel electron crescent (Burch et al., 2016b) can be seen in Fig. 7c, and bi-directional distributions as reported in Figs. 6 and 7 of Burch and Phan (2016) are qualitatively similar to our Fig. 7k and m.

6 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 ${\mathit{E}}_{J_{\mathrm{e}}}$ 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 eVDFs, 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 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 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 investigations, e.g. in the noon–midnight meridional plane. Electron velocity distribution functions generated by the model can be used to investigate, for example, 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 earthward. Li et al. (2020) observe electron Bernstein modes driven by perpendicular crescent distributions. As we have shown in Figs. 4 and 7, with sufficient resolution we can reproduce electron Bernstein waves and agyrotropic electron distributions. Thus, we are in a 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 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 indicate such problems. If such imbalances arise from a future model, some method of solving Gauss's law such as a Poisson solver should be implemented. A more detailed investigation into comparing electron eVDFs and dynamics with observations is expected in a future study.