the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Vlasov simulation of electrons in the context of hybrid global models: an eVlasiator approach
Markus Battarbee
Thiago Brito
Markku Alho
Yann PfauKempf
Maxime Grandin
Urs Ganse
Konstantinos Papadakis
Andreas Johlander
Lucile Turc
Maxime Dubart
Minna Palmroth
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 hybridkinetic simulations, which approximate electron dynamics as a chargeneutralizing fluid. We introduce a novel method for Vlasovsimulating electrons in the context of a hybridkinetic framework in order to examine the energization processes of magnetospheric electrons. Our extension of the Vlasiator hybridVlasov 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 singlecell and smallscale domains, and the solver successfully models Langmuir waves and Bernstein modes. We simulate a small testcase section of the nearEarth magnetotail plasma sheet region, reproducing a number of electron distribution function features found in spacecraft measurements.
Physical processes in nearEarth space are dominated by plasma effects such as nonthermal 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 pointlike 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 nearEarth 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 hybridkinetic models, where ions are treated as a kinetic selfconsistent species and electrons are a chargeneutralizing fluid. Successful approaches include hybridVlasov models (Palmroth et al., 2018) and hybridPIC (particleincell) 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 fullPIC (ions and electrons as interacting particles, e.g. Hesse et al., 2005), fullVlasov (ions and electrons as interacting distribution functions, e.g. Umeda et al., 2009; Schmitz and Grauer, 2006; Pezzi et al., 2019), hybridPIC electrons (dynamic electron particles, ions as a static background, e.g. Lapenta et al., 2007) and hybridVlasov 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 iontoelectron 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 selfheating 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 semiimplicit 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. FullPIC 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 smallscale threedimensional 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, MHDEPIC, is presented by Daldorff et al. (2014), with a twoway 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. MHDEPIC 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 smallscale 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 sixmoment multifluid fullMaxwell model. They note that they do not capture reconnection to an acceptable accuracy and have yet to publish global simulation results. Global 10moment 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 electronscale 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. Reconnectiondriven 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 electrontoproton 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 electrondriven waves and timedomain 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 iondetermined fields. The aim is to investigate how much of the global electron physics and distribution functions can be understood by utilizing iongenerated field as modelled by hybridkinetic codes, as opposed to a numerically unfeasible global fullkinetic approach. The paper is organized as follows. In Sect. 2, we introduce the ionkinetic hybridVlasov 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 ionkinetic 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.
Vlasiator (von Alfthan et al., 2014; Palmroth et al., 2018) is, at the present time, the only hybridVlasov code capable of simulating the global magnetospheric system of the Earth, accounting for ionkinetic 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 threedimensional 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 velocityspace grid has a resolution of 30 km s^{−1}, extending between $\pm \mathrm{4020}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$. To constrain computational cost and memory usage, those parts of the velocity distribution function which have a phasespace 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^{−17} and ${\mathrm{10}}^{\mathrm{15}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\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. Largescale global 3D runs will be made possible in the near future by adaptive mesh refinement (AMR), using nonuniform 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 selfconsistent simulation. However, until now, the electron population has been treated in the usual way of implementing it as a massless chargeneutralizing 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 selfconsistency enforced through plasma oscillation and current densities.
2.1 Solving the Vlasov equation
Vlasiator uses the hybridVlasov ion approach to model the nearEarth space plasma environment. The full sixdimensional (6D) phase space density ${f}_{\mathrm{s}}(\mathit{x},\phantom{\rule{0.125em}{0ex}}\mathit{v},\phantom{\rule{0.125em}{0ex}}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 μ_{0} and ε_{0} 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 divergencefree magnetic field initialcondition 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 staggeredgrid 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
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 righthand side. In the limit of slow temporal variations, the electron inertia term (the last term on the righthand side) also vanishes. The remaining two terms on the righthand 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 iontoelectron temperature ratio, or by neglecting the small electron pressure gradient term altogether. The standard ionhybrid 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 ionhybrid 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 semiLagrangian 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 SLICE3D 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 SLICE3D 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 Borispush 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 SLICE3D algorithm.
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 smallscale phenomena such as electrondominated 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 largescale fields and moments of a Vlasiator simulation of nearEarth space. In the eVlasiator approach, we read magnetic field vectors and proton plasma moments for the chosen simulation domain and apply userdefined 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 ionhybrid simulation proton moments, cellfaceaverage magnetic field components and celledgeaverage electric field components (the latter being used by the staggeredgrid 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. Remapping inputrun 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 cellaverage 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_{0} at the initial state) we have the 5D (2D3V) or 6D (3D3V) 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 halflength 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
we perform a semiLagrangian spatial advection operation using secondorder 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
we perform a semiLagrangian velocity space SLICE3D 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 SLICE3D scheme, the matrix is then decomposed into three shear motions, one along each Cartesian velocity coordinate axis, and performed using semiLagrangian fourthorder polynomial remapping, similar to how the regular Vlasiator Vlasov solver works. This approach is applicable as long as velocities are nonrelativistic. 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 hybridkinetic 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 SLICE3D 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 SLICE3D 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 SLICE3D 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
and the electron gyrofrequency is
In transformation matrix generation, substepping is constrained to a maximum of $\mathit{\delta}t\le min({\mathit{\tau}}_{\mathrm{pe}},\phantom{\rule{0.125em}{0ex}}{\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 cellvolumeaveraged electric field ${\mathit{E}}_{{J}_{\mathrm{e}}}$ which is itself derived from the smallscale electron oscillation. For each acceleration substep, we update electron motion V and the electric field ${\mathit{E}}_{{J}_{\mathrm{e}}}$ by performing two parallel fourthorder 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
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 lefthandside 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 righthand 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
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).
With each substep, the transformation matrix is evolved by compounding the following transformations:

Apply the acceleration δV_{e} derived from RK4substepped ${\mathit{E}}_{{J}_{\mathrm{e}}}$ acceleration.

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

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

Rotate the eVDF around the magnetic field direction for δtω_{ce}.

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

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 SLICE3D algorithm, which decomposes the transformation into three Cartesian shears and updates the eVDF.
3.4 Sample simulation setup
In this method introduction, we use a noon–midnight meridionalplane 2D3V Vlasiator simulation as our testcase input data. This 2D3V 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 AkhavanTafti 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{cm}}^{\mathrm{3}}$, and solar wind speed u_{sw} along the ${\widehat{e}}_{x}$ (Earth–Sun) direction with ${u}_{\mathrm{sw},\phantom{\rule{0.125em}{0ex}}x}=\mathrm{750}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.125em}{0ex}}\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}}\phantom{\rule{0.125em}{0ex}}\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}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ to ${x}_{+}=\mathrm{54.6}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$, from ${y}_{}=\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ to ${y}_{+}=+\mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$, and from ${z}_{}=\mathrm{6}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}$ to ${z}_{+}=+\mathrm{6}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\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 nearplasmasphere lobes. The electron gyroperiod τ_{ce} ranges from ∼14 ms in most of the lobes up to ∼770 ms at a tail current sheet Xline site.
The electron distributions are discretized onto eVlasiator velocity meshes, with the electron velocity mesh consisting of 400^{3} cells, extending from $\mathrm{4.2}\times {\mathrm{10}}^{\mathrm{7}}$ to $+\mathrm{4.2}\times {\mathrm{10}}^{\mathrm{7}}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ in each direction, resulting in an electron velocity space resolution of 210 km s^{−1}. The eVDF sparsity threshold was set to ${\mathrm{10}}^{\mathrm{21}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\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^{8} 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 iontoelectron 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 (cellaverage) 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.1 Singlecell stability of electron oscillation
To validate the performance of our electron solver, we performed singlecell tests, with resultant electron bulk velocities V_{e} and plasma oscillation electric fields ${\mathit{E}}_{{J}_{\mathrm{e}}}$ shown in Fig. 3. These singlecell 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}\phantom{\rule{0.125em}{0ex}}{\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},\phantom{\rule{0.125em}{0ex}}\mathrm{0}}=(\mathrm{100},\mathrm{150},\mathrm{200})\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$, close to but below our electron velocity resolution of $\mathrm{\Delta}v=\mathrm{210}\phantom{\rule{0.125em}{0ex}}\mathrm{km}\phantom{\rule{0.125em}{0ex}}{\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.
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 1Dsimulation setups 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}}\phantom{\rule{0.125em}{0ex}}{\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{\u0303}}{v}=\mathrm{1000}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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.
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 hybridVlasov 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, Bparallel, and Bperpendicular 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^{−3}. 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^{−9} and ${\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\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}}\phantom{\rule{0.125em}{0ex}}{\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 wellbehaved 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 maximumvalue 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.
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 substepassociated electron bulk frame. This transformation choice resulted in unstable growth of, in particular, ${\mathit{E}}_{{J}_{\mathrm{e}}}$, as could be expected (not shown).
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 endofsimulation temperature to initial temperatures. Heating is found in particular near the Xline 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 nongyrotropic. 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 ∇𝒫_{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}}\phantom{\rule{0.125em}{0ex}}{\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}_{\u27c2,\phantom{\rule{0.125em}{0ex}}\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 7a shows how temperature anisotropy ${T}_{\u27c2,\phantom{\rule{0.125em}{0ex}}\mathrm{e}}{T}_{\parallel ,\phantom{\rule{0.125em}{0ex}}\mathrm{e}}^{\mathrm{1}}$ indicates parallel energization in the lowdensity 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 highbeta 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 nearmagnetosphere 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 electrondriven waves in the PSBL (Onsager et al., 1993). Some increase in ${\mathit{E}}_{{J}_{\mathrm{e}}}$ is also seen at the Xline 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 nearEarth 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}\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{E}}<x<\mathrm{15}\phantom{\rule{0.125em}{0ex}}{R}_{\mathrm{E}}$ ($\mathrm{127}\times {\mathrm{10}}^{\mathrm{6}}<x<\mathrm{96}\times {\mathrm{10}}^{\mathrm{6}}\phantom{\rule{0.125em}{0ex}}\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 bidirectional 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 bidirectional 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 shiftedfootball 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^{−1} 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 ISEE1 and ISEE2 observations supporting their model. The formation mechanisms of eVDFs in Onsager et al. (1991) are listed as timeofflight, 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 ISEE2 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 phasespace 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 bidirectional distributions as reported in Figs. 6 and 7 of Burch and Phan (2016) are qualitatively similar to our Fig. 7k and m.
In this method paper we have presented a novel approach to investigating electron distribution function dynamics in the context of global ionhybrid field structures. Our method exploits global dynamics provided by hybridVlasov 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, bidirectional distributions, and perpendicular crescents.
Our model has several builtin limitations as it does not treat electrons as a fully selfconsistent 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 hybridVlasov simulation, it has a spatial resolution far below that required for resolving electronscale 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 longterm 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 nonuniform 3D3V Vlasiator meshes, allowing investigation of spatially threedimensional 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 selfconsistent 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.
Vlasiator (https://www.helsinki.fi/en/researchgroups/vlasiator, Palmroth, 2021) is distributed under the GPL2 opensource license at https://github.com/fmihpc/vlasiator/ (Palmroth and the Vlasiator team, 2020). Vlasiator uses a data structure developed inhouse (https://github.com/fmihpc/vlsv/, Sandroos, 2019). The Analysator software (https://doi.org/10.5281/zenodo.4462515; Battarbee and the Vlasiator team, 2020) was used to produce the presented figures. The run described here takes several gigabytes of disk space and is kept in storage maintained within the CSC – IT Center for Science. Data presented in this paper can be accessed by following the data policy on the Vlasiator web site.
MB wrote the paper and code description. MB and TB devised the solver method. MA assisted with data analysis, model development, and comparisons with observations. UG performed the dispersion tests. MP is the principal investigator of Vlasiator and leads the investigation. YPK, MG, KP, AJ, LT, MD, and MP participated in the discussion and finalization of the paper.
The authors declare that they have no conflict of interest.
We acknowledge the European Research Council for starting grant 200141QuESpace, with which the Vlasiator model (https://www.helsinki.fi/en/researchgroups/vlasiator, last access: 25 January 2021) was developed, and consolidator grant 682068PRESTISSIMO awarded for further development of Vlasiator and its use in scientific investigations. We gratefully acknowledge Academy of Finland grant numbers 309937TEMPO and 312351FORESAIL. PRACE (http://www.praceri.eu, last access: 25 January 2021) is acknowledged for granting us Tier0 computing time in HLRS Stuttgart, where Vlasiator was run in the HazelHen machine with project number 2014112573 and in the Hawk machine with project number 2019204998. The work of Lucile Turc is supported by the Academy of Finland (grant number 322544). The authors wish to thank the anonymous referees for their assistance in improving the approachability of the paper.
This research has been supported by the European Research Council (grant nos. 682068 and 200141) and the Academy of Finland, Luonnontieteiden ja Tekniikan Tutkimuksen Toimikunta (grant nos. 312351, 309937, and 322544).
Openaccess funding was provided by the Helsinki University Library.
This paper was edited by Wen Li and reviewed by two anonymous referees.
AkhavanTafti, M., Palmroth, M., Slavin, J. A., Battarbee, M., Ganse, U., Grandin, M., Le, G., Gershman, D. J., Eastwood, J. P., and Stawarz, J. E.: Comparative Analysis of the Vlasiator Simulations and MMS Observations of Multiple XLine Reconnection and Flux Transfer Events, J. Geophys. Res.Space, 125, e2019JA027410, https://doi.org/10.1029/2019JA027410, 2020. a
Artemyev, A. V., Baumjohann, W., Petrukovich, A. A., Nakamura, R., Dandouras, I., and Fazakerley, A.: Proton/electron temperature ratio in the magnetotail, Ann. Geophys., 29, 2253–2257, https://doi.org/10.5194/angeo2922532011, 2011. a
Artemyev, A. V., Petrukovich, A. A., Nakamura, R., and Zelenyi, L. M.: Profiles of electron temperature and B_{z} along Earth's magnetotail, Ann. Geophys., 31, 1109–1114, https://doi.org/10.5194/angeo3111092013, 2013. a
Artemyev, A. V., Walsh, A. P., Petrukovich, A. A., Baumjohann, W., Nakamura, R., and Fazakerley, A. N.: Electron pitch angle/energy distribution in the magnetotail, J. Geophys. Res.Space, 119, 7214–7227, https://doi.org/10.1002/2014JA020350, 2014. a
Artemyev, A. V., Angelopoulos, V., Liu, J., and Runov, A.: Electron currents supporting the nearEarth magnetotail during current sheet thinning, Geophys. Res. Lett., 44, 5–11, https://doi.org/10.1002/2016GL072011, 2017. a
Asano, Y., Nakamura, R., Runov, A., Baumjohann, W., McIlwain, C., Paschmann, G., Quinn, J., Alexeev, I., Dewhurst, J. P., Owen, C. J., Fazakerley, A. N., Balogh, A., Rème, H., and Klecker, B.: Detailed analysis of lowenergy electron streaming in the nearEarth neutral line region during a substorm, Adv. Space Res., 37, 1382–1387, https://doi.org/10.1016/j.asr.2005.05.059, 2006. a, b
Balsara, D. S.: Divergencefree reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics, J. Comput. Phys., 228, 5040–5056, https://doi.org/10.1016/j.jcp.2009.03.038, 2009. a, b
Battarbee, M. and the Vlasiator team: Analysator: python analysis toolkit, Zenodo, https://doi.org/10.5281/zenodo.4462515, 2020. a
Bessho, N., Chen, L.J. J., Shuster, J. R., and Wang, S.: Electron distribution functions in the electron diffusion region of magnetic reconnection: Physics behind the fine structures, Geophys. Res. Lett., 41, 8688–8695, https://doi.org/10.1002/2014GL062034, 2014. a
Bessho, N., Chen, L.J. J., and Hesse, M.: Electron distribution functions in the diffusion region of asymmetric magnetic reconnection, Geophys. Res. Lett., 43, 1828–1836, https://doi.org/10.1002/2016GL067886, 2016. a
Birdsall, C. K. and Langdon, A. B.: Plasma physics via computer simulation, Taylor and Francis, New York, 2005. a
Boris, J. P.: Relativistic plasma simulationoptimization of a hybrid code, Proceedings of Fourth Conference on Numerical Simulations of Plasmas, Naval Research Laboratory, Washington D.C., USA, 2–3 November 1970. a
Breuillard, H., Le Contel, O., Retino, A., Chasapis, A., Chust, T., Mirioni, L., Graham, D. B., Wilder, F. D., Cohen, I., Vaivads, A., Khotyaintsev, Y. V., Lindqvist, P.A., Marklund, G. T., Burch, J. L., Torbert, R. B., Ergun, R. E., Goodrich, K. A., Macri, J., Needell, J., Chutter, M., Rau, D., Dors, I., Russell, C. T., Magnes, W., Strangeway, R. J., Bromund, K. R., Plaschke, F., Fischer, D., Leinweber, H. K., Anderson, B. J., Le, G., Slavin, J. A., Kepko, E. L., Baumjohann, W., Mauk, B., Fuselier, S. A., and Nakamura, R.: Multispacecraft analysis of dipolarization fronts and associated whistler wave emissions using MMS data, Geophys. Res. Lett., 43, 7279–7286, https://doi.org/10.1002/2016GL069188, 2016. a
Burch, J. L. and Phan, T. D.: Magnetic reconnection at the dayside magnetopause: Advances with MMS, Geophys. Res. Lett., 43, 8327–8338, https://doi.org/10.1002/2016GL069787, 2016. a, b
Burch, J. L., Moore, T. E., Torbert, R. B., and Giles, B. L.: Magnetospheric Multiscale Overview and Science Objectives, Space Sci. Rev., 199, 5–21, https://doi.org/10.1007/s1121401501649, 2016a. a, b
Burch, J. L., Torbert, R. B., Phan, T. D., Chen, L.J., Moore, T. E., Ergun, R. E., Eastwood, J. P., Gershman, D. J., Cassak, P. A., Argall, M. R., Wang, S., Hesse, M., Pollock, C. J., Giles, B. L., Nakamura, R., Mauk, B. H., Fuselier, S. A., Russell, C. T., Strangeway, R. J., Drake, J. F., Shay, M. A., Khotyaintsev, Y. V., Lindqvist, P.A., Marklund, G., Wilder, F. D., Young, D. T., Torkar, K., Goldstein, J., Dorelli, J. C., Avanov, L. A., Oka, M., Baker, D. N., Jaynes, A. N., Goodrich, K. A., Cohen, I. J., Turner, D. L., Fennell, J. F., Blake, J. B., Clemmons, J., Goldman, M., Newman, D., Petrinec, S. M., Trattner, K. J., Lavraud, B., Reiff, P. H., Baumjohann, W., Magnes, W., Steller, M., Lewis, W., Saito, Y., Coffey, V., and Chandler, M.: Electronscale measurements of magnetic reconnection in space, Science, 352, aaf2939, https://doi.org/10.1126/science.aaf2939, 2016b. a, b
Burch, J. L., Dokgo, K., Hwang, K. J., Torbert, R. B., Graham, D. B., Webster, J. M., Ergun, R. E., Giles, B. L., Allen, R. C., Chen, L. J., Wang, S., Genestreti, K. J., Russell, C. T., Strangeway, R. J., and Le Contel, O.: HighFrequency Wave Generation in Magnetotail Reconnection: Linear Dispersion Analysis, Geophys. Res. Lett., 46, 4089–4097, https://doi.org/10.1029/2019GL082471, 2019. a
Cattell, C., Dombeck, J., Wygant, J., Drake, J. F., Swisdak, M., Goldstein, M. L., Keith, W., Fazakerley, A., André, M., Lucek, E., and Balogh, A.: Cluster observations of electron holes in association with magnetotail reconnection and comparison to simulations, J. Geophys. Res.Space, 110, A01211, https://doi.org/10.1029/2004JA010519, 2005. a
Daldorff, L. K. S., Tóth, G., Gombosi, T. I., Lapenta, G., Amaya, J., Markidis, S., and Brackbill, J. U.: Twoway coupling of a global Hall magnetohydrodynamics model with a local implicit particleincell model, J. Comput. Phys., 268, 236–254, https://doi.org/10.1016/j.jcp.2014.03.009, 2014. a
Daldorff, L. K. S., Tóth, G., Gombosi, T. I., Lapenta, G., Amaya, J., Markidis, S., and Brackbill, J. U.: Twoway coupling of a global Hall magnetohydrodynamics model with a local implicit particleincell model, J. Comput. Phys., 268, 236–254, https://doi.org/10.1016/j.jcp.2014.03.009, 2014. a
Daughton, W., Roytershteyn, V., Karimabadi, H., Yin, L., Albright, B. J., Bergen, B., and Bowers, K. J.: Role of electron physics in the development of turbulent magnetic reconnection in collisionless plasmas, Nat. Phys., 7, 539–542, https://doi.org/10.1038/nphys1965, 2011. a
Deca, J., Divin, A., Lembège, B., Horányi, M., Markidis, S., and Lapenta, G.: General mechanism and dynamics of the solar wind interaction with lunar magnetic anomalies from 3D particleincell simulations, J. Geophys. Res.Space, 120, 6443–6463, https://doi.org/10.1002/2015JA021070, 2015. a
Deca, J., Divin, A., Henri, P., Eriksson, A., Markidis, S., Olshevsky, V., and Horányi, M.: Electron and Ion Dynamics of the Solar Wind Interaction with a Weakly Outgassing Comet, Phys. Rev. Lett., 118, 205 101, https://doi.org/10.1103/PhysRevLett.118.205101, 2017. a
Deca, J., Henri, P., Divin, A., Eriksson, A., Galand, M., Beth, A., Ostaszewski, K., and Horányi, M.: Building a Weakly Outgassing Comet from a Generalized Ohm's Law, Phys. Rev. Lett., 123, 055101, https://doi.org/10.1103/PhysRevLett.123.055101, 2019. a
Dong, C., Wang, L., Hakim, A., Bhattacharjee, A., Slavin, J. A., DiBraccio, G. A., and Germaschewski, K.: Global TenMoment Multifluid Simulations of the Solar Wind Interaction with Mercury: From the Planetary Conducting Core to the Dynamic Magnetosphere, Geophys. Res. Lett., 46, 11584–11596, https://doi.org/10.1029/2019GL083180, 2019. a
Dungey, J. W.: Interplanetary magnetic field and the auroral zones, Phys. Rev. Lett., 6, 47–48, https://doi.org/10.1103/PhysRevLett.6.47, 1961. a
Ergun, R. E., Holmes, J. C., Goodrich, K. A., Wilder, F. D., Stawarz, J. E., Eriksson, S., Newman, D. L., Schwartz, S. J., Goldman, M. V., Sturner, A. P., Malaspina, D. M., Usanova, M. E., Torbert, R. B., Argall, M., Lindqvist, P. A., Khotyaintsev, Y., Burch, J. L., Strangeway, R. J., Russell, C. T., Pollock, C. J., Giles, B. L., Dorelli, J. J. C., Avanov, L., Hesse, M., Chen, L. J., Lavraud, B., Le Contel, O., Retino, A., Phan, T. D., Eastwood, J. P., Oieroset, M., Drake, J., Shay, M. A., Cassak, P. A., Nakamura, R., Zhou, M., AshourAbdalla, M., and André, M.: Magnetospheric Multiscale observations of largeamplitude, parallel, electrostatic waves associated with magnetic reconnection at the magnetopause, Geophys. Res. Lett., 43, 5626–5634, https://doi.org/10.1002/2016GL068992, 2016. a
Escoubet, C. P., Fehringer, M., and Goldstein, M.: Introduction – The Cluster mission, Ann. Geophys., 19, 1197–1200, https://doi.org/10.5194/angeo1911972001, 2001. a
Fargette, N., Lavraud, B., Øieroset, M., Phan, T. D., ToledoRedondo, S., Kieokaew, R., Jacquey, C., Fuselier, S. A., Trattner, K. J., Petrinec, S., Hasegawa, H., Garnier, P., Génot, V., Lenouvel, Q., Fadanelli, S., Penou, E., Sauvaud, J. A., Avanov, D. L. A., Burch, J., Chand ler, M. O., Coffey, V. N., Dorelli, J., Eastwood, J. P., Farrugia, C. J., Gershman, D. J., Giles, B. L., Grigorenko, E., Moore, T. E., Paterson, W. R., Pollock, C., Saito, Y., Schiff, C., and Smith, S. E.: On the Ubiquity of Magnetic Reconnection Inside Flux Transfer EventLike Structures at the Earth's Magnetopause, Geophys. Res. Lett., 47, e86726, https://doi.org/10.1029/2019GL086726, 2020. a
Grandin, M., Battarbee, M., Osmane, A., Ganse, U., PfauKempf, Y., Turc, L., Brito, T., Koskela, T., Dubart, M., and Palmroth, M.: HybridVlasov modelling of nightside auroral proton precipitation during southward interplanetary magnetic field conditions, Ann. Geophys., 37, 791–806, https://doi.org/10.5194/angeo377912019, 2019. a
Grigorenko, E. E., Kronberg, E. A., Daly, P. W., Ganushkina, N. Y., Lavraud, B., Sauvaud, J.A., and Zelenyi, L. M.: Origin of low protontoelectron temperature ratio in the Earth's plasma sheet, J. Geophys. Res.Space, 121, 9985–10,004, https://doi.org/10.1002/2016JA022874, 2016. a
Hada, T., Nishida, A., Teresawa, T., and Hones Jr., E. W.: Bidirectional electron pitch angle anisotropy in the plasma sheet, J. Geophys. Res.Space, 86, 11211–11224, https://doi.org/10.1029/JA086iA13p11211, 1981. a, b
Harris, E. G.: On a plasma sheath separating regions of oppositely directed magnetic field, Il Nuovo Cimento, 23, 115–121, https://doi.org/10.1007/BF02733547, 1962. a
Hesse, M., Kuznetsova, M., Schindler, K., and Birn, J.: Threedimensional modeling of electron quasiviscous dissipation in guidefield magnetic reconnection, Phys. Plasmas, 12, 100704, https://doi.org/10.1063/1.2114350, 2005. a, b
Hesse, M., Liu, Y. H., Chen, L. J., Bessho, N., Kuznetsova, M., Birn, J., and Burch, J. L.: On the electron diffusion region in asymmetric reconnection with a guide magnetic field, Geophys. Res. Lett., 43, 2359–2364, https://doi.org/10.1002/2016GL068373, 2016. a, b
Hoilijoki, S., Ganse, U., PfauKempf, Y., Cassak, P. A., Walsh, B. M., Hietala, H., von Alfthan, S., and Palmroth, M.: Reconnection rates and X line motion at the magnetopause: Global 2D3V hybridVlasov simulation results, J. Geophys. Res.Space, 122, 2877–2888, https://doi.org/10.1002/2016JA023709, 2017. a
Hoilijoki, S., Ergun, R. E., Schwartz, S. J., Eriksson, S., Wilder, F. D., Webster, J. M., Ahmadi, N., Le Contel, O., Burch, J. L., Torbert, R. B., Strangeway, R. J., and Giles, B. L.: ElectronScale Magnetic Structure Observed Adjacent to an Electron Diffusion Region at the Dayside Magnetopause, J. Geophys. Res.Space, 124, 10153–10169, https://doi.org/10.1029/2019JA027192, 2019a. a
Hoilijoki, S., Ganse, U., Sibeck, D. G., Cassak, P. A., Turc, L., Battarbee, M., Fear, R. C., BlancoCano, X., Dimmock, A. P., Kilpua, E. K. J., Jarvinen, R., Juusola, L., PfauKempf, Y., and Palmroth, M.: Properties of Magnetic Reconnection and FTEs on the Dayside Magnetopause With and Without Positive IMF B_{x} Component During Southward IMF, J. Geophys. Res.Space, 124, 4037–4048, https://doi.org/10.1029/2019JA026821, 2019b. a
Hoshino, M., Hiraide, K., and Mukai, T.: Strong electron heating and nonMaxwellian behavior in magnetic reconnection, Earth Planets Space, 53, 627–634, https://doi.org/10.1186/BF03353282, 2001. a, b
Huang, S. Y., Jiang, K., Yuan, Z. G., Sahraoui, F., He, L. H., Zhou, M., Fu, H. S., Deng, X. H., He, J. S., Cao, D., Yu, X. D., Wang, D. D., Burch, J. L., Pollock, C. J., and Torbert, R. B.: Observations of the Electron Jet Generated by Secondary Reconnection in the Terrestrial Magnetotail, Astrophys. J., 862, 144, https://doi.org/10.3847/15384357/aacd4c, 2018. a
Huang, Z., Tóth, G., van der Holst, B., Chen, Y., and Gombosi, T.: A sixmoment multifluid plasma model, J. Comput. Phys., 387, 134–153, https://doi.org/10.1016/j.jcp.2019.02.023, 2019. a
Janhunen, P., Palmroth, M., Laitinen, T., Honkonen, I., Juusola, L., Facsko, G., and Pulkkinen, T.: The GUMICS4 global MHD magnetosphereionosphere coupling simulation, J. Atmos. Sol.Terr. Phy., 80, 48–59, https://doi.org/10.1016/j.jastp.2012.03.006, 2012. a
Juusola, L., Hoilijoki, S., PfauKempf, Y., Ganse, U., Jarvinen, R., Battarbee, M., Kilpua, E., Turc, L., and Palmroth, M.: Fast plasma sheet flows and X line motion in the Earth's magnetotail: results from a global hybridVlasov simulation, Ann. Geophys., 36, 1183–1199, https://doi.org/10.5194/angeo3611832018, 2018a. a
Juusola, L., PfauKempf, Y., Ganse, U., Battarbee, M., Brito, T., Grandin, M., Turc, L., and Palmroth, M.: A possible source mechanism for magnetotail current sheet flapping, Ann. Geophys., 36, 1027–1035, https://doi.org/10.5194/angeo3610272018, 2018b. a
Karimabadi, H., Roytershteyn, V., Vu, H. X., Omelchenko, Y. A., Scudder, J., Daughton, W., Dimmock, A., Nykyri, K., Wan, M., Sibeck, D., Tatineni, M., Majumdar, A., Loring, B., and Geveci, B.: The link between shocks, turbulence, and magnetic reconnection in collisionless plasmas, Phys. Plasmas, 21, 062308, https://doi.org/10.1063/1.4882875, 2014. a
Kempf, Y., Pokhotelov, D., Von Alfthan, S., Vaivads, A., Palmroth, M., and Koskinen, H. E. J.: Wave dispersion in the hybridVlasov model: Verification of Vlasiator, Phys. Plasmas, 20, 1–9, https://doi.org/10.1063/1.4835315, 2013. a
Khotyaintsev, Y. V., Cully, C. M., Vaivads, A., André, M., and Owen, C. J.: Plasma Jet Braking: Energy Dissipation and Nonadiabatic Electrons, Phys. Rev. Lett., 106, 165001, https://doi.org/10.1103/PhysRevLett.106.165001, 2011. a
Kilian, P., Muñoz, P. A., Schreiner, C., and Spanier, F.: Plasma waves as a benchmark problem, J. Plasma Phys., 83, 707830101, https://doi.org/10.1017/S0022377817000149, 2017. a
Lapenta, G., Markidis, S., Marocchino, A., and Kaniadakis, G.: Relaxation of Relativistic Plasmas Under the Effect of WaveParticle Interactions, Astrophys. J., 666, 949–954, https://doi.org/10.1086/520326, 2007. a
Lapenta, G., Markidis, S., Divin, A., Goldman, M., and Newman, D.: Scales of guide field reconnection at the hydrogen mass ratio, Phys. Plasmas, 17, 082106, https://doi.org/10.1063/1.3467503, 2010. a
Lapenta, G., Markidis, S., Goldman, M. V., and Newman, D. L.: Secondary reconnection sites in reconnectiongenerated flux ropes and reconnection fronts, Nat. Phys., 11, 690–695, https://doi.org/10.1038/nphys340, 2015. a
Li, W. Y., Graham, D. B., Khotyaintsev, Y. V., Vaivads, A., André, M., Min, K., Liu, K., Tang, B. B., Wang, C., Fujimoto, K., Norgren, C., ToledoRedondo, S., Lindqvist, P.A., Ergun, R. E., Torbert, R. B., Rager, A. C., Dorelli, J. C., Gershman, D. J., Giles, B. L., Lavraud, B., Plaschke, F., Magnes, W., Contel, O. L., Russell, C. T., and Burch, J. L.: Electron Bernstein waves driven by electron crescents near the electron diffusion region, Nat. Commun, 11, 1–10, https://doi.org/10.1038/s4146701913920w, 2020. a
Lin, Y. and Wang, X. Y.: Threedimensional global hybrid simulation of dayside dynamics associated with the quasiparallel bow shock, J. Geophys. Res., 110, A12216, https://doi.org/10.1029/2005JA011243, 2005. a
Lin, Z. and Chen, L.: A fluid–kinetic hybrid electron model for electromagnetic simulations, Phys. Plasmas, 8, 1447–1450, https://doi.org/10.1063/1.1356438, 2001. a
Liu, Y.H. H., Daughton, W., Karimabadi, H., Li, H., and Roytershteyn, V.: Bifurcated Structure of the Electron Diffusion Region in ThreeDimensional Magnetic Reconnection, Phys. Rev. Lett., 110, 265004, https://doi.org/10.1103/PhysRevLett.110.265004, 2013. a
Londrillo, P. and Del Zanna, L.: On the divergencefree condition in Godunovtype schemes for ideal magnetohydrodynamics: the upwind constrained transport method, J. Comput. Phys., 195, 17–48, https://doi.org/10.1016/j.jcp.2003.09.016, 2004. a, b
Lu, S., Lin, Y., Angelopoulos, V., Artemyev, A. V., Pritchett, P. L., Lu, Q., and Wang, X. Y.: Hall effect control of magnetotail dawndusk asymmetry: A threedimensional global hybrid simulation, J. Geophys. Res.Space, 121, 11882–11895, https://doi.org/10.1002/2016JA023325, 2016. a
Lu, S., Artemyev, A. V., Angelopoulos, V., Lin, Y., Zhang, X.J., Liu, J., Avanov, L. A., Giles, B. L., Russell, C. T., and Strangeway, R. J.: The Hall Electric Field in Earth's Magnetotail Thin Current Sheet, J. Geophys. Res.Space, 124, 1052–1062, https://doi.org/10.1029/2018JA026202, 2019. a
Mozer, F. S., Agapitov, O. V., Artemyev, A., Drake, J. F., Krasnoselskikh, V., Lejosne, S., and Vasko, I.: Time domain structures: What and where they are, what they do, and how they are made, Geophys. Res. Lett., 42, 3627–3638, https://doi.org/10.1002/2015GL063946, 2015. a
Nakamura, R., Baumjohann, W., Fujimoto, M., Asano, Y., Runov, A., Owen, C. J., Fazakerley, A. N., Klecker, B., Rème, H., Lucek, E. A., Andre, M., and Khotyaintsev, Y.: Cluster observations of an ionscale current sheet in the magnetotail under the presence of a guide field, J. Geophys. Res.Space, 113, A07S16, https://doi.org/10.1029/2007JA012760, 2008. a
Ni, B., Thorne, R. M., Zhang, X., Bortnik, J., Pu, Z., Xie, L., Hu, Z.j., Han, D., Shi, R., Zhou, C., and Gu, X.: Origins of the Earth's Diffuse Auroral Precipitation, Space Sci. Rev., 200, 205–259, https://doi.org/10.1007/s1121401602347, 2016. a
Nunn, D.: Vlasov Hybrid Simulation – An Efficient and Stable Algorithm for the Numerical Simulation of Collision‐Free Plasma, Transport Theor. Stat., 34, 151–171, https://doi.org/10.1080/00411450500255518, 2005. a
Omidi, N., Phan, T., and Sibeck, D. G.: Hybrid simulations of magnetic reconnection initiated in the magnetosheath, J. Geophys. Res.Space, 114, A02222, https://doi.org/10.1029/2008JA013647, 2009. a
Onsager, T. G., Thomsen, M. F., Elphic, R. C., and Gosling, J. T.: Model of electron and ion distributions in the plasma sheet boundary layer, J. Geophys. Res. Space, 96, 20999–21011, https://doi.org/10.1029/91JA01983, 1991. a, b, c
Onsager, T. G., Thomsen, M. F., Elphic, R. C., Gosling, J. T., Anderson, R. R., and Kettmann, G.: Electron generation of electrostatic waves in the plasma sheet boundary layer, J. Geophys. Res.Space, 98, 15509–15519, https://doi.org/10.1029/93JA00921, 1993. a
Palmroth, M.: Vlasiator, available at: http://www.physics.helsinki.fi/vlasiator/, last access: 25 January 2021. a
Palmroth, M. and the Vlasiator team: Vlasiator: hybridVlasov simulation code, Github repository, https://doi.org/10.5281/zenodo.3640593, version 4.0 and the eVlasiator branch, 2020. a
Palmroth, M., Hoilijoki, S., Juusola, L., Pulkkinen, T., Hietala, H., PfauKempf, Y., Ganse, U., von Alfthan, S., Vainio, R., and Hesse, M.: Tail reconnection in the global magnetospheric context: Vlasiator first results, Ann. Geophys., 35, 1269–1274, https://doi.org/10.5194/angeo3512692017, 2017. a
Palmroth, M., Ganse, U., PfauKempf, Y., Battarbee, M., Turc, L., Brito, T., Grandin, M., Hoilijoki, S., Sandroos, A., and von Alfthan, S.: Vlasov methods in space physics and astrophysics, Living Reviews in Computational Astrophysics, 4, 1, https://doi.org/10.1007/s4111501800032, 2018. a, b, c, d, e, f, g, h
Paterson, W. R. and Frank, L. A.: Survey of plasma parameters in Earth's distant magnetotail with the Geotail spacecraft, Geophys. Res. Lett., 21, 2971–2974, https://doi.org/10.1029/94GL02105, 1994. a
Pezzi, O., Cozzani, G., Califano, F., Valentini, F., Guarrasi, M., Camporeale, E., Brunetti, G., Retinò, A., and Veltri, P.: ViDA: a Vlasov–DArwin solver for plasma physics at electron scales, J. Plasma Phys., 85, 905850506, https://doi.org/10.1017/S0022377819000631, 2019. a
Phan, T. D., Eastwood, J. P., Shay, M. A., Drake, J. F., Sonnerup, B. U. Ö., Fujimoto, M., Cassak, P. A., Øieroset, M., Burch, J. L., Torbert, R. B., Rager, A. C., Dorelli, J. C., Gershman, D. J., Pollock, C., Pyakurel, P. S., Haggerty, C. C., Khotyaintsev, Y., Lavraud, B., Saito, Y., Oka, M., Ergun, R. E., Retino, A., Le Contel, O., Argall, M. R., Giles, B. L., Moore, T. E., Wilder, F. D., Strangeway, R. J., Russell, C. T., Lindqvist, P. A., and Magnes, W.: Electron magnetic reconnection without ion coupling in Earth's turbulent magnetosheath, Nature, 557, 202–206, https://doi.org/10.1038/s4158601800915, 2018. a
Pritchett, P.: Particleincell simulations of magnetosphere electrodynamics, IEEE T. Plasma Sci., 28, 1976–1990, https://doi.org/10.1109/27.902226, 2000. a
Ricci, P., Lapenta, G., and Brackbill, J. U.: GEM reconnection challenge: Implicit kinetic simulations with the physical mass ratio, Geophys. Res. Lett., 29, 3–1–3–4, https://doi.org/10.1029/2002GL015314, 2002. a
Runov, A., Angelopoulos, V., Gabrielse, C., Liu, J., Turner, D. L., and Zhou, X.Z.: Average thermodynamic and spectral properties of plasma in and around dipolarizing flux bundles, J. Geophys. Res.Space, 120, 4369–4383, https://doi.org/10.1002/2015JA021166, 2015. a
Sandroos, A.: VLSV: file format and tools, Github repository, available at: https://github.com/fmihpc/vlsv/ (last access: 30 November 2020), 2019. a
Schmitz, H. and Grauer, R.: Kinetic Vlasov simulations of collisionless magnetic reconnection, Phys. Plasmas, 13, 092309, https://doi.org/10.1063/1.2347101, 2006. a
Sibeck, D. G., Omidi, N., Dandouras, I., and Lucek, E.: On the edge of the foreshock: modeldata comparisons, Ann. Geophys., 26, 1539–1544, https://doi.org/10.5194/angeo2615392008, 2008. a
Swisdak, M.: Quantifying gyrotropy in magnetic reconnection, Geophys. Res. Lett., 43, 43–49, https://doi.org/10.1002/2015GL066980, 2016. a, b
Tronci, C. and Camporeale, E.: Neutral Vlasov kinetic theory of magnetized plasmas, Phys. Plasmas, 22, 020704, https://doi.org/10.1063/1.4907665, 2015. a
Tóth, G., Jia, X., Markidis, S., Peng, I. B., Chen, Y., Daldorff, L. K. S., Tenishev, V. M., Borovikov, D., Haiducek, J. D., Gombosi, T. I., Glocer, A., and Dorelli, J. C.: Extended magnetohydrodynamics with embedded particleincell simulation of Ganymede's magnetosphere, J. Geophys. Res.Space, 121, 1273–1293, https://doi.org/10.1002/2015JA021997, 2016. a
Tóth, G., Chen, Y., Gombosi, T. I., Cassak, P., Markidis, S., and Peng, I. B.: Scaling the Ion Inertial Length and Its Implications for Modeling Reconnection in Global Simulations, J. Geophys. Res.Space, 122, 10336–10355, https://doi.org/10.1002/2017JA024189, 2017. a
Umeda, T., Togano, K., and Ogino, T.: Twodimensional fullelectromagnetic Vlasov code with conservative scheme and its application to magnetic reconnection, Comput. Phys. Commun., 180, 365–374, https://doi.org/10.1016/j.cpc.2008.11.001, 2009. a
von Alfthan, S., Pokhotelov, D., Kempf, Y., Hoilijoki, S., Honkonen, I., Sandroos, A., and Palmroth, M.: Vlasiator: First global hybridVlasov simulations of Earth's foreshock and magnetosheath, J. Atmos. Sol.Terr. Phy., 120, 24–35, https://doi.org/10.1016/j.jastp.2014.08.012, 2014. a, b
Wang, C.P., Gkioulidou, M., Lyons, L. R., and Angelopoulos, V.: Spatial distributions of the ion to electron temperature ratio in the magnetosheath and plasma sheet, J. Geophys. Res.Space, 117, A08215, https://doi.org/10.1029/2012JA017658, 2012. a, b
Wang, J., Huang, C., Ge, Y. S., Du, A., and Feng, X.: Influence of the IMF B_{x} on the geometry of the bow shock and magnetopause, Planet. Space Sci., 182, 104844, https://doi.org/10.1016/j.pss.2020.104844, 2020. a
Wang, L., Hakim, A. H., Bhattacharjee, A., and Germaschewski, K.: Comparison of multifluid moment models with particleincell simulations of collisionless magnetic reconnection, Phys. Plasmas, 22, 012108, https://doi.org/10.1063/1.4906063, 2015. a
Wilson, F., Neukirch, T., Hesse, M., Harrison, M. G., and Stark, C. R.: Particleincell simulations of collisionless magnetic reconnection with a nonuniform guide field, Phys. Plasmas, 23, 032302, https://doi.org/10.1063/1.4942939, 2016. a
Yamamoto, T. and Tamao, T.: Adiabatic plasma convection in the tail plasma sheet, Planet. Space Sci., 26, 1185–1191, https://doi.org/10.1016/00320633(78)900582, 1978. a
Zerroukat, M. and Allen, T.: A threedimensional monotone and conservative semiLagrangian scheme (SLICE3D) for transport problems, Q. J. Roy. Meteorol. Soc., 138, 1640–1651, https://doi.org/10.1002/qj.1902, 2012. a
Zhang, X., Angelopoulos, V., Artemyev, A. V., and Liu, J.: Whistler and Electron Firehose Instability Control of Electron Distributions in and Around Dipolarizing Flux Bundles, Geophys. Res. Lett., 45, 9380–9389, https://doi.org/10.1029/2018GL079613, 2018. a
Zhou, H., Tóth, G., Jia, X., Chen, Y., and Markidis, S.: Embedded Kinetic Simulation of Ganymede's Magnetosphere: Improvements and Inferences, J. Geophys. Res.Space, 124, 5441–5460, https://doi.org/10.1029/2019JA026643, 2019. a