Ion distributions upstream and downstream of the Earth ’ s bow shock : first results from Vlasiator

A novel hybrid-Vlasov code, Vlasiator, is developed for global simulations of magnetospheric plasma kinetics. The code is applied to model the collisionless bow shock on scales of the Earth’s magnetosphere in two spatial dimensions and three dimensions in velocity space retrieving ion distribution functions over the entire foreshock and magnetosheath regions with unprecedented detail. The hybrid-Vlasov approach produces noise-free uniformly discretized ion distribution functions comparable to those measured in situ by spacecraft. Vlasiator can reproduce features of the ion foreshock and magnetosheath well known from spacecraft observations, such as compressional magnetosonic waves generated by backstreaming ion populations in the foreshock and mirror modes in the magnetosheath. An overview of ion distributions from various regions of the bow shock is presented, demonstrating the great opportunities for comparison with multi-spacecraft observations.


Introduction
Global modeling of the Earth's magnetosphere became feasible with rapidly increasing computational power.Naturally, the computationally least demanding magnetohydrodynamic (MHD) approach evolved first with a number of mature global MHD codes now available for space weather research and operational forecast.Unfortunately, many large-scale magnetospheric processes cannot be successfully reproduced in MHD models.One such example is a complex waveparticle interaction pattern observed in the ion foreshock region in front of the Earth's bow shock, which emerges due to resonant interactions between the solar wind and backstreaming particles that are reflected and energized by the collisionless bow shock (Scholer et al., 1993;Eastwood et al., 2005a).Collisionless energization and particle reflection in the bow shock include non-MHD kinetic processes, and the turbulent magnetosheath formed behind the shock front also has fluctuations coupled to the kinetic scales.Thus a beyond-MHD kinetic approach is required to model the solar wind-bow shock interaction on a global scale.
Two distinct approaches are currently used to implement kinetic physics into global models.The first approach is based on coupling kinetic and MHD models for respective magnetospheric regions (e.g., Tóth et al., 2012), while the second approach is focused on solving the kinetic equations throughout the entire system.The former approach was successfully applied in modeling the inner magnetosphere and ring current region (Glocer et al., 2013) where collisional processes and multi-ion plasma composition become significant.The latter approach is typically based on hybrid particle-in-cell (PIC) simulations, where electrons are modeled as fluid while ions are modeled as macroparticles for which plasma kinetic equations are solved (Winske et al., 2003).Such simulations have been successful in modeling the Earth's bow shock and ion foreshock in two spatial dimensions (Lin, 2003;Omidi et al., 2005;Blanco-Cano et al., 2006;Omelchenko and Karimabadi, 2012) and occasionally in three spatial dimensions (Lin and Wang, 2005) as well as in modeling solar wind interaction with other solar system bodies (Brecht and Ferrante, 1991;Kallio and Janhunen, 2002).However, the hybrid-PIC approach may yield solutions with an undesirable level of numerical noise, especially when looking at ion velocity distribution functions.

D. Pokhotelov et al.: Ion distributions near bow shock
Numerical noise can be reduced by increasing the number of ions launched in the simulation, but this leads to a rapid increase in computational demands.Also, most hybrid-PIC simulations were applied to small-scale planetary magnetospheres or used a downscaled Earth's magnetic dipole to reduce computational demands.
An alternative approach to the hybrid-PIC modeling is to use Vlasov's equation and model directly the evolution of the six-dimensional ion distribution function (three in ordinary space, and three in velocity space) while treating electrons as a fluid, yielding a hybrid-Vlasov approach.Hybrid-Vlasov simulations require massive amounts of memory and computations to propagate the ion distribution function and are, thus, even more computationally challenging than the hybrid-PIC approach.The main benefit of the hybrid-Vlasov approach is the ability to produce noise-free uniformly discretized ion distribution functions comparable to those measured in situ by spacecraft (e.g., Kis et al., 2007).Due to the computational complexity, the global hybrid-Vlasov approach has not received a lot of attention in the past, while a number of local simulations have been implemented (Valentini et al., 2007;Valentini et al., 2010;Eliasson and Shukla, 2007).The newly developed global hybrid-Vlasov code, Vlasiator (Sandroos et al., 2013;Palmroth et al., 2013), has been utilized in this study to model the kinetic behavior of the collisionless bow shock, mainly focusing on ion distribution functions upstream and downstream of the Earth's bow shock.

Numerical model
Vlasiator (Sandroos et al., 2013;Palmroth et al., 2013) is a novel space plasma simulation code designed for global magnetospheric simulations.It is based on a selfconsistent hybrid-Vlasov model describing ions with a distribution function while treating electrons as a massless charge-neutralizing fluid.The fundamental description of charged particle motion in an electromagnetic field is given by Vlasov's equation where r and v are the spatial and velocity coordinates, f = f (r, v, t) is the six-dimensional phase-space density of a particle species with mass m and charge q, E is the electric field and B is the magnetic field.The bulk parameters of the plasma, such as the ion charge density ρ q , are obtained as velocity moments of the ion velocity distribution function.
Vlasov's equation is coupled with Maxwell's equations in which the displacement current is neglected.The system is closed by Ohm's law that models the effect of the massless electron fluid.It describes the dependence of the electric field on the magnetic field, and is needed when updating the magnetic field using Faraday's law.Vlasiator currently uses an ideal Ohm's law, modeling electrons as a massless charge-neutralizing fluid: where V i is the ion bulk velocity.The electric field used when propagating fields does not include the Hall term j × B/ρ q , where j is the total current density.This is a reasonably good approximation for the global simulations presented in this study, since we are not aiming to resolve the inertial ion scales and thus the Hall term can be neglected.In the Lorentz force, however, E is given by where the second term on the right-hand side is the Hall term.
Here the Hall term needs to be included, since otherwise momentum conservation is violated and no bulk force is exerted on the ions (Karimabadi et al., 2004).
Vlasiator uses a second-order accurate finite volume method to solve Vlasov's equation and supports full sixdimensional cases (three spatial and three velocity coordinates).The method uses Strang splitting (Strang, 1968) to separate propagation in ordinary and velocity spaces, and is based on a three-dimensional wave-propagation algorithm (Langseth and LeVeque, 2000;LeVeque, 1997).The field solver is a second-order accurate upwind constrained transport method, which is divergence-conserving by construction (Londrillo and del Zanna, 2004).
Global hybrid-Vlasov simulations require massively parallel computations on supercomputers to propagate the full sixdimensional distribution functions.Even five-dimensional simulations (two spatial coordinates, three velocity coordinates) described in this study consist of the order of 10 12 potential phase-space cells.To enable these simulations the Vlasiator code has been parallelized using a two-level strategy with the Message Passing Interface (MPI) and OpenMP.The current code scales well to tens of thousands of cores, and the simulations presented here have been computed on 16 384 cores on a Cray XE6 supercomputer.To reduce the computational load and memory usage the ion distribution function is described by a sparse representation where effectively empty phase-space cells are neither propagated nor stored.This reduces the number of phase-space cells in these simulations to 1 % of its original value.

Global magnetospheric simulation
The global magnetospheric simulation is set up to cover the region of near-Earth space from the inner boundary (defined at 6 R E radial distance) to the solar wind, thus covering the dayside part of the Earth's magnetosphere, the magnetosheath, the bow shock and the foreshock region.The Geocentric Solar Ecliptic System (GSE) is used with its x axis pointing from the Earth towards the Sun, its y axis in the ecliptic plane pointing towards dusk and its z axis parallel to the ecliptic pole.The simulation box is shifted towards the Sun, extending from −20 to 40 R E along the x axis and from −67 to 52 R E along the y axis.The grid dimensions in two-dimensional ordinary space are 450 and 900 along the x and y axis, respectively, thus making the spatial resolution 0.13 R E .Along the z axis the system is periodic.The velocity space has a resolution of 20 km s −1 , and its extent in all three velocity dimensions is −2000 km s −1 to 2000 km s −1 .The Earth's magnetic dipole is directed along the z axis, and has a strength of 8.0 × 10 22 Am 2 .
At the solar wind boundary on the edge facing the Sun a static Maxwellian ion velocity distribution is set in all cells, corresponding to typical solar wind conditions (V SW = 500 km s −1 , n i = 10 6 m −3 , T i = 10 5 K).The simulation is set up with an IMF of 5 nT pointing at a 45 • angle with respect to the x axis (IMF B y = −B x = 3.5 nT).To set the IMF we add a constant magnetic field to the background field in all cells including the boundary cells.At the inner boundary the distribution function is set to a stationary state (n i = 10 6 m −3 , T i = 10 5 K) and the perturbed magnetic field is set to zero.This inner boundary is far from being realistic with respect to the actual inner magnetosphere, and can be viewed as a first implementation that enables studies of phenomena, such as the ion foreshock, which are located far from the inner boundary.At the outflow boundaries each cell copies from its nearest normal neighbor cell the distribution function, and its magnetic field value.This allows plasma to flow out of the system.
In this simulation the ion thermal scales are well resolved in the velocity space with the ion thermal velocities of ∼ 65 km s −1 and 200 km s −1 in the ion foreshock and in the magnetosheath, respectively, compared with the velocity space resolution of 20 km s −1 .However, no attempt was made to resolve the thermal ion gyroscales of ∼ 200 km.The impact of spatial resolution on hybrid-Vlasov simulations is assessed in more detail in a separate study using local simulation with Vlasiator (Kempf et al., 2013).
After the initialization stage of ∼ 200 s required for the solar wind to reach the Earth's magnetic dipole, the collisionless bow shock forms in front of the Earth.At the 45 • IMF cone angle the bow shock consists of two distinct regions: the quasi-perpendicular shock region in the dusk sector and the quasi-parallel shock region in the dawn sector. Figure 1 shows perturbations of the ion density throughout the computational domain at t = 1027 s.The ion foreshock boundary is clearly seen as a line separating the regions of quasi-parallel and quasi-perpendicular bow shock.The quasi-perpendicular shock is unable to reflect particles into the solar wind and thus the ion population in the dusk sector consists only of the Maxwellian solar wind flow.In contrast, the quasi-parallel bow shock is characterized by reflected ion populations accelerated by the shock and streaming back against the incoming solar wind.
In the presented simulation, backstreaming ion distributions appearing across the ion foreshock region can be classified as follows.Gyrating ring-shaped ion distributions with relatively narrow energy spectra appear in close vicinity of the ion foreshock boundary (panel b in Fig. 1).The shape of the ion distribution function is shown in Fig. 1b as a 2-D slice in the xy velocity plane and also as a 3-D isosurface.Downstream from the ion foreshock boundary (panel c in Fig. 1) the backstreaming ion populations appear to be cap-shaped as shown in Fig. 1c, with the energy spectra getting broader and more diffused further downstream (see Fig. 1d).
The bulk velocity of backstreaming ion populations steadily decreases towards the inner region of the ion foreshock.Backstreaming ion velocities V b normalized by solar wind velocity range from V b /V SW = 2.0-2.2near the ion foreshock boundary to V b /V SW = 0.8-0.9deep inside the ion foreshock.The backstreaming ion populations seen across the ion foreshock are accompanied by quasi-monochromatic compressional oscillations with periods ∼ 30 s and wavelengths varying from ∼ 2 R E near the ion foreshock boundary to less than 1 R E deep inside the ion foreshock.The wave magnitude is typically 5-10 % of the background value both for the magnetic field and density, confirming the compressional nature of the waves, which is consistent with theoretical predictions for the fast magnetosonic perturbations in the MHD limit (e.g., Le and Russell, 1994).
Properties of ion distributions reproduced in the simulation are generally consistent with known single-and multispacecraft observations in the foreshock as summarized by Bonifazi and Moreno (1981); Fuselier (1995); Eastwood et al. (2005a).Narrow gyrating, often gyrophase-bunched, distributions are typically observed near the ion foreshock boundary.Cap-shaped (crescent-shaped in 2-D) ion distributions, also known as intermediate distributions (Fuselier, 1995), are observed deeper in the ion foreshock, with more diffused ion distributions appearing further downstream.The intermediate distributions are associated with quasimonochromatic foreshock waves known as 30 s sinusoidal waves (Le and Russell, 1994) believed to be generated by the backstreaming shock-energized ion population via ion/ion resonance interaction (Gary, 1991).Multi-spacecraft analysis using Cluster satellites shows these waves to be compressional magnetosonic modes with periods close to 30 s and wavelengths of ∼ 1 R E propagating obliquely to the background magnetic field (Eastwood et al., 2005b) in agreement with the picture reproduced in our simulation.Variations in the energy of backstreaming populations across the ion foreshock seen in Fig. 1 can be understood in terms of the efficiency of collisionless shock acceleration (which is expected to be highest near the ion foreshock boundary) and is also consistent with the observational statistics.For instance, Bonifazi and Moreno (1981) reported that V b /V SW are around 2, 1.7, and 1.2 for field-aligned, intermediate, and diffuse distributions, respectively, which is also consistent with more recent multi-spacecraft statistics (Eastwood et al., 2005b).
The magnetosheath is formed behind the bow shock and is composed of highly anisotropic shock-energized plasmas.Ion velocity distributions change from horseshoe-shaped in the close vicinity of the bow shock to gyrotropic bi-Maxwellian distributions deep inside the magnetosheath (panel a in Fig. 1).Such bi-Maxwellian ion distributions with high-temperature anisotropy transverse to the background magnetic field are known to be a subject to mirror mode instability as well as to ion-cyclotron instability (Southwood and Kivelson, 1993).Profiles of plasma parameters across the magnetosheath (not shown) clearly demonstrate the anticorrelation between plasma density and magnetic field enhancements that is characteristic of perturbations induced by mirror mode instability.Such mirror mode structures with spatial scales of a few ion inertial lengths are routinely observed across the Earth's magnetosheath (Soucek et al., 2008).The loss cone seen in the magnetosheath ion distribution is also a characteristic feature of the mirror modes and is due to ion trapping between magnetic mirrors (Soucek and Escoubet, 2011).

Conclusions and summary
We present the first results of self-consistent global hybrid-Vlasov simulations of the magnetospheric bow shock and ion foreshock.We demonstrate that the new code, Vlasiator, is able to reproduce the key features of solar windmagnetosphere interactions such as ion reflection and energization at the collisionless bow shock and to simulate associated wave-particle interaction processes taking place in the ion foreshock and in the magnetosheath.In the quasiparallel region of the ion foreshock, reflected and energized solar wind ions form backstreaming ion distributions of gyrating/beam type near the ion foreshock boundary and capshaped (intermediate) type distributions deep in the ion foreshock.Quasi-monochromatic compressional magnetosonic waves with periods of ∼ 30 s propagating obliquely to the background magnetic field appear across the ion foreshock region clearly associated with the backstreaming ion populations.Characteristics of the backstreaming ion populations and associated electromagnetic waves are in agreement with the properties of ion velocity distributions and compressional magnetosonic waves typically observed in the Earth's ion foreshock region.In the magnetosheath region the horseshoe-shaped and gyrotropic bi-Maxwellian ion distributions appear near the bow shock and deep in the magnetosheath, respectively.The bi-Maxwellian ion distributions are associated with large-scale standing structures appearing across the magnetosheath.Anti-correlation between density and magnetic fields across the magnetosheath structures as well as the shape of ion distributions suggests the mirror mode instability as a generation mechanism for these structures, which is also consistent with multi-spacecraft observations of mirror modes in the Earth's magnetosheath.
The global simulation presented here demonstrates qualitative agreement with earlier five-dimensional simulations of the collisionless bow shock using hybrid-PIC codes (Omidi et al., 2005;Blanco-Cano et al., 2006).In the quasiperpendicular region the bow shock resembles a fast magnetosonic shock similar to that reproduced in hybrid-PIC simulations.In the quasi-parallel region the bow shock demonstrates rather complex spatio-temporal structure with shock rippling and reformation, but still retains the features of a fast magnetosonic shock in contrast to hybrid-PIC simulations by Omidi et al. (2005), where the quasi-parallel region appears as a magnetosonic wave followed by a rotational discontinuity.The most notable difference with respect to hybrid-PIC simulations is that Vlasiator ion velocity distributions appear as uniformly discretized functions similar to those seen in the experimental data, in contrast to statistically noisy distributions derived from hybrid-PIC simulations.Hybrid-PIC simulations of the ion foreshock reproduced similar quasimonochromatic 30 s wave structures.However, in the hybrid-PIC simulations (Blanco-Cano et al., 2006) the ion foreshock waves are mostly non-compressive Alfvénic oscillations, which appears to contradict the observed statistical characteristics of the 30 s ion foreshock waves known to be mostly compressive magnetosonic and oblique (Eastwood et al., 2005b).Other types of wave structures known to exist in the ion foreshock, such as steepened nonlinear shocklets (Le and Russell, 1994), have been reproduced earlier in global hybrid-PIC simulations (Omidi et al., 2005;Blanco-Cano et al., 2006) but do not appear in our simulations, possibly due to the simulation box not being large enough for the highamplitude nonlinear structures to evolve before the foreshock structures get advected by the solar wind flow into the bow shock.

D.
Fig. 1.Ion distributions (in units s 3 m −6 ) in the magnetosheath (A), near the ion foreshock boundary (B), and in the ion foreshock (C and D).