# Vlasov methods in space physics and astrophysics

Minna Palmroth, Urs Ganse, Yann Pfau-Kempf, Markus Battarbee, Markku Alho, Joonas Nättilä, Ivan Zaitsev, Giulia Cozzani, Konstantinos Papadakis, Leo Kotipalo, et al.

► To cite this version:

Minna Palmroth, Urs Ganse, Yann Pfau-Kempf, Markus Battarbee, Markku Alho, et al. Living Reviews in Computational Astrophysics, 2025, 11, 10.1007/s41115-025-00024-0. insu-05537449

HAL Id: insu-05537449

https://insu.hal.science/insu-05537449v1

Submitted on 5 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Vlasov methods in space physics and astrophysics

Minna Palmroth · Urs Ganse · Yann Pfau-Kempf · Markus Battarbee · Markku Alho · Joonas Nättilä et al. [full author details at the end of the article]

Received: 9 May 2025 / Accepted: 8 September 2025 / Published online: 14 November 2025
© The Author(s) 2025

## Abstract

This paper reviews Vlasov-based numerical methods used to model plasma in space physics and astrophysics. Plasma consists of collectively behaving charged particles that form the major part of baryonic matter in the Universe. Many physical concepts ranging from our own planetary environment to the Solar system and beyond can be understood in terms of kinetic plasma physics, represented by the Vlasov equation. We introduce the physical basis for the Vlasov systems, and then outline the associated numerical methods that are typically used. A particular application of the Vlasov system is Vlasiator, the world's first global hybrid-Vlasov simulation for the Earth's magnetic domain, the magnetosphere. We introduce the design strategies for Vlasiator and outline its numerical concepts ranging from solvers to coupling schemes. We review Vlasiator's parallelisation methods and introduce the used high-performance computing (HPC) techniques. A short review of verification, validation, and physical results is included. The purpose of the paper is to present the Vlasov equation and its use in numerical modelling, introduce an example implementation, and to illustrate that even with massive computational challenges, an accurate description of physics is highly rewarding as it allows to significantly advance our understanding of complex plasma systems. As an update to the previous paper in this series, we especially highlight latest trends in HPC including graphics processing units (GPUs) and emerging applications in astrophysical plasmas.

This article is a revised version of https://doi.org/10.1007/s41115-018-0003-2.

**Change summary:** Major revision, updated and expanded.

**Change details:** We have noticed the first serious attempts to use Vlasov-approach in astrophysics, and in anticipation of followers, we have here emphasized the astrophysics setup by inviting a new author from this field, converting the equations to the form that can be used in gravitational plasmas, and added examples and updates accordingly. Another new thing relative to the 2018 version is the heterogeneous supercomputers, so we have emphasized GPU parallelization. We have also briefly mentioned quantum and AI algorithms that are coming, and subgrid-methods, which are also coming. Otherwise we have kept the structure similar (not the same), and updated where things have changed. The number of references has increased from 240 to 432.
<table><tbody><tr><td>5.3</td><td>Solvers and time-integration</td><td>43</td></tr><tr><td>5.3.1</td><td>Vlasov solver</td><td>43</td></tr><tr><td>5.3.2</td><td>Field solver</td><td>44</td></tr><tr><td>5.3.3</td><td>Time stepping</td><td>45</td></tr><tr><td>5.4</td><td>Boundary conditions</td><td>46</td></tr><tr><td>5.4.1</td><td>Domain outer boundaries</td><td>46</td></tr><tr><td>5.4.2</td><td>Inner boundary</td><td>46</td></tr><tr><td>5.5</td><td>Initial conditions</td><td>47</td></tr><tr><td>5.6</td><td>Parallelisation strategies</td><td>49</td></tr><tr><td>5.6.1</td><td>Parallelisation on accelerated hardware</td><td>50</td></tr><tr><td>5.6.2</td><td>Parallel input / output</td><td>51</td></tr><tr><td>5.7</td><td>Verification and test problems</td><td>51</td></tr><tr><td>5.8</td><td>eVlasiator</td><td>52</td></tr><tr><td>5.9</td><td>Physics results</td><td>53</td></tr><tr><td>5.9.1</td><td>Terrestrial foreshock</td><td>53</td></tr><tr><td>5.9.2</td><td>Terrestrial magnetosheath</td><td>54</td></tr><tr><td>5.9.3</td><td>ULF wave transmission</td><td>54</td></tr><tr><td>5.9.4</td><td>Magnetic reconnection at the terrestrial magnetopause and magnetotail</td><td>55</td></tr><tr><td>5.9.5</td><td>Results based on global 3D-3V simulations</td><td>55</td></tr><tr><td>5.10</td><td>Future avenues</td><td>58</td></tr><tr><td>6</td><td>Conclusions and outlook</td><td>59</td></tr><tr><td></td><td>References</td><td>62</td></tr></tbody></table>

# 1 Introduction

While physical understanding of the natural world is inherently based on empirical evidence, numerical simulation tools have become an integral part of the majority of fields within physics. When tested against observations, numerical models can strengthen or invalidate existing theories and quantify the degree to which the theories have to be improved. Simulation results can also complement observations by giving them a larger context. In space physics, a satellite measures only one point at one time in the vast volume of space, and thus discerning spatial phenomena from temporal changes is difficult. This is a shortcoming that has also led to the use of satellite constellations, such as the European Space Agency Cluster mission (Escoubet et al. 2001), and the National Aeronautics and Space Administration (NASA) THEMIS and MMS missions (Angelopoulos 2008; Burch et al. 2016). Simulations are considerably more cost-effective than spacecraft, and they can be adopted to address physical systems that cannot be reached by *in situ* experiments, such as distant galaxies. Finally, and most importantly, predictions of physical environments under varying conditions are always based on modelling and numerical simulations. Accurate predictions include two types of modelling efforts: those targeting rapid solutions for forecast, and those which test and improve the current understanding together with top-tier experiments. This paper concerns the latter approach.

Predicting the physical conditions in the vicinity of the Earth in particular has become increasingly important because of the expensive satellite assets that form a critical infrastructure for the modern way of life (e.g. Eastwood et al. 2017; Koskinen et al. 2017). The space environmental conditions that threaten space- or ground-based
technology or human life are commonly termed as *space weather*. Energetic particles from the Sun alter the communication conditions especially at high latitudes, affecting radio broadcasts, aircraft communication with air traffic control, and radar signals. Sudden changes in the magnetic field induce currents in terrestrial long conductors, such as gas pipelines, railways, and power grids that can sometimes be disrupted (e.g. Wik et al. 2008). Increasing numbers of satellites are being launched, vulnerable to sudden events in the geospace. In 2022, the Starlink company lost 38 satellites due to relatively moderate space weather activity (Zhang et al. 2022). Overall, some estimations show that in the worst case, an extreme space weather event could induce huge economic costs (Oughton et al. 2017). Predicting the geospace conditions is ultimately based on modelling.

In physics of collisionless plasmas, the dominant physical interactions are caused by long-range interactions, most importantly the electromagnetic forces over a collection of charged particles. There are three main approaches to model plasmas: 1) the fluid approach (e.g., magnetohydrodynamics, MHD), 2) the fully kinetic approach, and 3) hybrid approaches combining the first two. The fluid approach leads to computationally lighter simulations, while fully kinetic simulations require significantly larger computing power. As a result, due to the size of the simulation domain, the first simulations of the near-Earth plasma system were based on MHD (e.g. Janhunen et al. 2012). However, single-fluid MHD models are basically scaleless in that they assume that plasmas have a single temperature approximated by a Maxwellian distribution. Therefore, they provide a limited context to space observations that are almost always multi-temperature and hence non-Maxwellian, as there are no collisions to thermalise the particle velocity distributions. The second approach uses a kinetic formulation as represented by the Vlasov theory (Vlasov 1961). In this approach, plasma species are treated as velocity distribution functions in a six-dimensional phase space consisting of a three-dimensional ordinary space (3D) and a three-dimensional velocity space (3V). The majority of kinetic simulations model the Vlasov theory by a particle-in-cell (PIC) method (Nishikawa et al. 2021), where a large number of particles are propagated within the simulation, and the distribution function is constructed from particle statistics in space and time. The fully kinetic PIC approach means that both electrons and ions are treated as particles within the simulation. Such simulations in 3D are computationally extremely costly, and while the first such endeavours were using local geometries such as current sheets (e.g. Daughton et al. 2011), the increased computational resources have recently enabled fully-kinetic PIC simulations of small magnetised planets (e.g., Ballerini et al. 2024; Lavorenti et al. 2022).

A hybrid approach in the kinetic simulation regime usually means that electrons are treated with a fluid description, but protons and possibly heavier ions are treated kinetically. Again, the vast majority of simulations use a hybrid-PIC approach (Nishikawa et al. 2021), which have recently been extended to 3D (e.g. Lin et al. 2017; Fatemi et al. 2024). This paper does not discuss the details of the PIC approach, but instead concentrates on the hybrid-Vlasov method, where the ion velocity distribution is discretised and modelled with a 3D-3V grid. The difference to hybrid-PIC is that in hybrid-Vlasov the distribution functions are evolved in time as an entity, and not constructed from particle statistics. The main advantage is therefore
that the distribution function is noiseless. This can be important for the problem at hand, because the distribution function is in many respects the core of plasma physics as the majority of the plasma parameters and processes can be derived from it. Hybrid-Vlasov methods were previously used mostly in local geometries, because the 3D-3V requirement implies a large computational cost. A global approach, which in space physics means simulation box sizes exceeding thousands of ion inertial lengths or gyroradii per dimension, was not possible until recently. The world's (so far) only global magnetospheric hybrid-Vlasov simulation, the massively parallel Vlasiator (Palmroth et al. 2018a), is therefore the prime application in this article.

This paper is organised as follows: Section 2 introduces the typical plasma systems in astrophysics and space physics, and the relevant processes that motivate the development of Vlasov-based methods. Section 3 introduces the Vlasov theory, while Sect. 4 introduces the numerical representations of the Vlasov theory. Section 5 describes Vlasiator in detail and justifies the decisions made in the design of the code to aid those who would like to design their own (hybrid-)Vlasov system. There are no standard verification cases for a (hybrid-)Vlasov system, but we describe the test cases used for Vlasiator. The physical findings are then illustrated briefly, showing that Vlasiator has made a paradigm change in space physics, emphasising the role of scale coupling in large-scale plasma systems. After the publication of the previous edition of this review (Palmroth et al. 2018a), Vlasov-based approaches have been increasingly applied to astrophysical systems, and hence we pay more attention to astrophysical plasmas in this version of the paper. Further, the computational architectures have changed to include also graphics processing units (GPUs), affecting the parallelisation strategies. Finally, new ideas have been brought up involving quantum algorithms and machine learning, which are also briefly covered here.

## 2 Kinetic physics in space and astrophysical plasmas

Interactions between charged particles and electromagnetic fields follow the same basic rules throughout the universe, but the applicability of simplified theories and the relevant spatial, temporal, and virial scales vary greatly between different scopes of research. In this section, we present an overview of regions of interest and the phenomena found within them. We briefly introduce the physics encountered within astrophysical setups (Sect. 2.1), Solar system (Sect. 2.2), focussing to near-Earth space and planetary environments (Sect. 2.3). The purpose is to provide examples of systems which can be studied with Vlasov simulations.

### 2.1 Astrophysical media and objects

Prime examples of themes requiring modelling are e.g., the dynamics of hot, cold, and dark matter in an expanding universe with unknown boundaries. The birth of the universe connects the rapid expansion and cooling of baryonic matter with quantum fluctuation anisotropies that eventually lead to the formation of galactic superclusters. Astrophysical simulations of the universe should naturally account for expansion of
space-time and associated effects of general relativity, and modelling of high-energy phenomena should correctly account for special relativity due to velocities approaching the speed of light. A forerunner in modelling the universe is EAGLE (Schaye et al. 2015), which utilises smoothed particle hydrodynamics, with subgrid modelling providing feedback of star formation, radiative cooling, stellar mass loss and feedback from stars and accreting black holes. These simulations operate on very much larger scales compared to the Vlasov equation for ions and electrons, yet they depend strongly on knowledge of processes at smaller length and time scales. Due to the majority of the universe consisting of the mostly empty interstellar and intergalactic media, the energy content of turbulent space plasmas must be understood. This has been investigated through the Vlasov equation (see, e.g., Weinstock 1969). Conversely, turbulent behaviour at large scales can act as a model for extending power laws to smaller scales (Maier et al. 2009). An alternative if less common approach for modelling galactic dynamics is to describe the distribution of stars as a Vlasov-Poisson system (see Sect. 3.3), with gravitational force terms instead of electromagnetic effects (Guo and Li 2008). This approach highlights the use of the Vlasov equation also on large spatial scales.

## 2.2 Solar plasmas

Plasma simulations of the solar system are mostly concerned with the modelling of solar activity and its influence on the heliosphere. Solar activity can be divided into two components: the solar wind, consisting of particles escaping continuously from the solar corona due to its thermal expansion; and transient phenomena such as flares and coronal mass ejections, during which energy and plasma are released explosively from the Sun into the heliosphere. Topics of active research in solar physics include for example the acceleration and expansion of the solar wind (for a review, see, e.g., Verscharen et al. 2019, and references therein), coronal heating (De Moortel and Browning 2015; Cranmer et al. 2017), details of solar dynamo and convection (Käpylä et al. 2023), and flux emergence (Schmieder et al. 2014). The latter is particularly important for transient solar activity, as flares and coronal mass ejections are due to the destabilisation of coronal magnetic structures through magnetic reconnection. Of great interest is also the propagation of the solar wind and solar transients into the heliosphere, in particular for studying their interaction with Earth and other planetary environments. Because of the large scales of the systems considered, solar and heliospheric simulations are generally based on MHD, indicating that the currently existing theories of the Sun and the solar eruptions concern thermalised plasmas. Applying the Vlasov approach to near-Earth physics, having important analogies to the solar plasmas, may therefore provide important for solar theories as well. Kinetic, non-thermal physics have been mentioned among the future aspirations in solar modelling (e.g., Gordovskyy et al. 2019)

## 2.3 Near-Earth space

Figure 1 illustrates the near-Earth space plasma system. Before interacting on Earth's magnetosphere, the supermagnetosonic solar wind is slowed down by the *bow shock*.
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Vlasiator 2D-3V modelling of the magnetosphere in the noon-midnight meridian plane viewed from the morning sector. The colour scale is illustrative only to depict the annotated different physical regions, which are outlined in the text. The solar wind and IMF directions are shown as arrows. The inset shows the coupling to the ionosphere domain in the form of field-aligned currents, from Ganse et al. (2025)</figcaption></figure>

(e.g. Balogh and Treumann 2013), and the region of shocked plasma downstream is the *magnetosheath* (Lucek et al. 2005). The interplanetary magnetic field (IMF) intensifies at the shock, increasing the magnetic field strength to roughly four-fold compared to that in the solar wind (e.g. Spreiter and Stahara 1994). The magnetosphere is the region dominated by the Earth's magnetic field (in green and purple in Fig. 1). It is compressed on the sunward side due to the pressure from the solar wind, and extends in a long magnetotail on the nightside. The magnetosphere is electromagnetically coupled to the conducting *ionosphere* (e.g. Palmroth et al. 2021).

The bow shock–magnetosheath system hosts highly variable and turbulent conditions. One of the most important factors controlling the level of variability in the magnetosheath is the angle between the bow shock normal and the IMF direction (e.g. Czaykowska et al. 2001; Dimmock 2017). At the side of the bow shock where the IMF and the shock normal are quasi-parallel (termed *quasi-parallel shock*), some particles reflect and propagate back upstream causing instabilities and waves in the *foreshock* that forms upstream of the bow shock (e.g. Hoppe et al. 1981; Eastwood et al. 2005; Wilson 2016; Turc et al. 2025). On the quasi-perpendicular side of the
bow shock, the downstream magnetosheath is less turbulent, but can exhibit large-scale waves originating from anisotropies in the ion distribution function (e.g. Génot et al. 2011; Soucek et al. 2015; Hoilijoki et al. 2016). The coupling of the foreshock–bow shock–magnetosheath system is under active research because magnetosheath plasma is ultimately the driver of the dynamics of near-Earth space. The most important open questions about the system include the processes that determine properties and evolution of the magnetosheath plasma in space and time. The system has been studied numerically with MHD models, which are useful for inferring average properties of plasma in the dayside system (e.g. Palmroth et al. 2001; Chapman and Cairns 2003; Dimmock and Nykyri 2013; Mejnertsen et al. 2018). However, they are unable to take into account plasma processes that require a non-Maxwellian ion distribution function, i.e., particle reflection, kinetic waves, and, therefore, the formation of the foreshock and the quasi-parallel and -perpendicular magnetosheath.

The earthward boundary of the magnetosheath is called *the magnetopause*, a current layer exhibiting large gradients in plasma parameters, which separates Earth’s magnetic field from the magnetosheath. Energy, momentum, and mass exchange between the upstream plasma and the magnetosphere occurs at the magnetopause (Palmroth et al. 2003, 2006; Pulkkinen et al. 2006; Anekallu et al. 2011; Daughton et al. 2014; Nakamura et al. 2017). Therefore, magnetopause processes are important in determining the amount of energy driving the space weather phenomena.

The magnetosphere is coupled to the ionised upper atmosphere called *the ionosphere* via the magnetic field lines. The Earth's ionosphere is a weakly ionised medium, divided into three regions – named *D* (60–90 km in altitude), *E* (90–150 km), and *F* (>150 km) – corresponding to three peaks in the electron density profile (Hargreaves 1995). From the magnetospheric point of view, the ionosphere represents a conducting layer closing currents flowing between the ionosphere and magnetosphere (Lysak 1990), reflecting waves (Wright and Russell 2014), and depositing precipitating particles (e.g., Rodger et al. 2013). Further, the ionosphere is a source of cold electrons (Cran-McGreehin and Wright 2005) and heavier ions (e.g., Peterson et al. 1981). These cold ions of ionospheric origin may affect local processes in the magnetosphere, such as magnetic reconnection at the magnetopause (André et al. 2010; Toledo-Redondo et al. 2016). Traditionally, global MHD models use an electrostatic module as basis for the ionosphere, coupled to the magnetosphere by currents, precipitation and electric potential (e.g., Goodman 1995; Janhunen et al. 2012). New paradigms on coupling the ionosphere to the magnetosphere have also been presented (e.g., Laundal et al. 2024).

The ionosphere itself is modelled either empirically or based on first principles: The International Reference Ionosphere (IRI) model describes the ionosphere empirically from 50 km to 1500 km altitude (Bilitza and Reinisch 2008; Bilitza et al. 2022), while for instance, the Sodankylä Ion and Neutral Chemistry model solves the photochemistry of the *D* region, taking into account several hundred chemical reactions (Verronen et al. 2005, and references therein). At higher altitudes, transport processes become important, and models such as the IRAP Plasmasphere–Ionosphere Model (IPIM; Marchaudon and Blelly 2015) couple a kinetic model for the transport of suprathermal electrons with a fluid approach to resolve the
chemistry and transport of ions and thermal electrons in the convecting ionosphere. Neither empirical nor first-principles based models are using the Vlasov equation for ionospheric physics, which at the ionosphere concerns much finer scales.

## 2.4 Other planetary environments

In general, the interaction of the solar wind with solid bodies is characterised by 1) whether the body has an intrinsic magnetic field, 2) whether it has an atmosphere, and 3) how strong the solar wind driving is at the distance of the body. While the modelling of the magnetospheres of the outer giants is only achievable to date using fluid approaches, the relatively weak dipole and the small size of Mercury's magnetosphere has been targeted for global PIC simulations (Richer et al. 2012; Lapenta et al. 2022; Ballerini et al. 2024). Kinetic models are also a popular tool for investigating the plasma environment of non-magnetised bodies such as Mars, Venus, comets, and asteroids. In the case of weakly magnetised bodies with atmospheres or cometary comas, the solar wind interacts directly with the atmosphere or coma, and especially the ionised component thereof. These highly variable interactions depend on the extent and density of the neutral populations and the ionisation processes generating the electromagnetically reactive ions.

At one extreme, Venus' dense atmosphere is a compact ionospheric obstacle, while on the other extreme, cometary outgassing is both a neutral coma as well as an extended source for ions (see Brain et al. 2016, Fig. 1 for an overview). These environments provide an ion component with gyration length scales comparable to the system size, prompting the use of multi-fluid and multi-species hybrid plasma models (e.g. Ma et al. 2002; Modolo et al. 2016; Jarvinen et al. 2020). These environments are rich with kinetic waves (via e.g. cyclotron resonance) and shocks (Szegö et al. 2000). Rubin et al. (2014) provide an illustrative example of the differences of multi-fluid and ion-kinetic models. These topics have so far not been extensively studied with Vlasov models.

For weakly magnetized and atmosphere-less bodies, the interaction on the inflow side is mainly absorptive for the plasma and diffusive for the magnetic field, but the body will leave behind a plasma wake. For instance, Moon's 3D wake structure has been extensively observed (Whang and Ness 1970; Zhang et al. 2014), characterising kinetic physics, such as counterstreaming ion and electron populations (Halekas et al. 2011). Vlasov simulations have been used to investigate electrostatic instabilities in the wake (Tao et al. 2012). In particular, Umeda and Ito (2014) and Umeda and Fukazawa (2015) have studied the interaction of a weakly magnetized body with the solar wind by means of full-Vlasov simulations.

Solar system bodies without global magnetic fields may still possess localized magnetic regions, such as the magnetic anomalies of Mars and the Moon. In the Martian plasma environment, these crustal remanent fields are a significant influence for the system, enabling e.g. reconnection between the IMF and the crustal magnetic fields (Hara et al. 2022), with some modelling work incorporating the crustal fields in the simulation setups (Ma et al. 2002; Romanelli et al. 2018).
<table><caption>**Table 1** Typical plasma parameters and scales for selected solar-terrestrial and astrophysical environments</caption><thead><tr><th>Physical system</th><th>Near-Earth space</th><th>Solar system</th><th>Interstellar gas</th><th>Neutron star</th></tr></thead><tbody><tr><td>λ<sub>D</sub> (m)</td><td>10<sup>-3</sup> – 10<sup>2</sup></td><td>10<sup>-4</sup> – 10<sup>1</sup></td><td>10</td><td>10<sup>-3</sup></td></tr><tr><td>c/ω<sub>pi</sub> (m)</td><td>&lt;10<sup>5</sup></td><td>10<sup>3</sup> – 10<sup>5</sup></td><td>10<sup>7</sup></td><td>0.1</td></tr><tr><td>r<sub>Li</sub> (m)</td><td>10<sup>3</sup> – 10<sup>6</sup></td><td>10<sup>1</sup> – 10<sup>7</sup></td><td>10<sup>8</sup></td><td>10<sup>-11</sup></td></tr><tr><td>System size (m)</td><td>10<sup>9</sup></td><td>10<sup>11</sup></td><td>1 – 10<sup>2</sup> pc</td><td>10<sup>4</sup></td></tr><tr><td>Process time scales</td><td>1 s – 1 d</td><td>1 s – 1 mo</td><td>300 kyr</td><td>1 s</td></tr></tbody></table>

λ<sub>D</sub>: the Debye length is the characteristic of a plasma related to its ability to shield out the electric potentials applied to it; c/ω<sub>pi</sub>: the ion inertial length is the scale at which ions decouple from electrons; r<sub>Li</sub>: the ion Larmor radius is the radius at which the ion gyrates around the magnetic field

## 2.5 Scales and processes

The following processes are central in explaining plasma behaviour in the solar-terrestrial system and astrophysical domains: 1) magnetic reconnection enabling energy and mass transfer between different magnetic domains, 2) shocks forming due to supersonic relative flow speeds between plasma populations, 3) turbulence providing energy dissipation across scales, and 4) plasma instabilities transferring energy between the plasma and waves. All these processes contribute to particle acceleration, which is one of the most researched topics within space and astrophysics, and notorious in requiring understanding from local microphysics to global scales. Below, we introduce some examples of these processes. Simulations of non-thermal space plasmas encompass a vast range of scales, from the smallest ones (electron kinetic scales, ion kinetic scales) to local and even global structures. Table 1 lists typical ranges of a handful of plasma parameters encountered in different branches of space physics and astrophysics. Especially in a larger astrophysical context, simulations cannot directly encompass all relevant spatial and temporal scales. However, it is important to note that scientific results of kinetic effects can be achieved even without directly resolving all the spatial scales that may at first glance appear to be a requirement (Pfau-Kempf et al. 2018).

### 2.5.1 Magnetic reconnection

Magnetic reconnection is a fundamental plasma process enabling the reconfiguration of magnetic topology, mixing of plasma populations, and the conversion of magnetic energy into bulk flow, thermal energy, and accelerated particles (Burch and Nakamura 2025). In collisionless plasmas, Vlasov theory provides a powerful framework for understanding reconnection dynamics (e.g., Shuster et al. 2021). Within the magnetosphere, reconnection occurs in abundance (e.g., Fuselier et al. 2024). Actively researched topics focus on understanding the nature and location of reconnection as a function of driving conditions (e.g. Hoilijoki et al. 2014; Fuselier et al. 2017; Han et al. 2024; Nakamura et al. 2025).
The Earth’s magnetotail hosts a hot and dense *plasma sheet*, the home of perhaps most diligent scientific investigations within the domain of magnetospheric physics, even leading to the conception of dedicated satellite constellations such as the 5-spacecraft THEMIS mission. Especially in focus have been explosive magnetotail eruptions that accelerate particles and involve both reconnection and kinetic instabilities (e.g. Sergeev et al. 2012; Palmroth et al. 2023). Tail disruptions have important analogues in solar eruptions (e.g. Birn and Hesse 2009; Birn et al. 2017), and investigating them with global simulations together with *in situ* measurements may shed light into other astrophysical systems as well.

In recent years, significant strides have been made largely due to high-resolution observations from the Magnetosphere Multi-Scale (MMS) mission in synergy with both local and global simulations. Despite this progress, many outstanding questions remain, as discussed in Nakamura et al. (2025). Among these, one of the most under-explored yet fundamental problems is reconnection onset. Recent observations include thin electron-scale current sheets (Wang et al. 2018), divergent electron velocity patterns without magnetic topology changes (Motoba et al. 2022), and electron-scale magnetic islands and Hall currents (Genestreti et al. 2023). Some works suggest that electron-only reconnection may represent an initial phase preceding the onset of standard ion-coupled reconnection in the magnetotail (Hubbert et al. 2022, see also Sect. 4.11). These findings highlight the complexity of reconnection and the necessity of combining *in situ* observations with simulations and theory.

## 2.5.2 Collisionless shocks

Shocks form due to plasma populations flowing supersonically with respect to each other, redistributing flow energy into thermal energy and accelerating particles (e.g. Balogh and Treumann 2013; Marcowith et al. 2016). Shock fronts such as those found at supernova explosions are efficient particle accelerators (Fermi 1949). Diffusive shock acceleration (e.g., Axford et al. 1977; Krymskii 1977; Blandford and Ostriker 1978; Bell 1978) is the primary source of solar energetic particles, and occurs from the non-relativistic (e.g. Lee 2005) to the hyper-relativistic (Aguilar et al. 2015) energy regimes. Shock-particle interactions including kinetic effects have been modelled using various analytical and semi-empirical methods (see, e.g. Afanasiev et al. 2015, 2018; Hu et al. 2017; Kozarev and Schwadron 2016; Le Roux and Arthur 2017), but drastic approximations are usually required in order to model the whole acceleration process, and Vlasov methods have not yet been utilised (Perri et al. 2022).

The classic extension of hydrodynamic shocks into the MHD regime has been disproven by a number of hybrid models due to, e.g., shock reformation (Caprioli and Spitkovsky 2013; Hao et al. 2017) and anisotropic pressure and energy imbalances due to non-thermal particle populations (Chao et al. 1995; Génot 2009). Only a self-consistent treatment including kinetic effects is capable of describing diffusive shock acceleration accurately. Recent works coupling shocks and high-energy particle effects include, e.g., those by Schwartz et al. (2022); Matsukiyo and Matsumoto (2024). Challenges associated with simulating shocks include modelling
gyrokinetic scales for ions whilst allowing the simulation to cover the large spatial volume involved in the particle trapping and energisation process. Radially expanding shock fronts within strong magnetic domains result in a requirement for high resolution both spatially and temporally. Modern numerical approaches usually make some compromises, e.g. reducing the problem to less than three dimensions either in ordinary or velocity space. The Vlasov approach is especially interesting in probing the physics of particle injection, trapping, acceleration and escape.

Numerical challenges are also caused by the presence of a foreshock or precursor upstream of a quasi-parallel shock (see also Sect. 2.3). Waves generated in the foreshock impact the shock front, cause shock reformation and affect particle acceleration (Wu et al. 2015; Haggerty and Caprioli 2020; Liu et al. 2021; Xu and Lazarian 2022). The waves can also transmit through the shock, influencing the plasma properties in the downstream. At Earth, they are observed to propagate into the magnetosphere (Turc et al. 2025). In addition to waves, transient mesoscale structures, with sizes ranging from a few ion gyroradii to several planetary radii, can emerge in the foreshock (Zhang et al. 2022; Kajdič et al. 2024). Observations indicate that these foreshock transients play an important role in particle acceleration near the shock, with both ions and electrons reaching relativistic energies (Wilson et al. 2016; Turner et al. 2018; Liu et al. 2019; Raptis et al. 2025). They may contribute to providing a seed population required for diffusive shock acceleration. Simulating foreshock processes, including foreshock transients, to understand their role in shock physics and their impact on the global dynamics of the system is therefore a challenging cross-scale problem.

### 2.5.3 Turbulence

Turbulence studies investigate how energy cascades and dissipates from large to small scales and vice versa. It is an outstanding question, which can be addressed using kinetic simulations (Bruno and Carbone 2013; Edyvean et al. 2024; Marino and Sorriso-Valvo 2023). Hybrid-Vlasov simulations (Valentini et al. 2010; Verscharen et al. 2012; Perrone et al. 2013; Arrò et al. 2024; Califano and Cerri 2023) have in particular been utilised to study the fluctuations around the ion inertial scale, which is of particular importance as it marks the transition between the fluid and kinetic scales. A recent breakthrough in solar wind turbulence is related to the so-called helicity barrier effect (Meyrand et al. 2021). Recent works have also focused on how turbulence transmits through shocks (Trotta et al. 2022, 2023) and how solar wind turbulence affects planetary plasma environments (Behar et al. 2022).

### 2.5.4 Instabilities

Plasma instabilities arise when a source of free energy in the plasma allows a wave mode to grow non-linearly. Plasma instabilities are ubiquitous in the universe. In Solar-Terrestrial physics, instabilities play an important role in, for example, the Kelvin-Helmholtz instability, which transfers solar wind plasma to the magnetosphere (e.g., Nakamura et al. 2017). In astrophysical media, instabilities are important for instance in accretion disks, where turbulence is driven by the
magnetorotational instability (e.g., Balbus and Hawley 1998). Vlasov models have been applied to many instabilities, such as the Rayleigh–Taylor instability (Umeda and Wada 2016, 2017), Weibel-type instabilities (Inglebert et al. 2011; Ghizzo et al. 2017), the Kelvin–Helmholtz instability (Umeda et al. 2010b; Settino et al. 2020; Tarvus et al. 2024), and kinetic current driven ion drift instability (Palmroth et al. 2023).

## 3 Modelling with the Vlasov equation

Various astrophysical and space environments such as those outlined above require a kinetic approach to describe plasma phenomena. In this framework, the plasma medium can be described as a system composed of particle populations. A comprehensive representation of the system implies that, for a given time instance, the spatial coordinate $\vec{x}$ and momentum $\vec{p}$ are defined for each individual particle. The numerical simulation of such a system involves solving the equation of motion for every single particle. When we encounter large-scale phenomena, the sheer abundance of particles – and, consequently, the number of equations to solve – becomes so immense that their direct numerical solution is impossible, even with numerical methods and the utilisation of the largest supercomputers available today. In some cases, when considering the macroscopic behaviour of a plasma system, a fully kinetic description may be redundant, and it becomes reasonable to treat plasma as a fluid using the MHD approximation, which is suitable for modelling e.g., large-scale transport and plasma average properties. However, the vast majority of non-stationary processes – such as reconnection, instabilities, and shocks – include non-Maxwellian plasmas, which are described using the Vlasov approach.

### 3.1 The Vlasov equation

In plasmas, as in neutral gases, the dynamical state of every constituent particle can be described by its position $\vec{x}$ and momentum $\vec{p}$ at a given time $t$. From a statistical perspective, the state of the system can be described by introducing a distribution function for species $s$ (electrons, protons, helium ions, etc) by $f_s = f(\vec{x}, \vec{p}, t)$, which represents the probability density in a six-dimensional phase space. The distribution function characterises the probability of finding a given number of particles of species $s$ within the spatial domain $d^3\vec{x}$ and having momentum within the range $d^3\vec{p}$. Hereinafter, we denote the species $s$ only if necessary.

By analogy with a thermal gas, the evolution of the distribution function in a plasma can be described by the Boltzmann equation:

$$ \frac{\partial f}{\partial t} + \vec{v} \cdot \frac{\partial f}{\partial \vec{x}} + \vec{F} \cdot \frac{\partial f}{\partial \vec{p}} = \left( \frac{\partial f}{\partial t} \right)_{\text{coll}} \quad (1) $$

where the terms represent the time evolution of the distribution function, the advection term, the external force $\vec{F}$ term, and particle interactions described by the collision integral $\left(\frac{\partial f}{\partial t}\right)_{\text{coll}}$. Collisions can be understood as particle-particle
interactions or, more broadly, as any interaction altering the momentum distribution. In highly luminous systems, encountered in astrophysics and laser-plasmas, an intense radiation field can also be a source for the collisions (see e.g., Nättilä 2024).

In majority of cases in space and astrophysical plasmas, where the plasma density is relatively low, collisions play a minor role. Instead, the behaviour of individual particles is predominantly governed by the collective fields generated by the ensemble of surrounding particles within a *sphere of influence*. This phenomenon, known as *collective behaviour*, allows us to neglect the collision integral. Assuming that external forces are described by the Lorentz force, the collisionless Boltzmann equation reduces to the Vlasov equation (Vlasov 1961):

$$ \frac{\partial f}{\partial t} + \vec{v} \cdot \frac{\partial f}{\partial \vec{x}} + q(\vec{E} + \vec{v} \times \vec{B}) \cdot \frac{\partial f}{\partial \vec{p}} = 0. \quad (2) $$

Here $q$ is the electric charge of the particle species considered. For small (coordinate) velocities, $|\vec{v}| \ll c$, where $c$ is the speed of light, the studied system is said to be Newtonian and the momentum reduces to $\vec{p} \rightarrow m\vec{v}$ (where $m$ is the mass of the particle species considered) and the distribution function simplifies to $f(\vec{x}, \vec{p}, t) \rightarrow f(\vec{x}, \vec{v}, t)$. If a significant part of the plasma acquires high enough kinetic energy, the relativistic effects start to become important: then, $\vec{p} = \gamma m \vec{v}$, where $\gamma = (1 - v^2/c^2)^{-1/2}$ is the Lorentz factor. The Vlasov equation (3) is Lorentz-invariant and, therefore, is generally valid also in the relativistic case (Thomas 2016). There are only very few numerical applications related to space or astrophysics that directly solve the relativistic Vlasov equation (see, e.g., Nättilä 2022), but it is more common in other contexts, such as laser-plasma interaction applications (Martins et al. 2010; Inglebert et al. 2011). By using a frame that is propagating at relativistic speeds (a Lorentz-boosted frame) the smallest time or space scales to be resolved become larger and the plasma length shrinks due to Lorentz contraction, indicating that the simulation execution times are accelerated.

## 3.2 Coupling with electromagnetic fields

The *self-consistent* simulation of a collisionless plasma environment requires coupling the Vlasov equation with the equations governing the force fields to form a closed system. The Vlasov equation describes the temporal evolution of the distribution function, using force fields – such as electromagnetic and gravitational fields – as input data. These fields are sourced by current, charge, and mass densities, all of which must be updated according to the evolving distribution function.

Depending on the spatio-temporal scales of the processes under study, an appropriate field model must be chosen. There are two primary strategies for closing the system of equations: the electrostatic framework, which enables the resolution of electron-scale dynamics while assuming ions to be immobile, and the electromagnetic framework, which is required for addressing cross-scale, complex problems. Additionally, in astrophysical contexts, the inclusion of the gravitational potential becomes essential when modelling large-scale structures such as galaxy formation.
### 3.2.1 The Vlasov–Maxwell system

A comprehensive description of a collisionless plasma where the gravitational potential is negligible is provided by solving the Vlasov equation

$$ \frac{\partial f_s}{\partial t} + \vec{v} \cdot \frac{\partial f_s}{\partial \vec{x}} + \frac{q_s}{m_s} (\vec{E} + \vec{v} \times \vec{B}) \cdot \frac{\partial f_s}{\partial \vec{v}} = 0 \quad (3) $$

for each particle species s. In Eq. (3), $\vec{E}$ and $\vec{B}$ are the electric and magnetic fields, respectively; $\vec{v}$ is the velocity coordinate, $q_s$ is the charge, and $m_s$ is the mass of the particle. Vlasov equation is coupled with the Maxwell equations for the electromagnetic fields. The sources for the fields are charge density $\rho_q$ obtained as the first-order moment of the distribution function for all species of electric charge $q_s$

$$ \rho_q = \sum_s q_s \int_{V_v} f_s(\vec{x}, \vec{v}, t) d^3\vec{v} \quad (4) $$

and current density $\vec{j}$

$$ \vec{j} = \sum_s q_s \int_{V_v} \vec{v} f_s(\vec{x}, \vec{v}, t) d^3 \vec{v}. \quad (5) $$

Maxwell's equations relate the fields to their sources, allowing us to determine the electric and magnetic fields. Specifically, the charge density contributes to the electric field $\vec{E}$ through Gauss's law:

$$ \nabla \cdot \vec{E} = \frac{\rho_q}{\epsilon_0}. \quad (6) $$

where $\epsilon_0$ is vacuum permittivity. The current density influences the magnetic field $\vec{B}$ according to Ampère's law:

$$ \nabla \times \vec{B} = \mu_0 \vec{j} + \mu_0 \epsilon_0 \frac{\partial \vec{E}}{\partial t}, \quad (7) $$

where $\mu_0$ is the vacuum permeability.

Faraday's law establishes the relationship between the electric and magnetic fields:

$$ \nabla \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}. \quad (8) $$

Numerical schemes typically discretise these three equations. The divergence-free condition ($\nabla \cdot \vec{B} = 0$) on the magnetic field is either enforced through a divergence-free solver or a numerical divergence-cleaning mechanism.

The Vlasov–Maxwell system is crucial for simulating a wide range of plasma phenomena, where electromagnetic fluctuations are driven by kinetic effects. These include the dynamics of charged particle beams, kinetic instabilities, magnetic reconnection, and turbulence in both space and astrophysical plasmas. The Vlasov–Maxwell system is also essential for understanding space weather phenomena, such
as solar flares and cosmic ray propagation, where both electromagnetic fields and the detailed shape of particle distributions must be considered simultaneously.

### 3.2.2 Hybrid-Vlasov systems

Hybrid-Vlasov systems retain solving the Vlasov equation for the non-Maxwellian species, utilising simplifications for other species that participate in the system dynamics. For example, an ion-kinetic hybrid model solves the Vlasov equation for ions (e.g. protons), while electrons can be modelled using a fluid approach, and ensuring charge neutrality. This has the advantage that the model is not required to resolve the short temporal and spatial scales in which electrons influence the electromagnetic fields. Conversely, an electron-kinetic hybrid-Vlasov simulation would solve the Vlasov equation for electrons, relegating ions to the role of a static background charge distribution. This approach is suitable, for instance, for evaluating electron-driven small-scale instabilities.

In hybrid-Vlasov systems, the most important quantities are the number density $n_s$, the bulk velocity $\vec{V}_s$, and the pressure tensor $P_s$, which are obtained as follows:

$$ 
\begin{aligned}
n_s &= \int_{V_v} f_s d^3 \vec{v}, \\
\vec{V}_s &= \frac{1}{n_s} \int_{V_v} \vec{v} f_s d^3 \vec{v}, \\
P_s &= \int_{V_v} m_s \vec{w} \otimes \vec{w} f_s d^3
\end{aligned} 
\quad (9) 
$$

where $\vec{w}$ is the thermal velocity $\vec{w} = \vec{v} - \vec{V}_s$, and $\otimes$ refers to a dyadic product.

Typically, an ion-kinetic hybrid-Vlasov system is closed by taking moments of the Vlasov equation and making approximations pertinent to the simulation system at hand. Integrating (3) over the velocity space, one gets the continuity equation or the zeroth moment of the Vlasov equation:

$$ \frac{\partial n_s}{\partial t} + \nabla \cdot (n_s \vec{V}_s) = 0, \quad (10) $$

Multiplying Eq. (3) by the phase-space momentum $m_s \vec{v}_s$ and integrating over the velocity space, one obtains the equation of motion or the first moment of the Vlasov equation:

$$ n_s m_s \left( \frac{\partial \vec{V}_s}{\partial t} + (\vec{V}_s \cdot \nabla) \vec{V}_s \right) - n_s q_s (\vec{E} + \vec{V}_s \times \vec{B}) + \nabla \cdot P_s = 0, \quad (11) $$

This process results in a hierarchy of moment equations, where each lower-order moment depends on the next higher-order one. To render the system tractable, a closure condition must be applied at a chosen level. In ion-kinetic hybrid-Vlasov models, a common approach is to close the system at the level of the first velocity moment by summing the electron and ion momentum equations. Terms proportional
to the small electron-to-ion mass ratio $m_e/m_i \ll 1$ are then neglected, leading to the generalised Ohm's law

$$ \vec{E} + \vec{V}_i \times \vec{B} = \frac{\vec{j}}{\sigma} + \frac{\vec{j} \times \vec{B}}{n_e e} - \frac{\nabla \cdot \mathcal{P}_e}{n_e e} + \frac{m_e}{n_e e^2} \frac{\partial \vec{j}}{\partial t}, \quad (12) $$

where $\vec{V}_i$ is the bulk velocity of the ions, $\sigma$ is the conductivity, $e$ is the elementary charge, $n_e$ is the electron number density, and $\vec{j}$ is given by Eq. (5). In the limit of slow temporal variations, the rightmost electron inertia term is typically neglected. Further, assuming high conductivity, one is left with the Hall term $\vec{j} \times \vec{B}/n_e e$ and the electron pressure gradient term $\nabla \cdot \mathcal{P}_e/n_e e$ on the right-hand side. The electron pressure term can be handled in a number of ways, such as using isobaric or isothermal assumptions. If electrons are assumed to be completely cold ($\mathcal{P}_e = 0$), the equation can be written as the Hall MHD Ohm's law

$$ \vec{E} + \vec{V}_i \times \vec{B} = \frac{\vec{j} \times \vec{B}}{n_e e}. \quad (13) $$

Thus, the ion-kinetic hybrid-Vlasov system of equations retains the Vlasov equation (3) for ions and the Maxwell–Ampère and Maxwell–Faraday equations (7) and (8) but replaces the Gauss equation by a generalised Ohm's equation (12) with appropriate approximations. If rapid field fluctuations are excluded from the solution, the displacement current from Ampère's law can be omitted, yielding

$$ \nabla \times \vec{B} = \mu_0 \vec{j}, \quad (14) $$

Equation (14) implies that the currents are divergence-free ($\nabla \cdot \vec{j} = 0$). The Darwin approximation makes the equation system more easily tractable, as then the Hall term becomes

$$ \frac{(\nabla \times \vec{B}) \times \vec{B}}{\mu_0 n_e e}. \quad (15) $$

In the case of purely transverse (divergence-free) electric fields, omitting the displacement current is equivalent to the *Darwin approximation* (see Sect. 4.5.3).

### 3.2.3 Unconventional closure

Depending on the physics to be included in the simulation, other, more exotic combinations of equations can be constructed to provide closure of the system. For example, the eVlasiator approach (Battarbee et al. 2021; Alho et al. 2022) models a subset of electron dynamics on top of a background of static ions and global static magnetic fields. In that modinkel, electron motion is solved through the Vlasov equation, but field updates are limited to electrostatic oscillations.
### 3.3 Coupling with gravity: The Vlasov-Poisson system

Considering the conservative force field $\vec{F} = -\vec{\nabla}U$, the Vlasov equation forms a closed system with the Poisson equation which defines the self-consistent collisionless gravitational potential $U$:

$$ \left\{ \begin{aligned} \frac{\partial f}{\partial t} + \vec{v} \cdot \frac{\partial f}{\partial \vec{x}} - \frac{\partial U}{\partial \vec{x}} \cdot \frac{\partial f}{\partial \vec{v}} = 0, & \nabla^2 U = 4\pi G \int_{V_v} f \, dv, \end{aligned} \right. \quad (16) $$

where $G$ is the gravitational constant. In an astrophysical context, this applies to the gravitational potential in a collisionless stellar system, such as galaxies, where stars interact via their collective gravitational field.

If the effects of magnetic field variations can be neglected and only phenomena on spatial scales much smaller than the characteristic ion scales are considered, it is sufficient to focus on the self-consistent evolution of electrostatic fields.

Thus, in non-relativistic systems, the Vlasov equation (3) takes the form:

$$ \frac{\partial f_s}{\partial t} + \vec{v} \cdot \frac{\partial f_s}{\partial \vec{x}} + \frac{q_s}{m_s} \vec{E} \cdot \frac{\partial f_s}{\partial \vec{v}} = 0 \quad (17) $$

for species $s$, and the system is closed by the Poisson equation:

$$ \nabla^2 \Phi + \frac{\rho_q}{\epsilon_0} = 0, \quad (18) $$

where $\Phi$ is the electric potential ($\vec{E} = -\vec{\nabla}\Phi$). The total charge density $\rho_q$ is given by Eq. (4).

The Vlasov-Poisson system is particularly applicable to problems where the dominant interactions arise from electrostatic forces, such as Langmuir wave dynamics and Landau damping, or to galactic dynamics problems where only the gravitational potential needs to be considered.

### 3.4 Fundamental properties of the Vlasov equation

The Vlasov equation has a number of useful properties helping to solve it numerically. Fundamentally, the Vlasov equation is a six-dimensional advection equation in phase-space, expressing that the total (material) derivative of phase-space density vanishes along the trajectories in phase-space. In the absence of any source terms, finding a solution at any given point in time requires to determine the motion of the phase-space density. One particularly convenient property follows from Liouville's theorem, from which the Vlasov equation is derived. It states that *phase-space density is constant along the trajectories of the system*. This means that a solution of the system at one point in time can be followed forward or backward in time arbitrarily, as long as the trajectories of phase-space elements, characterising the Vlasov equation, are known.

One consequence of Liouville's theorem is that initial density perturbations tend to form smaller and smaller structures as trajectories of the phase-space regions with different densities converge over time. In physical reality, this so-called filamentation
has a natural cutoff at scales where diffusive scattering effects become important, however, this is not part of the mathematical description of the Vlasov equation. Therefore, numerical implementations need to address this issue either by explicit filtering steps, or inherent numerical diffusivity.

A fundamental consideration in any physical modelling is the conservation of certain quantities, like mass, momentum, energy, and electric charge. This applies to Vlasov-based plasma modelling as well. Variational approaches have been used to derive the Vlasov-Maxwell and Vlasov-Poisson systems of equations as well as reduced forms (e.g. Marsden and Weinstein 1982; Ye and Morrison 1992; Brizard 2000; Brizard and Tronko 2011). Care has to be taken when developing numerical solutions of the Vlasov equation that quantities relevant to the problem to be solved are conserved adequately by the method (e.g. Filbet et al. 2001; Crouseilles et al. 2010; Cheng et al. 2013, 2014; Becerra-Sagredo et al. 2016).

# 4 Numerical modelling and HPC aspects

Finding solutions to the Vlasov equation in order to model physical systems typically involves simulations. Therefore, the phase space density $f(\vec{x}, \vec{v}, t)$ needs to be numerically represented, which strongly influences the choice of solvers and the resulting simulation code. This section is structured around different numerical representations of $f$.

A problem common to all numerical approaches solving the Vlasov equation is the curse of dimensionality – to fully reproduce all physical behaviour, the simulation domain must be 6-dimensional, with all 6 dimensions bg of sufficient resolution or fidelity for the desired physical system. This considerably impacts the size of the phase space and hence the computational burden of the algorithm.

## 4.1 Eulerian approach

In a straightforward manner, the velocity distribution function $f(\vec{x}, \vec{v}, t)$ can be discretised on a Eulerian grid, which can be operated on by different kinds of solvers (see Fig. 2a). The structure of the Vlasov equation allows both finite volume (FV) as well as semi-Lagrangian solvers to be employed, and all of these have been implemented with some success (Arber and Vann 2002). Discretisation of velocity space to a finite grid size $\Delta v$ also automatically imposes a lower limit for phase-space filamentation (see Sect. 3.4), at which the grid will naturally smooth out the structure. In some cases this is a purely numerically diffusive process, whereas others use explicit smoothing, filtering or subgrid modelling (e.g. Klimas 1987; Klimas and Farrell 1994). The 6-dimensional structure of the phase space, along with the physical scales and resolutions imposed by the underlying physical system (see Sect. 2) make a Eulerian representation in a Cartesian grid computationally impractical in the vast majority of cases concerning a large simulation volume.

Let us consider as an example a simulation of the entire Earth’s magnetosphere using a full 3D-3V Eulerian hybrid-Vlasov model. Let it extend up to the lunar orbit ($x \sim \pm 60 R_E$ in every direction) resolving approximately the solar wind ion inertial
length for bulk velocities ($\Delta x \sim 100$ km), and let the velocity space encompass typical solar wind velocities with some safety margin ($v \sim \pm 2000$ km/s) while resolving the solar wind thermal speed ($\Delta v \sim 30$ km/s). In this case, the resulting phase space would contain a total of $10^{18}$ cells. If each of them were represented by a single-precision floating point value, a minimum of 4 EiB of memory would be required to represent it!

Fortunately, there are many possibilities for simplification of the computational grid size:

* Reduction of phase-space dimension, if the physical system under consideration allows it, is an easy and efficient way to reduce computational load. Simulations of longitudinal wave instabilities (Jenab and Kourakis 2014; Shoucri 2008) and fundamental studies of filamentation have been performed in a 1D-1V setup, whereas laser wakefield and wave interaction simulations tend to be modelled in 2D-2V or 2D-3V setups (Besse et al. 2007; Sircombe et al. 2004; Thomas 2016). Another possibility here is to globally reduce the number of grid points by introducing a sparse grid representation, where the grid may be uneven with respect to Cartesian coordinates, while remaining static during runtime (Kormann and Sonnendrücker 2016; Guo and Cheng 2016). This is sometimes referred to as a sparse grid representation.

* Gyrokinetic simulations reduce the velocity space by dropping the azimuthal velocity dimensions perpendicular to the magnetic field, thus assuming complete gyrotropy of the distribution functions (e.g. Görler et al. 2011).

* Adaptively refined grids can be employed to reduce resolution and thus computational expense in areas of phase space that are considered less important for the physical problem at hand (Besse et al. 2008; Wettervik et al. 2017; Ganse et al. 2023). For mesh refinement to dynamically change during simulation run-time, an automatic classification and determination of required resolution in each part of the simulation domain is required. Depending on the physical process, formulating such a classification criterion requires a considerable amount of

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Different ways of numerically representing the phase-space density $f(\vec{x}, \vec{v}, t)$: a) In a Eulerian grid, every grid cell stores the local value of phase-space density, which is transported across cell boundaries. b) Spectral representations (shown here: Fourier-space in $\vec{v}$) allow for some update steps of phase-space density to be performed locally. c) In a tensor train representation, phase-space density is represented as a sum of tensor products of single coordinates' distribution functions which get transported individually</figcaption></figure>
physical insight. Kotipalo et al. (2024) compared spatial gradient-based and local kinetic physics measures for refinement of magnetospheric features.

* In many physical scenarios, large parts of phase space contain an extremely low, if not zero, density, so that their contribution to the overall dynamics of the system is negligible to the overall dynamic development. Suitable pruning of the phase-space grid can thus be performed to obtain a data structure that dynamically removes grid elements during runtime and keeps them only in regions deemed relevant for the physical system dynamics. The computational speed-up gained through the reduction of phase-space volume thus obtained can in some cases be a tradeoff against physical accuracy, and needs to be carefully considered. We have implemented this option in Vlasiator, and call it the dynamic sparse phase space representation, discussed more in Sect. 5.2.2. This method is not to be confused with the static sparse grid methods (Kormann and Sonnendrücker 2016; Guo and Cheng 2016) that are fundamentally dimension reduction techniques.

In space plasmas, where Lorentz force interactions are dominant over any collisional or gravitational terms, the magnetic field $\vec{B}$ makes the particles gyrate while the electric field $\vec{E}$ causes them to accelerate and drift. It can be advantageous to take these characteristics of acceleration into consideration when choosing an appropriate discretisation approach for the numerical phase-space representation. Common ideas include:

* A polar velocity coordinate system aligned with the magnetic field and centred around the drift velocity, $\vec{v} = (v_\|, v_r, v_\phi)$, in which the gyrophase coordinate $v_\phi$ has a much lower resolution than $v_\|$ and $v_r$. This can be employed in cases where the velocity distribution functions are known not to deviate strongly from gyrotropy, i.e., to exhibit cylindrical symmetry with respect to the magnetic field direction. However, the disadvantage of a polar velocity space over a Cartesian one is the more complex coordinate transformation required for transport into the spatial neighbours in any but the field-aligned direction.

* A Cartesian representation of velocity space, in which the coordinate axes co-rotate with the local gyration at every given spatial cell. Such a setup has the advantage that no transformation of velocity space due to the magnetic field will have to be performed at all, and no numerical diffusion due to gyration motion will occur. It does however come at the cost of more complicated spatial updates, since neighbouring spatial domains do no longer have collinear velocity space coordinate axes. A suitable interpolation and/or reconstruction has to take place when numerically solving the spatial transport problem, which again introduces numerical diffusivity problems potentially negating the advantages.

For the actual process of solving the Vlasov equation, a fundamental decision has to be made in the structure of the code, whether the solution steps are to be performed in a proper 6D manner (e.g. Vogman et al. 2018), or whether a Strang-splitting approach will be taken (Strang 1968; Cheng and Knorr 1976; Mangeney et al. 2002), in which the position and velocity space solution steps are performed independently of each
other. Due to the large number of involved dimensions, and thus computational context of each solver step, the latter approach tends to have significant performance benefits, whilst still achieving convergence (Einkemmer and Ostermann 2014). Alternative time-splitting methods based on Hamiltonians have also been proposed (e.g. Crouseilles et al. 2015; Casas et al. 2017).

If a Cartesian velocity grid is employed in the simulation, multiple families of solvers are available for it. Filbet and Sonnendrücker (2003). In all cases, the effects of diffusivity of the solver need to be considered. Especially uncontrolled velocity space diffusion manifests itself as numerical heating, as the distribution function tends to broaden over time. Higher orders of solvers and reconstruction methods, as well as explicit formulations in which moments of the distribution function are conserved are therefore advisable (Balsara 2017).

The choice of a Eulerian representation of phase space brings certain implementation details for HPC aspects with it. The relative ease of spatial decomposition into independent computational domains, which communicate through ghost cells, can be employed readily for Eulerian Vlasov simulations, providing a straightforward path to parallel implementations. On the other hand, the inherent limitations of Eulerian schemes (such as conditions for time steps, see Sect. 5.3.1) limit their overall numerical efficiency, and the high-dimensional nature of phase space can lead to challenges in appropriately represented and resolved numerical grids. As so often in HPC, design decisions have to be based on the specific properties of the physical system under investigation.

### 4.1.1 Finite volume solvers

As the Vlasov equation (3) is fundamentally a hyperbolic conservation law in 6D, it can be solved using the well-established methods of Finite Volumes (FV, Leveque 2002). In this approach, the phase-space fluxes are calculated through each interface of a discrete simulation volume (or phase-space cell) by reconstructing the phase-space density distribution through an appropriate interpolation scheme. The characteristic velocities at both sides of this interface are determined and the Riemann problem (Toro 2014) at each of these interfaces is solved to update the material content in each cell.

If Strang splitting is used to perform separate spatial and velocity-space updates, it is noteworthy that the state vector only contains a single scalar quantity (the phase-space density) and each cell interface update only needs to take a single characteristic velocity into consideration. For the update in a spatial direction, the characteristic is given by the corresponding cell's velocity space coordinates, whereas in the velocity space update step, the acceleration due to magnetic, electric and external field forces is homogeneous throughout each spatial cell. The Riemann problem for the Vlasov update does therefore not require the solution or approximation of an eigenvalue problem, which significantly simplifies its solution in comparison to hydrodynamic or MHD FV solvers. This property also enables the efficient use of semi-Lagrangian solvers (discussed more in Sect. 4.4).

As will be shown in Sect. 5, versions of the Vlasiator code until 2014 employed a FV formulation of its phase space update (von Alfthan et al. 2014) and numerous
other implementations exist (Banks and Hittinger 2010; Wettervik et al. 2017). A comprehensive introduction to the implementation and thorough analysis of the behaviour of a fully 6D implementation of a FV Vlasov solver is given by Vogman et al. (2018).

### 4.1.2 Finite difference solvers

While the Vlasov equation (3) could in principle be solved by directly employing finite difference methods, this approach does not seem to be favoured, and its applications in the literature appear to be limited to fundamental theory studies only (e.g., Schaeffer 1998; Holloway 1995). The largest issue with finite difference formulations is the lack of explicit conservation of the moments of the distribution function and related quantities. While high-order methods can still maintain suitable approximate conservation properties, their computational demands and/or diffusivity make them impractical.

## 4.2 Spectral solvers

Instead of a direct discretisation of the phase-space density $f(\vec{x}, \vec{v}, t)$ in its $\vec{x}$ and $\vec{v}$ coordinates, a change of basis functions can be employed, each coming with benefits and limitations. The transformation of $f$ into a different basis can be performed in the velocity coordinates only (cf. Fig. 2b), or in both spatial and velocity coordinates, depending on the physical application.

If a splitting scheme is employed, where velocity and real space advection updates are treated separately, the advection in a Fourier-transformed coordinate can be completely performed in Fourier space, as the transform of any coordinate $x \to k_x$ results in the differential advection operator $v_x \nabla_x$ turning into a simple multiplication:

$$v_x \nabla_x f(x) \xrightarrow{\text{Fourier}} i v_x k_x \tilde{f}(k_x). \qquad (19)$$

However, for the acceleration update, this transformation brings in the additional complication that the acceleration $\vec{a}$ would have to be independent of $\vec{v}$, which is true for the electrostatic Vlasov-Poisson system, but incorrect in Vlasov-Maxwell scenarios, due to the $v$-dependence of the Lorentz force. In order to accommodate velocity-dependent acceleration, solving a system in such a way typically requires multiple forward and backward Fourier transforms within one time step (Klimas and Farrell 1994).

In a thus-represented velocity space, the question arises of which maximum velocity space frequency $\vec{k}_{v,\text{max}}$ is available, and the filamentation problem itself becomes a boundary problem at the maximum extents of velocity $\vec{k}$-space (Eliasson 2011). However, stability issues of this scheme remain under discussion (Figua et al. 2000; Klimas et al. 2017).
Finally, a full Fourier-space representation of $\tilde{f}(\vec{k}_r, \vec{k}_v, t)$, in which also the real space coordinate $\vec{x} \to \vec{k}_r$ is transformed is a possibility. However, it further complicates the treatment of configuration and velocity space boundaries (Eliasson 2001). When used with periodic spatial boundary conditions, such a setup can be quite efficient for the study of kinetic plasma wave interactions.

Apart from the Fourier basis, other orthogonal function systems can be used as the basis for description of phase-space densities. A popular choice is presented by Hermite functions (Koshkarov et al. 2021; Delzanno 2015; Camporeale et al. 2016), whose $L^2$ convergence behaviour closely matches that of physical velocity distribution functions, and whose property of being eigenfunctions to the Fourier transform can be used as a numerical advantage. Since the zeroth Hermite function $H_0(\vec{v})$ is simply a Maxwellian particle distribution, a hybrid-Vlasov code with Hermitian basis functions should replicate MHD behaviour in this limit. Adaptive inclusion of higher-order Hermite functions then allows an increasing amount of kinetic physics to be numerically represented (Vencels et al. 2015).

A common problem of any kind of spectral method, be it Fourier-based, or using any other choice of nonlocalised basis functions, is the formulation of boundary conditions. While simulations of wave or scattering behaviour can usually assume periodic boundary conditions, macroscopic systems will require modelling boundaries where the plasma interacts with solid or gaseous matter (such as planetary or stellar surfaces), inflow conditions as well as outflow boundaries. Due to the unavailability of suitable spectral formulations for these boundaries, spectral-domain solvers have not gained foothold in modelling efforts of macroscopic astrophysical systems.

In any nonlocal choice of basis function for the phase-space representation, be it Fourier-, Hermite- or wavelet-based (Besse et al. 2008), extra thought has to be put into scalability of parallel solvers for it. If a change of basis function (such as a switch from a real-space to a Fourier space representation) is required as part of the simulation update step, this will typically not scale beyond hundreds of cores in supercomputing environments. In a recent advance, Koshkarov et al. (2021) presented a method based on a Hermite basis for velocity space discretisation that employs a direct Cartesian real space formulation with a discontinuous Galerkin solver approach and displays excellent scaling properties on HPC hardware.

## 4.3 Tensor train

An entirely separate class of numerical representations for the phase-space density is provided by the tensor train formalism (Kormann 2015) illustrated in Fig. 2c. The idea behind this approach is inspired by Strang-splitting solvers, in which spatial and velocity dimensions are treated in individual and consecutive solver steps. The distribution function $f(x_1, x_2, \dots, x_n)$ is represented as a tensor product of component basis functions,
$$f(x_1, x_2, \dots, x_n) = \prod_{k=1}^{n} f_k(x_k) \qquad (20)$$

in which each $f_k(x_k)$ is only dependent on a single coordinate $x_k$, and thus only affected by a single dimension’s update step. The generalised formulation is called the *Tensor Train* of ranks $r_1 \dots r_n$ (see Fig. 2c),

$$f(x_1, x_2, \dots, x_n) = \sum_{\alpha_1=1}^{r_1} \dots \sum_{\alpha_n=1}^{r_n} \prod_{k=1}^{n} Q_k(\alpha_k - 1, x_k, \alpha_k), \qquad (21)$$

in which the distribution function is entirely formulated in terms of sums of products of $Q_k$, which themselves only depend on a single coordinate $x_k$. Transport can be performed by individually affecting each $Q_k$, followed by a rounding step to keep the tensor train compact.

While this approach has so far only been employed in low-dimensional approaches and for feasibility studies, it has gained momentum through the extensive use of tensor train formulations in machine learning applications (where the same concept is alternatively referred to as “Matrix product states” or “Low Rank Approximations”, LoRA) (Oseledets 2011; Orús 2014) and attempts at formulating and optimizing Vlasov solver implementations are an active area of research (Ehrlacher and Lombardi 2017; Guo and Qiu 2022). Similarly, dynamical low-rank formulations of the Vlasov and Boltzmann equations promise a mathematically similar and numerically highly practical approach (Einkemmer and Joseph 2021; Einkemmer and Lubich 2018). Large numerical simulations using tensor train models similar low-rank approximations in space physics have not yet been performed, efforts to integrate them into existing codebases are ongoing.

## 4.4 Semi-Lagrangian and fully Lagrangian solvers

As a consequence of Liouville’s theorem (cf. Sect. 3.4), numerical solutions of the Vlasov equation can be elegantly formulated in Lagrangian and semi-Lagrangian ways, by following the characteristics in phase space. Since the spatial velocity of any point in phase space is simply given by its velocity space coordinates, and its acceleration due to Lorentz’ force is provided by the local electromagnetic field quantities, a unique characteristic for each point in phase space is easily obtained in a simulation (cf. Sect. 4.1.1).

As the simulation progresses, the distribution of these sample points will shift, maintaining their initial phase-space density values, and the volumes in between them obtain phase-space density values through interpolation. If necessary, new samples can be created, or existing ones merged, where filamentation requires it. In essence, fully Lagrangian simulation codes (Kazeminezhad et al. 2003; Nunn 2005; Jenab and Kourakis 2014) track the motion of samples of density through phase space, stepping forward in time, resulting in an updated phase-space distribution. This is illustrated in Fig. 3a. Sometimes, these methods are referred to as Lagrangian particle methods, as each phase-space sample can be modelled as a macroparticle.
Much more common than the fully Lagrangian formulation of Vlasov solvers is the family of semi-Lagrangian solvers (Sonnendrücker et al. 1999). In these, the phase-space samples to be propagated are obtained at every time step from a Eulerian description of phase space, their transport along the characteristics is calculated within the time step, and the resulting updated phase-space density is sampled back into a Eulerian grid (which can be either structured or unstructured, see Besse and Sonnendrücker 2003). This process can be performed either forwards in time (Crouseilles et al. 2009, see Fig. 3b), in which the source grid points are scattered into the target locations, or backwards in time (Sonnendrücker et al. 1999; Pfau-Kempf 2016), where each target grid point performs a gather operation, spatially interpolating inside the previous state of the time steps (Fig. 3c). Backwards semi-Lagrangian methods are sometimes also referred to as Flux Balance Methods (see Filbet et al. 2001). Either way, an interpolation step will be involved which may again lead to significant numerical diffusion, unless methods are used to minimise it. Some of the more common interpolation procedures used are cubic splines and Hermite reconstruction because they produce smooth results with reasonable accuracy and are less dissipative than other methods using continuous interpolations (Filbet and Sonnendrücker 2003). Lagrange interpolation methods produce more accurate results but require higher-order polynomials and large stencils to limit diffusion. The high-order discontinuous Galerkin method for spatial discretisation, along with a semi-Lagrangian time stepping method, has also been used in Vlasov-Poisson systems providing an improvement in accuracy compared to previously used techniques (Rossmanith and Seal 2011). The flexibility of combining different approaches is also seen in a recent particle-based semi-Lagrangian method for solving the Vlasov-Poisson equation (Cottet 2018).

Cheng and Knorr (1976) were the first authors to employ semi-Lagrangian updates of a Vlasov-Poisson problem in a Strang-splitting setup, which they refer to as the time-splitting scheme, in which they independently treated advection due to temporal and spatial updates. Mangeney et al. (2002) later formulated a Strang-splitting scheme for the Vlasov-Maxwell equation. As such a splitting

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Illustration of Lagrangian and semi-Lagrangian approaches: a) In a full Lagrangian solver, phase-space samples are moved along their characteristic trajectories, never mapped to a grid and the phase-space values in between are obtained by interpolation. b) In a forward semi-Lagrangian method, the same process is used, but the phase-space values are redistributed to a grid at regular intervals. c) A backward semi-Lagrangian method follows the characteristic of every target phase-space point backwards in time and obtains the source value by interpolation</figcaption></figure>
scheme performs acceleration and translation steps separately, the phase-space trajectories of any simulation point approximates their physical behaviour in a staircase-like dimension-by-dimension pattern. Further, Fatone et al. (2019) investigated the application of a semi-Lagrangian spectral method to solve the Vlasov-Poisson system, employing Legendre polynomials and Hermite functions for approximating the velocity distribution function. The method exhibits good conservation properties, validated through numerical tests on standard benchmark simulations like the two-stream instability and Landau damping.

## 4.5 Field solvers

The Vlasov equation does not stand alone in describing the physical system in consideration, but requires a further prescription of the fields introducing the force terms. In the vast majority of cases in computational astrophysics, these will be electromagnetic forces self-consistently produced through the motion of charged particles in plasma, although there have been examples of Vlasov-gravity simulations (Guo and Li 2008), in which the Poisson equation was solved based on the simulation's mass distribution. Also in a few cases, the fields affecting phase-space distributions are considered entirely an external simulation input, with no feedback from the phase-space density onto the fields (Palmroth et al. 2013; Pfau-Kempf et al. 2018), which can be called “test-Vlasov” in analogy to test-particle simulations. These are particularly useful as test cases before the fully operational code can be run.

A key requirement for any field solver is to preserve the solenoidality of the magnetic field $\vec{B}$ expressed by equation (6). There are two main avenues used to achieve this goal (e.g., Tóth 2000; Balsara and Kim 2004; Zhang and Feng 2016, and references therein). Either the field reconstruction is divergence-free by design, such as the one used in Vlasiator (see Sect. 5.3), or a procedure is needed to periodically clean up the divergence of $\vec{B}$ arising from numerical errors.

In the following sections, different solvers for electromagnetic fields (and their simplifications) will be discussed in relation to Vlasov simulations. These are fundamentally very similar in structure to the field solvers used in other simulation methods, such as PIC and MHD, and can in many cases be adapted directly from these with little change.

### 4.5.1 Electrostatic solvers

If modelling a physical system in which the interaction of plasma with magnetic fields is of little importance (such as electrostatic wave instabilities, dusty plasmas, surface interactions Chane-Yook et al. 2006 and other typically local phenomena), the magnetic field force ($q\vec{v} \times \vec{B}$) part of the Vlasov equation can be neglected, and a purely electrostatic system remains.

Neglecting the effects of $\vec{B}$ completely leads to the Vlasov-Poisson system of equations (17) and (18), for which the field solver needs to find a solution to the Poisson equation at every time step. Being an elliptic differential equation that is
solved instantaneously in time, no time step limit arises from the field solver itself. Typically, solvers use approximate iterative approaches, multigrid methods or Fourier-space solutions (Birdsall and Langdon 2004).

Another option, if an initial solution for the electric field has been found (or happens to be trivial), is to update it in time by using Ampère's equation in the absence of $\vec{B}$,

$$ \frac{\partial \vec{E}}{\partial t} = -\epsilon_0 \vec{J} \qquad (22) $$

in either an explicit finite-difference manner, or using more advanced implicit formulations (Cheng et al. 2014). However, special care should be taken to prevent the violation of the Gauss law (6) by using appropriate numerical methods.

### 4.5.2 Full electromagnetic solvers

If the kinetic physics of both electrons and ions is to be considered, and particularly if radio wave or synchrotron emissions are intended outcomes of the system, one must use the full set of Maxwell's equations. A popular and well-established family of electromagnetic field solvers is the finite difference time-domain (FDTD) approach, which has a longstanding history in electrical engineering applications. In formulating the finite differences for the $\partial \vec{E}/\partial t \sim \nabla \times \vec{B}$ and $\partial \vec{B}/\partial t \sim \nabla \times \vec{E}$ terms of Maxwell's equations, it is often advantageous to use a staggered-grid approach, in which the electric field and magnetic field grids are offset from one another by half a grid spacing in every direction (Yee 1966). In this setup, every component of the electric field vector is surrounded by magnetic field components and vice versa, so that the finite difference evaluation of the rotation can be performed without any need for interpolation.

Care should be taken when employing FDTD solvers for studies of wave propagation at high frequencies or wave numbers, as the numerical dispersion relations of waves are deviating from their physical counterparts for high $k$, and this effect in particular is anisotropic in nature, and most strongly pronounced in cases of diagonal propagation, due to the intrinsic differences in the manner by which grid-aligned and non-grid-aligned computations are handled. Kärkkäinen et al. (2006) and Vay et al. (2011) present a thorough analysis of this problem in the case of PIC simulations, and provide suggestions for mitigating their effects. The largest disadvantage of FDTD solvers is their stringent requirement to resolve the propagation of fields at the speed of light, thus leading to extremely short time step lengths. In order to simulate anything at non-microscopic timescales, other methods will need to be used. Fourier-space solvers of Maxwell's equations are advantageous in this respect, as they do not come with fundamental time step limitations. This is weighed up by the fact that their parallelisation is more difficult, and the formulation of appropriate boundary conditions is not always possible (cf. Sect. 4.2).

The solvers discussed so far can be used to solve the kinetic dynamics of both electrons and ions. Examples of such fully kinetic codes are Pezzi et al. (2019a);
Hakim and Juno (2020). They are bound to resolve electron-kinetic length scales and time scales determined by their light crossing time.

### 4.5.3 Hybrid solvers

If large-scale phenomena with timescales much larger than the local light crossing time are being investigated, FDTD Maxwell solvers quickly lose their appeal. If magnetic field phenomena are still to be considered self-consistently in the simulation, appropriate modifications of the electrodynamic behaviour have to be taken, so that modelling these phenomena with longer time steps becomes feasible.

One common way to get rid of the speed of light as a limiting factor is by getting rid of the electromagnetic mode as a solution to Maxwell's equation altogether, in a process called the *Darwin approximation* (Schmitz and Grauer 2006a; Bauer and Kunze 2005) for systems where displacement currents and relativistic effects can be neglected. In this process, the electric field is decomposed into its longitudinal and transverse components $\vec{E} = \vec{E}_L + \vec{E}_T$, with $\nabla \times \vec{E}_L = 0$ and $\nabla \cdot \vec{E}_T = 0$. Only $E_T$ is allowed to participate in the temporal update of $\vec{B}$, which decouples the Maxwell–Ampère and Maxwell–Faraday equations and rules the electromagnetic mode out of the simulated physical system. As a result, the fastest remaining wave mode in the system becomes the Alfvén wave, and the maximum time step rises significantly, by a factor of $c/v_A$.

For closure of the Maxwell equations, moments of the distribution functions of all species, including electrons, is required. In the absence of a full kinetic treatment of electrons, the moment closure is achieved through approximative terms in the modified Ohm's law (Eq. 12), providing a relation between the electric and magnetic field. The level of complexity of Ohm's law directly influences the simulation results as it immediately affects the kinetic physics described by the model.

## 4.6 Other solvers

There are some emerging approaches to solving the Vlasov equation for plasmas. Engel et al. (2019) designed a quantum algorithm for the Vlasov–Maxwell system at the limit of electrostatic Landau damping. In principle, such quantum algorithms could provide exponential speed-ups, but Engel et al. (2019) estimate measurement errors remaining a hindrance for extracting accurate results, and a general solution for the Vlasov–Maxwell system will require considerable effort.

As machine learning progresses and data availability increases, researchers are investigating the potential of neural networks in handling complex systems governed by high-dimensional partial differential equations (PDEs). Within the realm of scientific machine learning, physics-informed neural networks (PINNs, Raissi et al. 2019) are leading this effort, utilising the capabilities of deep learning to restrict the predictions of the underlying physical laws and invariances in the system. In a recent study, Zhang et al. (2023) demonstrated the effectiveness of a PINN implemented by a multi-layer perceptron (MLP) to model solutions for the Vlasov–Poisson system. They evaluated solution accuracy by simulating the two-stream instability using a
PIC code and training the PINN with randomly sampled data in time, spatial and velocity space. By employing a loss function which encoded the physical constraints of the Vlasov-Poisson system, they assessed the PINN's accuracy by comparing the reconstructed electron distribution function to the PIC solution, showing errors within the range of 3% - 5%. We foresee that these efforts will be greatly emphasised in the future.

## 4.7 Boundary conditions

Boundary conditions are important in the stability of simulations of physical systems. For example, problems that involve wave propagation (such as those dictated by the Vlasov equation) it is important that waves are able to leave the simulation domain without generating spurious reflections at the simulation walls. While common boundary conditions like Dirichlet (prescribed values at the boundary) and Neumann (prescribed derivatives at the boundary) are often employed, they do not always work well for open boundaries.

A popular technique for treating open boundary conditions in numerical wave propagation simulations is the implementation of a Perfectly Matched Layer (PML, Berenger 1994). This method is widely employed in time domain solvers such as computational electromagnetics (PML, Sofi et al. 2021) and acoustic wave propagation (PML, Li 2010), and in PIC simulations (e.g., Pasik et al. 1999). The PML serves as an artificial absorbing boundary condition (ABC) that enables the transmission of outgoing waves (at any incident angle and frequency) out of the computational domain while suppressing spurious reflections at the boundaries. PMLs are implemented by extending the domain with either a complex coordinate stretching or by adding an anisotropic absorbing layer normal to the boundary that ensures impedance matching at the interface with the computational domain, minimising reflections.

As such, the PML may seem like an ideal candidate for an ABC for Vlasov simulations. However, the PML boundary condition is derived under the assumption that the wave equation in the absorbing layer matches the wave equation in the simulation domain. In homogeneous, isotropic media, this match can be made exact and allow for reflection-less dumping of the outgoing waves. However, in Vlasov-based plasma, which is by definition anisotropic and inhomogeneous (permittivity and permeability are not constant), this matching cannot be achieved at the interface of the PML and the numerical domain. This leads to side effects such as unstable growth of other wave modes, reflections, and artifacts generated at the interface of the two domains. rendering the simulation invalid.

## 4.8 Coupling schemes

A reduction of the computational burden of a model can be achieved by coupling different schemes in order to focus the use of the costlier kinetic model on the region (s) of interest while solving other parts of the system with less intensive algorithms. This is also a means to extend the simulation domain where one system is taken as
the boundary condition of the other. Various classes of coupled models exist, depending on the coupling interface chosen.

One strategy is to define a spatial region of interest in which the expensive kinetic model is applied, embedded in a wider domain covered by a significantly more computation-efficient fluid model. Rieke et al. (2015) present a method for modelling electrons and protons by a full Vlasov method in parts of the simulation domain, while retaining to a 2-fluid code in areas of less interest. The multiscale model utilising a hierarchy of Vlasov models combined with a multi-moment fluids description was recently implemented in the Murphy2 code (Allmann-Rahn et al. 2024). There are also approaches coupling a fluid code to a PIC code both using a static and adaptive mesh refinement (e.g. Tóth et al. 2016; Chen et al. 2017; Li et al. 2024), used also in reconnection studies (e.g. Usami et al. 2013; Zhou et al. 2020). The disadvantage of this strategy is that scale coupling cannot be addressed as the kinetic effects do not spread into the fluid regime, and smaller-scale physics can only affect the solution in the domain at which the kinetic physics is in force.

Another strategy consists in defining the regions of interest in velocity space, that is coupling a fluid scheme describing the large-scale behaviour of a system with a Vlasov model handling suprathermal populations introducing kinetics into the model. Again, this is a development for which a certain amount of theoretical work and testing on small cases has been done (e.g. Tronci et al. 2014) but not yet extended to larger-scale applications.

One feature of discretised Vlasov schemes is that positivity and energy conservation may be broken at the numerical level. In Allmann-Rahn et al. (2022), a Vlasov solver was coupled to a multi-moment fluid method, wherein a portion of the velocity distribution function was corrected at each step based on the multi-moment approach, ensuring energy conservation. This approach was verified using test-cases of Landau damping, two-stream instability, and reconnection based on the Geospace Environmental Modeling (GEM) Reconnection Challenge (Birn et al. 2001).

## 4.9 Subgrid models

Especially when modelling physics with large scale separations, it may be unfeasible to use sufficient resolution to capture all fine-scale details of the problem at hand (Dubart et al. 2020). Subgrid models can produce approximate effects of unresolved physics, and they are widely used in atmospheric sciences in modelling of turbulence to account for dissipation effects at small scales (Deardorff 1985). Similarly, subgrid models for turbulence in magnetised plasmas have been investigated (e.g. Teaca et al. 2021). Lately, Miura and Hamba (2022) have implemented a Hall MHD subgrid model effectively doubling their resolution with respect to plasma turbulence. Further, effects of small-scale instabilities like the tearing mode can be described as hyper-resistivity, or anomalous electron viscosity (Strauss 1986). This may have special relevance in studies of magnetic reconnection, where the details of such dissipative processes that can be modelled as hyper-resistive terms may influence the behaviour of the reconnecting systems (Aunai et al. 2013). Similar hyper-resistive
terms can also be used in plasma simulations to introduce stabilizing numerical diffusion (Maron et al. 2008; Aunai et al. 2024).

In addition to modelling dissipation of energy in fluids or fields, numerical simulations can include the effects of other under-resolved or completely unresolved physical phenomena through subgrid models. Solving the Vlasov equation can benefit from sub-grid modelling of the behaviour of the VDFs, such as the pitch-angle diffusion model developed by Dubart et al. (2022). This model alleviates the lack of ElectroMagnetic Ion Cyclotron (EMIC) wavemodes and mirror modes in 3D-3V Vlasiator simulations where spatial resolution is not sufficient to properly model these plasma instabilities (see Sect. 5.7). In effect, it adds an extra term to the Vlasov equation (3), which can act akin to a collisional term as seen in equation (1). A case study of implementing and testing a subgrid model to counteract under-resolving these features is shown in Dubart et al. (2023), indicating that the sub-grid model is capable of sufficiently correcting temperature anisotropy evolution in magnetosheath plasma, resulting in a more accurate pressure balance in the global model.

## 4.10 Computational burden of Vlasov simulations

Representing numerically the complete phase space of a kinetic plasma system including all required physical processes is computationally intensive, and a large amount of data needs to be stored and processed. Different possible representations of the phase-space distribution and solution methods and their expected scaling shall be given in this section. Computational requirements for equivalent PIC simulations are estimated in comparison, although due to their different tuneable parameters, a rigorous comparison is difficult and beyond the scope of this work.

As shown in Sect. 4.1, a blunt Eulerian discretisation of a magnetospheric simulation without any velocity space sparsity results in $10^{18}$ sample points or a minimum of 4 EiB memory requirement, which is unrealistic on current and next-generation architectures. A first approach is to reduce the dimensionality from a full 3D space to a 2D slice, which results in a reduction of sample points of the order of $10^4$. Obviously a further reduction to 1D yields a similar gain. With a sparse velocity space strategy as used in Vlasiator (see Sect. 5.2.2 below) a further reduction by a factor of $10^2 - 10^3$ sample points can be achieved. Typically modern large-scale kinetic simulations both with Vlasov-based methods (e.g. Palmroth et al. 2023) and PIC methods (e.g. Ren et al. 2024) reach an order of magnitude of $10^{11} - 10^{12}$ sample points.

Beyond the number of sample points to be treated, the length of the propagation time step relative to the total simulation time aimed for is a crucial component of the computational burden of a model. Certain classes of solvers are limited in that respect as they cannot allow a signal to propagate more than one sampling interval or discretisation cell per time step (see Sect. 5.3.1). With respect to hybrid models using the Darwin approximation, the inclusion of electromagnetic (light) waves in the model description results in a reduction of the allowable time step by a factor of $10^3$ or more. Eulerian solvers typically have similar limitations which can impact the time step by a factor of approximately $10^2$ due to the Larmor motion in velocity space in
the presence of a magnetic field. Subcycling strategies and the use of Lagrangian algorithms are common approaches to alleviate these issues, at the potential cost of some physical detail, however.

A choice of basis function for the representation of velocity space other than Eulerian grids (like spectral or Hermite bases) can in many cases be beneficial to limit the memory requirements for reasonable approximations of the velocity space morphology. Care must however be taken that non-local transformations from one basis to another, such as a Fourier transform, tend to have unfavourable scaling behaviour in massively parallel implementations. Tensor-train formulations appear to be a promising avenue for the representation of phase space densities that have suitable computational properties, but large-scale space plasma applications have not been demonstrated yet.

Higher computational requirements are expected if physics of multiple particle species (especially electrons) are essential for the system under investigation. The need to represent multiple separate distribution functions multiplies the memory and computation requirements. The relative mass ratios of these species also have an effect on the kinetic time and length scales that need to be resolved. Going from a purely ion-kinetic hybrid-Vlasov to a full-Vlasov simulation, in which electrons are included as a kinetic species shortens the Larmor times by a factor of $m_p/m_e = 1836$ and depending on the employed solver may require resolution of the plasma's Debye length. The latter factor means that, with respect to the reference hybrid simulation considered above, which approximately resolves the ion kinetic scales, a spatial resolution increase of the rough order of $10^5$ would be required (see Table 1), amounting to a $10^{15}$ more sampling points. In order to reduce this considerable overhead, a common approach is to rescale physical quantities such as the electron-to-proton mass ratio and/or the speed of light (e.g. Hockney and Eastwood 1988), while aiming to preserve the essential features of kinetic plasma behaviour.

Most of these scaling relations apply likewise in PIC. In these, however, the parameter most strongly affecting the computational burden of the phase-space representation is the particle count. As a rule of thumb, a PIC simulation with a particle count similar to the sampling point count of an equivalent Vlasov simulation will have a similar overall computational cost. For many physical scenarios, this particle count can be chosen to be significantly lower (on the order of 100-1000 particles/cell), especially if noisy representations of the velocity spaces are acceptable. In simulations with high dynamic density gradients with depleted regions, and in setups where a minimisation of sampling noise is essential (such as investigations of nonlinear wave phenomena), PIC and Vlasov simulations are expected to reach a break-even point in terms of computational burden.

Even if the simulated system is to be re-scaled, simplified, or under-resolved, there is a possibility that the number of sampling points can reach very large values. To efficiently model such a system for hundreds of thousands of time steps, a distributed computing system is likely required. In order to ensure efficient use of computing resources in such a system, profiling of the simulation code is of utmost importance, pinpointing regions where optimisations can be performed. For such tasks,
simulation codes can be coupled with performance analysis toolkits (e.g., Coti et al. 2024).

### 4.10.1 Parallelisation on accelerated hardware

As solving the Vlasov equation requires computing phase-space updates for a very large amount of grid points, parallelisation of the computational solution on accelerated hardware is a very enticing scenario. Porting simulation codes to heterogeneous and/or parallel architectures such as computing clusters with general purpose graphics processing units (GPGPUs, or GPUs for short) (Battarbee et al. 2025) can improve performance through both more computations being performed in parallel, but also through access to high-bandwidth memory (HBM).

Vlasov methods often involve a grid of phase-space values the geometry of which is known in advance, which should lead to quasi-trivial parallelisation on accelerated hardware. Early hardware implementations had very low double-precision floating point (FP64) performance, steering implementations into less accurate single-precision (FP32) calculations. Latu (2011) evaluated this accuracy trade-off for an early 1D-1V local spline approach, with a more recent study on the topic being of the semi-Lagrangian discontinuous Galerkin scheme by Einkemmer (2020). Other early adoptions of GPGPUs included Rocha Filho (2013) (1D-1V, including long-range interactions), Mehrenberger et al. (2013) (1D-1V, reformulating the semi-Lagrangian approach to utilise a circulant matrix method with FFTs), and Sandroos et al. (2013) (3D-3V finite-volume approach). All these early adoptions were written directly in CUDA for NVIDIA accelerator hardware.

As availability and adoption of accelerated hardware at computing clusters increased, more codes were written or ported with these architectures in mind. Many still utilised CUDA for efficient but vendor-locked code, such as Malkov and Kudryavtsev (2020) (2D-2V, conservative fifth-order WENO scheme) and Einkemmer and Moriggl (2023) (extending their semi-Lagrangian discontinuous Galerkin scheme to 3D-3V and up to 1536 GPGPUs). With the advent of performance-portable acceleration libraries (such as Kokkos, Trott et al. 2022; RAJA, Beckingsale et al. 2019; alpaka, Matthes et al. 2017) and increasing maturity of standardized implementations (OpenACC and OpenMP), portable code updates have become more commonplace. Examples of Kokkos include a C++ mini-app implementing parts of the Gysela code (Asahi et al. 2020), the BSL6D 3D-3V semi-Lagrangian solver of the Vlasov-Poisson system, and the 2D-2V multi-GPU Vlasov-Poisson solver VCK-GPU (Ho and Vogman 2024). Other approaches of hardware-accelerated Vlasov codes include an OpenACC implementation of FARSIGHT, an 1D-1V Barycentric Lagrange treecode with adaptive mesh refinement (Sandberg et al. 2025) and a python-based discontinuous Galerkin approach solving semi-discrete equations in a tensor-product form through existing CUDA-Python-libraries. A further option is to write performance-portable code which is converted to CUDA or HIP on-demand, as seen in the initial GPU port of Vlasiator (Battarbee et al. 2025).
## 4.11 Achievements in Vlasov-based modelling

The progress of available scientific computing capabilities towards and beyond petascale in the last decade has driven the interest in and applicability of Vlasov-based methods to multidimensional space and astrophysical plasma problems. Tables 2 and 3 compile existing research work using direct solutions of the Vlasov

<table><caption>**Table 2** Space and astrophysical applications of Vlasov-based plasma simulation methods</caption><thead><tr><th>Application</th><th>Model characteristics</th><th>References</th></tr></thead><tbody><tr><td rowspan="5">Kelvin-Helmholtz instability</td><td>SL 2D-2V full Vlasov</td><td>Umeda et al. (2010a)</td></tr><tr><td>SL 2D-2V full Vlasov</td><td>Umeda et al. (2014, 2016); Umeda (2020)</td></tr><tr><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Settino et al. (2020)</td></tr><tr><td>FD 1D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Malara et al. (2018); Guzzi et al. (2021); Malara et al. (2022)</td></tr><tr><td>FV 2D-2V full Vlasov-Poisson</td><td>Vogman et al. (2020); Vogman and Hammer (2021)</td></tr><tr><td>Rayleigh-Taylor instability</td><td>SL 2D-2V full Vlasov</td><td>Umeda and Wada (2016, 2017)</td></tr><tr><td rowspan="2">Collisionality in the solar wind</td><td>FD 1D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Pezzi et al. (2016)</td></tr><tr><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Pezzi (2017); Pezzi et al. (2019c, 2017a, 2017b)</td></tr><tr><td>Magnetic holes</td><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Arrò et al. (2024)</td></tr><tr><td rowspan="2">Kinetic Alfvén Waves</td><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Pucci et al. (2016)</td></tr><tr><td>FD 1D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Maiorano et al. (2020)</td></tr><tr><td rowspan="2">Collisionless shocks</td><td>DG 1D-2V full-Vlasov</td><td>Juno et al. (2021)</td></tr><tr><td>DG 1D-3V full-Vlasov</td><td>Juno et al. (2023)</td></tr><tr><td>VLF radio emissions and whistler waves in the Earth's radiation belts</td><td>L 1D-3V hybrid-Vlasov (e<sup>-</sup>)</td><td>Nunn et al. (1997); Nunn (2005); Nunn et al. (2009); Omura and Nunn (2011); Nunn and Omura (2012); Harid et al. (2014); Nunn et al. (2021)</td></tr><tr><td rowspan="2">Solar wind interaction with unmagnetised or weakly magnetised bodies</td><td>SL 1D-3V hybrid-Vlasov (e<sup>-</sup>)</td><td>Gibby et al. (2008)</td></tr><tr><td>SL 2D-3V full Vlasov</td><td>Umeda et al. (2011); Umeda (2012); Umeda et al. (2013); Umeda and Ito (2014); Umeda and Fukazawa (2015)</td></tr><tr><td>Charge and potential distribution around a spacecraft</td><td>3D Vlasov-Poisson (iterative relaxation algorithm) and Vlasov-Laplace (Lagrangian)</td><td>Chane-Yook et al. (2006); Vauchelet et al. (2010)</td></tr><tr><td rowspan="4">Weibel instabilities</td><td>SL 1D-2V hybrid-Vlasov (e<sup>-</sup>)</td><td>Inglebert et al. (2011)</td></tr><tr><td>SL 2D-2V hybrid-Vlasov (e<sup>-</sup>)</td><td>Ghizzo et al. (2017, 2023, 2024)</td></tr><tr><td>SL 2D-2V and 2D-3V hybrid-Vlasov (e<sup>-</sup>)</td><td>Sarrat et al. (2017); Ghizzo and Del Sarto (2020)</td></tr><tr><td>SL 1D-2V full Vlasov</td><td>Ghizzo and Del Sarto (2021); Del Sarto et al. (2024)</td></tr></tbody></table>

<footnote>FD finite difference. FV finite volume. L fully lagrangian. SL semi-lagrangian. DG discontinuous galerkin. FC flux conservative. e⁻, H⁺, He⁺⁺: kinetic species (electrons, protons, helium ions) in a hybrid setup. Note that studies using Vlasiator are reviewed separately in Sect. 5.9</footnote>
<table><caption>**Table 3** Space and astrophysical applications of Vlasov-based plasma simulation methods: reconnection and turbulence</caption><thead><tr><th>Application</th><th>Model characteristics</th><th>References</th></tr></thead><tbody><tr><td rowspan="4">Magnetic reconnection</td><td>SL 2D-3V full Vlasov</td><td>Umeda et al. (2009, 2010a, b); Zenitani and Umeda (2014)</td></tr><tr><td>FC 2D-3V full Vlasov</td><td>Schmitz and Grauer (2006a, 2006b, 2008)</td></tr><tr><td>FD 2D-3V full Vlasov</td><td>Pezzi et al. (2019b)</td></tr><tr><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Califano et al. (2001); Finelli et al. (2021)</td></tr><tr><td>Turbulence-induced magnetic reconnection</td><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Cerri and Califano (2017); Franci et al. (2017); Fadanelli et al. (2021); Manzini et al. (2023); Foldes et al. (2024)</td></tr><tr><td rowspan="2">Electron-only Reconnection</td><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Califano et al. (2020); Arró et al. (2020)</td></tr><tr><td>FD 3D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Granier et al. (2024)</td></tr><tr><td rowspan="4">Solar wind turbulence</td><td>FD 3D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Servidio et al. (2015); Cerri et al. (2017b, 2018, 2019); Sisti et al. (2021)</td></tr><tr><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>)</td><td>Valentini et al. (2010, 2011); Servidio et al. (2012, 2014); Valentini et al. (2014); Vásconez et al. (2014, 2015); Valentini et al. (2016); Cerri et al. (2016); Leonardis et al. (2016); Cerri et al. (2017a); Valentini et al. (2017); Perrone et al. (2018); Pezzi et al. (2018, 2019d); Hu et al. (2020); Finelli et al. (2021); Pezzi et al. (2021); Sisti et al. (2021); Vásconez et al. (2021)</td></tr><tr><td>FD 2D-3V hybrid-Vlasov (H<sup>+</sup>, He<sup>++</sup>)</td><td>Perrone et al. (2013, 2014, 2014); Sorriso-Valvo et al. (2018)</td></tr></tbody></table>

FD finite difference. FV finite volume. L fully lagrangian. SL semi-lagrangian. DG discontinuous galerkin. FC flux conservative. e⁻, H⁺, He⁺⁺: kinetic species (electrons, protons, helium ions) in a hybrid setup. Note that studies using Vlasiator are reviewed separately in Sect. 5.9

equation in plasma physics. They focus on works with direct application to space and astrophysics, meaning that purely theoretical studies, as well as research from adjacent fields (e.g. nuclear fusion and laser-plasma interaction) have been omitted from this list. Studies based on the hybrid-Vlasov code Vlasiator are placed in Sect. 5.9 so they are not included in this list. The Vlasov equation has been used in space plasma physics and plasma astrophysics to model a variety of processes and plasma environments including turbulence and magnetic reconnection (see Table 3), plasma instabilities, collisionless shocks, magnetic holes, collisions in the solar wind, the interaction of the solar wind with the Earth and other bodies, radio emissions and whistler chorus waves in near-Earth space, as well as the charge distribution around spacecraft (see Table 2).
Among the various applications, turbulence has emerged as one of the most extensively studied processes using Vlasov-based approaches. Vlasov simulations of turbulent plasmas have been designed to mimic the conditions of the solar wind and the Earth's magnetosheath, enabling direct comparisons with *in situ* observations, particularly from MMS. A comprehensive and recent review of achievements in the modelling of plasma turbulence using Vlasov codes is presented in Califano and Cerri (2023). Magnetic reconnection represents another fundamental process thoroughly investigated using both hybrid (Califano et al. 2001; Finelli et al. 2021) and fully kinetic Vlasov models (Umeda et al. 2009; Schmitz and Grauer 2006b; Pezzi et al. 2019b). The most recent developments in this area are reviewed in Shay et al. (2025). Here, we provide further details about two processes that have received significant attention in recent years in term of Vlasov modelling, namely electron-only magnetic reconnection and Kelvin-Helmholtz instability (see also Sect. 2.5).

### 4.11.1 Electron-only reconnection

Recent hybrid-Vlasov simulations of plasma turbulence (Valentini et al. 2007) have provided key insights into the phenomenon of electron-only reconnection, where magnetic reconnection proceeds via electron-scale dynamics without ion participation. Using 2D-3V hybrid-Vlasov simulations, Califano et al. (2020) demonstrated that when turbulence is driven at scales approaching the ion kinetic range, electron-scale current sheets dominate the dissipation process, giving rise to electron-only reconnection events that qualitatively resemble those observed in the Earth's magnetosheath by MMS (Phan et al. 2018). Arró et al. (2020) examined the statistical signatures of turbulence in the presence of electron-only reconnection, finding no significant deviations in intermittency or energy cascade properties when compared to turbulence mainly populated by standard ion-coupled reconnection. These results suggest that the turbulent cascade at sub-ion scales is insensitive to the specific mechanism associated with magnetic reconnection. Extending the investigation of electron-only reconnection in turbulence to fully three-dimensional simulations, Granier et al. (2024) performed 3D-3V simulations with compressible, quasi-isotropic injection near ion scales. This study confirms the occurrence of electron-only reconnection in high-beta regimes and highlights the critical role of the ion gyroradius relative to the injection scale in enabling an electron-MHD regime. Moreover, they reported significant ion heating mediated by turbulent processes, even in the absence of ion outflows typically associated with reconnection. These studies underscore the ability of Vlasov-based methods to capture subtle kinetic features of electron-only reconnection in turbulent, collisionless plasmas.

### 4.11.2 Kelvin-Helmholtz Instability

A growing body of simulations based on the Vlasov approach has investigated the dynamics of the Kelvin-Helmholtz instability (KHI) in collisionless plasmas, Fully kinetic Vlasov simulations have demonstrated the ability of full electromagnetic Vlasov solvers to resolve both the primary vortex formation and the emergence of
secondary instabilities Umeda et al. (2010a, 2014, 2016); Umeda (2020). These results highlighted the role of anisotropic pressure and heat flux in modifying the linear instability growth and the formation of nonlinear structures. These studies showed that a kinetic approach is essential for evaluating gyroviscous terms and higher moments, which are seen to deviate significantly from analytical approximations. Also, they indicate that results from kinetic simulations and fluid simulations including the finite-Larmor-radius (FLR) term disagree when the FLR effects are included only in the equation of motion of the fluid model, thus challenging the applicability of FLR-MHD closures. Vogman et al. (2020); Vogman and Hammer (2021) implemented high-order continuum Vlasov solvers with exact two-species equilibria to systematically analyse the effects of magnetisation, diamagnetic drifts, and charge separation in low-beta, inhomogeneous plasmas. Their simulations showed substantial departures from Hall-MHD predictions in both growth rates and transport properties. They also revealed that distribution functions from kinetic simulations deviate substantially from Maxwellian distributions, indicating the importance of kinetic effects, especially during the nonlinear phase of the instability.

The Kelvin–Helmholtz instability has also been explored with hybrid codes. Several studies demonstrated the importance of accurate equilibrium modelling in capturing the full dynamical range of KHI evolution. Settino et al. (2020) employed the exact stationary equilibrium developed by Malara et al. (2018) to initialise a sheared flow configuration, showing that this choice is critical for enabling a robust turbulent cascade and for the formation of strongly non-Maxwellian ion distributions near shear layers. In contrast, simplified initial conditions based on Maxwellian distribution functions tend to suppress turbulence. Earlier work by Cerri et al. (2013) proposed approximate equilibria incorporating FLR corrections, while an exact equilibrium was derived by Malara et al. (2018) and later extended by Guzzi et al. (2021) and Malara et al. (2022) to include inhomogeneities in density, temperature, and magnetic field orientation, providing more realistic and stable configurations suitable for exploring kinetic instabilities at shear boundaries. Overall, these studies underscore the benefits of constructing stationary equilibrium configurations when investigating shear-driven instabilities in inhomogeneous plasmas.

## 5 Vlasiator

This section considers the choices and approaches made for the Vlasiator code describing the near-Earth space with hybrid ion-kinetic physics. The evolution of the phase-space density of kinetic ions is solved with the Vlasov equation (3), while the evolution of electromagnetic fields described through Faraday's law (8), Gauss' law and the solenoid condition (6), and Ampère's law without the displacement current (14). Electrons are modelled as a massless charge-neutralising fluid. Closure is provided via the generalised Ohm's law (12) under the assumptions of high conductivity, slow temporal variations, and polytropic electrons, i.e. the Hall MHD Ohm's law (13) with the electron pressure term. The source code of Vlasiator is available at https://github.com/fmihpc/vlasiator (Pfau-Kempf et al. 2024) according to the Rules of the Road mapped out at http://helsinki.fi/vlasiator.
## 5.1 Background

Vlasiator has its roots in the discussions within the global MHD simulation community around 2005. It was becoming evident that while global MHD simulations are important, their capabilities, especially in the inner magnetosphere, are limited. The inner magnetosphere consists of spatially overlapping plasma populations of different temperatures (e.g., Baker 1995) and therefore the environment cannot be satisfactorily modelled with MHD to a degree allowing, for example, environmental predictions for societally critical spacecraft or as a context for upcoming missions, like the Van Allen Probes (e.g., Fox and Burch 2014). To this end, two strategies emerged, including either coupling a global MHD simulation with an inner magnetospheric simulation (e.g., Huang et al. 2006), or going beyond MHD with the then newly introduced hybrid-PIC approach (e.g., Omidi and Sibeck 2007). Coupling different codes carries a risk that the effects of the coupling scheme dominate over the improved physics. On the other hand, while hybrid-PIC simulations had produced important breakthroughs (e.g. Omidi et al. 2005), the velocity distributions computed through binning are noisy due to the limited number of launched particles, which could compromise physical conclusions. Further, due to the limited number of particles, the hybrid-PIC simulations could not provide sharp gradients, which would become a problem especially in the magnetotail, where the lobes surrounding the dense plasma sheet are almost empty. As accurately modelling magnetotail physics is critical in the global description of the magnetosphere, the idea about a hybrid-Vlasov simulation emerged.

The objective was simple, just to go beyond MHD by introducing protons as a kinetic population modelled by a distribution function and thus avoiding the numerical noise inherent to the PIC approach. Several challenges were identified. First, if one neglects electrons as a kinetic population, one will, for example, not be able to model electron-scale instabilities that can be important in magnetotail physics (e.g., Pritchett 2005). Second, a global hybrid-Vlasov approach is still an extreme computational challenge even with a coarse ion-scale resolution, since it must be carried out in six dimensions. Further, doubts existed about whether grid resolutions achievable with the computational resources available at that time would facilitate ion kinetic physics. However, with a new approach without historical heritage, one could utilise the latest HPC methods and new computational architectures, provided that the code would always be portable to the latest technology. The computational resources were still obeying the “Moore law”, and petascale systems had just become operational (Kogge 2009). With these prospects in mind, Vlasiator was proposed to the newly established European Research Council in 2007, which solicited new ideas with a high risk – high gain vision (Palmroth 2022).

Since 2018, the above notions are still valid, and supercomputer technology has continued to change and diversify. Vlasiator’s strategy was always to target to upcoming infrastructure, so that when each architectural revolution took place, new hardware and features could be immediately utilised. In recent years, this has led to the need of adapting the code to harness heterogeneous supercomputers containing a large amount of graphics processing units.
## 5.2 Grid discretisation

Vlasiator discretises the 6D space of Vlasov’s equation separately into the 3D position space and the 3V velocity space.

### 5.2.1 Spatial grid discretisation

The position space is discretised on a Cartesian grid of cubic cells. Each cell holds the volume-averaged values of variables that are either being propagated or reconstructed (e.g. the electric and magnetic fields and the ion velocity distribution moments) as well as housekeeping variables. In addition, the three-dimensional uniform Cartesian velocity space grid (see Sect. 5.2.2) is stored in each spatial cell. For position space Vlasiator uses the Distributed Cartesian Cell-Refinable Grid library (DCCRG; Honkonen et al. 2013). The library can distribute the grid over a large supercomputer using the domain decomposition approach (see Sect. 5.6 for details on the parallelisation strategies).

The position space domain can be decomposed with a uniform resolution throughout, as has been done extensively in particular with 1D or 2D domains with the transverse directions assuming periodicity. Recent developments (Ganse et al. 2023) have enabled 3D simulations to be run with static mesh refinement in predefined volumes of the simulation domain (e.g. Palmroth et al. 2023), focusing resolution and thus computations on regions of scientific interest and importance in overall dynamics. Further developments yet allow for adaptive mesh refinement (AMR) to track dynamically regions of interest following a set of fluid-based and kinetic criteria (i.e., local differences of macroscopic variables such as particle and magnetic flux density, and the ratio of current density to the magnetic field component perpendicular to current Kotipalo et al. 2024). Both static and adaptive mesh refinements are only applied for the grid containing variables utilised by the Vlasov equation solvers, and run-time dynamic load-balancing strategies are applied on this grid (see Sect. 5.6).

The position space is decomposed separately on a uniform grid spanning the whole domain for the purposes of field solving, using another dedicated library (https://github.com/fmihpc/fsgrid). This grid is decomposed uniformly across all computational tasks since the field solver has a uniform computational cost, unlike the Vlasov equation part of the solvers. Its cells are congruent with the highest-resolution level cells of DCCRG grids applying mesh refinement. Both the Vlasov solver and field solver grids are coupled at every time step, mapping plasma moments to the field solver grids prior to the update of the electric and magnetic fields, and mapping the fields back to the Vlasov solver grid thereafter (Papadakis et al. 2022).

### 5.2.2 Velocity grid discretisation

The velocity space grid in Vlasiator is purpose built for managing non-Maxwellian distribution functions which can vary in size and complexity throughout the simulation domain. A major performance gain in terms of memory and computation is achieved by storing and propagating the volume average of the distribution
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Illustration of the velocity space grid and block resolution as well as the phase-space density threshold $f_{min}$ for the case of a Maxwellian distribution computed for protons of density $n = 10^6 \text{ m}^{-3}$, temperature $T = 0.5$ MK and a drift velocity of 500 km s$^{-1}$. Similar plots can be reproduced with the script `tools/vdf_params.py` included in the Vlasiator source code (Pfau-Kempf et al. 2024)</figcaption></figure>

function $f$ in every cell at position $\vec{x}$ only in those velocity space cells where $f$ exceeds a given density threshold $f_{min}$. In order to accurately model propagation and acceleration, a buffer layer is maintained by modelling also cells that are adjacent in position or velocity space. The principle is illustrated in Fig. 4 for an artificial 1D case where a Maxwellian distribution is plotted along with the velocity space resolution and $f_{min}$ as well as in Fig. 5 for the case of a distribution function extracted from a magnetosphere simulation. This density threshold $f_{min}$ can be constant or scaled linearly with the ion number density. For each ion population, the maximal velocity space extents and the resolution are set by the user. This so-called sparse velocity space strategy allows the user to simulate the distribution function only in the regions where there is significant phase-space density instead of wasting computational resources covering the full extents of reachable velocity space. It is however important to set the value of $f_{min}$ carefully in order to conserve the moments of $f$ (density, pressure, etc.) to the desired accuracy and in order to include all expected features such as low-density high-energy populations. See Fig. 1 and related discussion by von Alfthan et al. (2014) for more details on the effect of sparsity on conservation of moments. Furthermore, a detailed discussion of the effects of the grid discretisation parameters on the simulation of a collisionless shock was reported by Pfau-Kempf et al. (2018).

A further point of adjustment is the granularity of activating or de-activating regions of velocity space. Vlasiator supports bookkeeping of phase-space cells in units (blocks) of $4 \times 4 \times 4$ cells (such as shown in Figs. 4 and 5), or with a recent addition (Battarbee et al. 2025), of $8 \times 8 \times 8$ cells. This unit defines the computational unit which can be utilized by the Vlasov solvers and thus affects parallelization strategies. A larger block size is also able to manage bookkeeping for a similar distribution function with less bookkeeping overhead, at the cost of potentially simulating larger buffer regions. This trade-off is exemplified in Fig. 6.
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Illustration of the sparse velocity space. The colour shows the phase-space density in arbitrary units. Left: full extent of velocity space including a population drifting at $V_x \approx -500$ km/s. Middle: Cut through the population. Right: Slice showing the cells with $f > f_{min}$ in full colour and retained neighbours greyed out. 1: block fully above $f_{min}$. 2: block partially above $f_{min}$. 3: retained neighbouring block. 4: disregarded block without neighbours above $f_{min}$. See Sect. 5.6 and Fig. 6 for details on the block structure. Figure from Pfau-Kempf (2016)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6** Illustration of the block-based velocity mesh. Panel a): Velocity domain consisting of blocks of 4x4x4 cells. A red block (hidden) is flagged as containing significant plasma density, causing neighbouring blocks (face neighbours in dark blue, edge and corner neighbours in fuchsia) to be retained in the sparse mesh structure. b): Exploded view of the velocity mesh, showing the red content block. c): Exploded view of the same physical velocity distribution, represented on a mesh consisting of 8x8x8 cell blocks instead of 4x4x4 cell blocks. d) Velocity domain consisting of blocks of 8x8x8 cells. Increasing block size reduces the amount of metadata required for representing the velocity grid, but on average increases the amount of velocity cells used to represent a given velocity distribution function</figcaption></figure>
## 5.3 Solvers and time-integration

The structure of the hybrid-Vlasov set of equations leads to the logical split into a solver for the Vlasov equation and a solver for the electric and magnetic field propagation.

### 5.3.1 Vlasov solver

In advancing the Vlasov equation, Vlasiator utilises Strang splitting (Umeda et al. 2009, 2011, and references therein), where updates of the particle distribution functions are performed, separately using a spatial translation operator $S_T$ for advection; $S_T = (\vec{v} \cdot \frac{\partial f}{\partial \vec{v}})$, and an acceleration operator $S_A = (\frac{g}{m}\vec{E} \cdot \frac{\partial f}{\partial \vec{v}})$ including rotation $(\frac{g}{m}(\vec{v} \times \vec{B}) \cdot \frac{\partial f}{\partial \vec{v}})$ for each phase-space volume average. The splitting is performed using a standard leapfrog scheme, where

$$ \tilde{f}^{N+1} = S_A(\frac{\Delta t}{2})S_T(\Delta t)S_A(\frac{\Delta t}{2})\tilde{f}^N. $$

Acceleration over timestep $\Delta t$ is thus calculated based on the field values computed at the mid-point of each acceleration step, i.e. at each actual time step as used for translation.

A global time step is defined, with time advancement calculated for distribution functions and fields in separate yet linked computations. Earlier versions of Vlasiator have used finite volume (FV) Vlasov solvers. In the earliest versions of the code, a FV method based on the solver proposed by Kurganov and Tadmor (2000) was used (Palmroth et al. 2013). A Riemann-type FV solver (LeVeque 1997; Langseth and LeVeque 2000) was used in subsequent works (Kempf et al. 2013; Pokhotelov et al. 2013; Sandroos et al. 2013; von Alfthan et al. 2014; Sandroos et al. 2015; Kempf et al. 2015). For these solvers, the classical Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928) for maximal allowable time steps when calculating fluxes from one phase-space cell to another is

$$ \Delta t < \min \left( \frac{\Delta x_i}{\max(|v_i|)}, \frac{\Delta v_i}{\max(|a_i|)} \right), $$

where $a_i$ is the acceleration $\partial v_i / \partial t$ and $i$ is indexed over three dimensions. In previous versions, the CFL condition was found to be very limiting. Instead, current versions of Vlasiator utilize a semi-Lagrangian scheme (SLICE-3D, Zerroukat and Allen 2012; lib-slice3d 2016), in which mass conservation is ensured by a conservative remapping from a Eulerian to a Lagrangian grid. Note however that the sparse velocity space strategy implemented in Vlasiator (see Sect. 5.2.2) breaks the mass conservation (see Pfau-Kempf et al. 2018, for a discussion of the effect of the phase space density threshold on mass conservation). The specificity of the SLICE-3D scheme is to split the full 3D remapping into successive 1D remappings, which reduces the computational cost of the spatial translation and facilitates its parallel implementation.
The velocity space update due to acceleration $S_A(\frac{\Delta t}{2})$ will generally be described by an offset 3D rotation matrix (due to gyration around $\vec{B}$). As every offset rotation matrix can be decomposed into three shear matrices $S = S_x S_y S_z$, each performing an axis-parallel shear into one spatial dimension (Chen and Kaufman 2000), a numerically efficient semi-Lagrangian acceleration update using the SLICE-3D approach is possible. Before each shear transformation, the velocity space is rearranged into a single-cell column format parallel to the shear direction in memory, each of which requires only a one-dimensional remapping with a high reconstruction order (in Vlasiator, 5th order reconstruction is typically employed for this step). These column updates are optimised to make full use of vector instructions. This update method comes with a maximum rotation angle limit due to the shear decomposition of about 22°, which imposes a further time step limit. For larger rotational angles per time step (caused by stronger magnetic fields), the acceleration can be subcycled.

The position space update $S_T(\Delta t)$ will generally be described by a translation matrix with no rotation, and the same SLICE-3D approach lends itself to it in a similar vein as for velocity space. The main difference is the typical use of 3rd order reconstruction in order to keep the stencil width at two. The use of a semi-Lagrangian scheme allow the implementation of a time step limit

$$ \Delta t < \min \left( \frac{\Delta x_i}{\max(|v_i|)} \right) $$

based on spatial translation only. This condition constrains the spatial translation of any volume averages to a maximum value of $\Delta x_i$ in direction $i$, accounting for only those velocities within phase-space which have populated active cells (see the sparse grid implementation, Sect. 5.2). This is employed not due to stability requirements, but rather to decrease communication bandwidth by ensuring that a single ghost cell in each spatial direction is sufficient.

### 5.3.2 Field solver

The field solver in Vlasiator is based on the upwind constrained transport algorithm by Londrillo and Del Zanna (2004) and uses divergence-free reconstruction of the magnetic fields (Balsara 2009). It utilizes a second-order Runge-Kutta algorithm including the interpolation method demonstrated by Valentini et al. (2007) to obtain the intermediate moments of $f$ needed to update the electric and magnetic fields (von Alfthan et al. 2014).

The algorithm is subject to a CFL condition such that the fastest-propagating wave mode cannot travel more than half a spatial cell per time step. As the field solver was extended to include the Hall term in Ohm's law, the CFL limit severely impacts the time step of the whole propagation in regions of high magnetic field strength or low plasma density. If the imbalance between the time step limits from the Vlasov solver and from the field solver is too strong, the computation of the electric and magnetic fields is subcycled such as to retain an acceptable global time step length (Pfauf-Kempf 2016).
The field solver equations can also be set to include an approximation of the electron pressure gradient term. Although Vlasiator does not simulate electrons, it is possible to estimate the electron pressure at any point within the simulation if some reference values and a suitable polytropic index are provided. For example, by providing solar wind inflow electron densities and temperature, an inflow pressure can be evaluated. If electrons are considered, for example, to behave in an adiabatic fashion (polytropic index $k = 5/3$), and quasi-neutrality is enforced, the local electron density $n_e$ can be evaluated as being equal to the charge density of ions $\rho_q$ divided by the elementary charge $e$. Then the local (isotropic) electron pressure can be solved from the electron equation of state, through evaluation of the equation

$$P_e n_e^{-k} = \text{constant.} \qquad (23)$$

Solving the electron pressure for each simulation cell, local gradients can be evaluated and fed into the Ohm's law as the electron pressure gradient term $-\frac{\nabla \cdot P_e}{n_e e}$.

### 5.3.3 Time stepping

The leapfrog scheme of the propagation is initialised by half a time step of acceleration. If the time step needs to change during the simulation due to the dynamics of the system, $f$ is accelerated backwards by half an old time step and forwards again by half a new time step and the algorithm resumes with the new global time step. The complete sequence of the time propagation in Vlasiator is depicted in Fig. 7, including a synthetic version of the equations used in the different parts.

<figure><img src="image_9.png" /><figcaption>Fig. 7 Time stepping in Vlasiator. The translation and acceleration of $f$ are leapfrogged following the Strang-splitting method. The algorithm is initialised by half a time step of acceleration (step 0. in red). Then 1. $f$ is translated forward by one step $\Delta t$ (possibly subcycled, see text). 2. E, B are stepped forward by $\Delta t$ (possibly subcycled, see text). 3. $f$ is accelerated forward by $\Delta t$. The sequence is repeated (4.–6.)</figcaption></figure>
## 5.4 Boundary conditions

While Vlasiator can simulate completely periodic setups, boundary conditions are needed for more complex setups such as magnetospheric simulations.

### 5.4.1 Domain outer boundaries

The reference frame used in Vlasiator is defined as follows: the origin is located at the centre of the Earth, the *x* axis points towards the Sun, the *z* axis points northward perpendicular to both the *x* axis and the ecliptic plane, and the *y* axis completes the right-handed set. This coordinate system is equivalent to the Geocentric Solar Ecliptic (GSE) frame which is commonly used when studying the near-Earth space.

The solar wind enters the simulation domain at the +*x* boundary, while copy conditions (i.e., homogeneous Neumann boundary conditions obtained by copying the velocity distribution functions and magnetic fields from the boundary cells to their neighbouring ghost cells) or static conditions (i.e., initial values are kept constant throughout one boundary plane) can be used for the -*x* boundary and for the boundaries perpendicular to the solar wind flow. In 2D-3V runs, as well as in the limited 3D run used by Pfau-Kempf et al. (2020), periodic conditions are applied in the out-of-plane direction (i.e., +*z* and -*z* for the ecliptic runs and +*y* and -*y* for the polar runs).

The outer boundary conditions are allowed to be time-dependent (Zhou et al. 2022). These conditions are specified using a time series of solar wind parameters: (*t*, *n*<sub>s</sub>, *T*<sub>s</sub>, *U*<sub>s</sub>, *B*), where *t* is time, *B* is the solar wind magnetic field, and *n*<sub>s</sub>, *T*<sub>s</sub>, and *v*<sub>s</sub> are the ion number density, scalar temperature, and bulk velocity for ion species *s*. The time-varying inflow conditions are implemented as piece-wise linear functions, configured by defining solar wind and IMF parameter values at specific times in the simulation. At preset time intervals, Vlasiator updates the phase space distributions in the boundary cells to match the updated conditions. To achieve this, a linear interpolation between successive data points in the solar wind input parameter list is performed if the simulation time falls between two time stamps in the input data. If the simulation time exceeds the range of the input data, the values from the last time stamp are used. The input data is structured such that each file corresponds to a single ion species, facilitating straightforward support for multi-species simulations.

### 5.4.2 Inner boundary

The simulation also requires an inner boundary around the Earth, in order to screen the origin of the terrestrial dipole. In the inner magnetosphere, the magnetic field strength increases dramatically with decreasing distance to the Earth's centre, resulting in very small time steps, which would significantly slow down the whole simulation. Also, close to the Earth, the ionospheric plasma can no longer be described as a collisionless and fully ionised medium, and another treatment would be required in order to properly simulate this region. The inner boundary is therefore typically located at 30, 000 km (about 4.7 R<sub>E</sub>) from the Earth's centre, yielding a timestep length of the order of 0.01 s which is compatible with available
computational resources and the overall computational burden of the systems to simulate (see Sect. 4.10).

Until version 5.2 of Vlasiator, the inner boundary in 2D and 3D setups (thus including the first 3D simulation studied by Palmroth et al. 2023) was constraining tangential and radial components of the magnetic field (akin to emulating a perfect conductor behaviour), but presented issues occasionally leading to numerical instability of the model. This was alleviated by implementing a scheme applying Neumann conditions to the perturbed magnetic field at the inner boundary radius, or optionally keeping the field static in the boundary layer, which has proven to be more stable in 2D setups. In all of these versions, the distribution functions in the boundary cells retain their initial Maxwellian distributions throughout the run (see Sect. 5.5 for details on the initial conditions).

The above inner-boundary conditions nevertheless did not provide a representation of the feedback from the ionosphere to the simulation domain, thus limiting the capability to study the nightside transition region and inner-magnetospheric processes. To remediate this issue, the latest step in improving the fidelity of the inner boundary behaviour in 3D modelling with Vlasiator from version 5.2 onwards has been the implementation of a coupled ionosphere model (Ganse et al. 2025), which builds on the heritage of global magnetohydrodynamic models (e.g., Goodman 1995; Janhunen et al. 2012). The ionosphere module maps plasma parameters and field-aligned currents along the geomagnetic field lines from near the inner boundary down to an infinitely thin ionospheric layers at an altitude of 100 km above the ground. These values are used to obtain altitude-resolved electron production rates, from which the ionospheric Hall and Pedersen conductivities are calculated. The conductivities are then height-integrated to obtain the corresponding conductances and the current continuity equation is solved on a spherical grid at 100 km altitude to obtain the electric potential of the ionosphere. This potential is mapped back along the geomagnetic field lines to the inner boundary of the magnetospheric Vlasov domain, where its gradient is used to determine an $\vec{E} \times \vec{B}$ drift that is applied to the boundary cells' ion distribution function. This allows a more realistic modelling of global magnetospheric circulation compared to previous versions of Vlasiator or models applying constant ionospheric conductivities. The ionosphere coupling is detailed in the work by Ganse et al. (2025).

## 5.5 Initial conditions

In addition to defining boundary conditions, the phase-space cells within the simulation box must be initialised to some reasonable values, after which the magnetic and gas pressures and flow conditions cause the state of the simulation to change, and to converge towards a valid description of the magnetospheric system as the box is flushing. The usual method employed in Vlasiator is to initialise the velocity space in each cell within the simulation (excluding the region inside the inner boundary) to match values sampled from a Maxwellian distribution in agreement with the inflow boundary solar wind density, temperature, and bulk flow
direction. The inner boundary is initialised with a constant proton temperature and number density with no bulk flow.

The Earth's magnetic field closely resembles that of a magnetic dipole, and within the scope of Vlasiator the dipole has been approximated as being aligned with the z-axis. For ecliptic x-y plane simulations the dipole field can be used as-is, but for polar x-z plane simulations, a 2D line dipole (which scales as $r^{-2}$ rather than $r^{-3}$) must be used instead in order to prevent the occurrence of unphysical currents due to out-of-plane field curvature. When using this approach, one must calculate the dipole strength that represents reality in some chosen manner, and this is achieved by choosing a value that in turn reproduces the magnetopause at the realistic ~10 R_E standoff distance (a similar treatment as found in, e.g., Daldorff et al. 2014). This however prevents the direct mapping of positions in the inner magnetosphere to ionospheric altitudes in terms of geomagnetic latitude, thus limiting the capability to study magnetosphere-ionosphere coupling mechanisms in this setup, as discussed in Grandin et al. (2019).

As the dipole magnetic field is not included in the inflow boundary, it has to be enforced that the boundary-perpendicular magnetic field component is zero in order to respect the solenoid condition. For ecliptic runs, the dipole field component is aligned with z and thus there is no component perpendicular to the inflow boundary. For polar runs, the dipole field component perpendicular to the inflow boundary must be removed to prevent magnetic field divergence. This is achieved by placing a (curl-free) mirror dipole identical to the Earth's dipole model at the position $(2 \cdot (X_1 - \Delta x), 0, 0)$, with $X_1$ the x-coordinate of the inflow outer boundary, that is at twice the distance from the origin to the edge of the final non-boundary simulation cell. For each simulation cell, the static background magnetic flux through each face is thus assigned as a combination of flux calculated from the chosen dipole field model, a mirror dipole if present, and the IMF vector. This background field, which is curl-free and divergence-free, is left static, and instead any calculations involving magnetic fields operate on a separate field which acts as a perturbation from this initial field.

For 3D simulations, instead of applying a mirror dipole, an initial magnetic field setup was developed that smoothly transitions from a pure dipole field in the inner domain to the pure IMF closer to the inflow plane. To ensure the initial setup is divergence-free, the magnetic field is defined through the magnetic vector potential, which for a regular (ideal) dipole would be

$$ \vec{B}'(\vec{r}) = \nabla \times \vec{A}'(\vec{r}) = \frac{\mu_0}{4\pi} \left[ \frac{3\vec{r}(\vec{m} \cdot \vec{r})}{r^5} - \frac{\vec{m}}{r^3} \right], \quad (24) $$

where $\vec{r}$ is a distance vector pointing from the source of the dipole to the point of measurement and $\vec{m}$ is the Earth's magnetic dipole moment vector. To scale this dipole field to zero near the inflow boundary, the vector potential

$$ \vec{A}(\vec{r}) = \frac{\mu_0 \vec{m} \times \vec{r}}{4\pi r^3} \quad (25) $$

is scaled gradually to zero (over a user-defined distance) and the magnetic field $\vec{B}(\vec{r})$
is then solved as the curl of $\vec{A}(\vec{r})$. To retain one curl-free background field with non-curl-free components stored in another separate (perturbation) field, the background field is constructed by summing equation (24) with the constant IMF value, and the difference of the inflow-to-zero-scaled field $\vec{B}(\vec{r})$ and the ideal dipole field $\vec{B}'(\vec{r})$ is stored in the perturbed field.

## 5.6 Parallelisation strategies

Given the curse of dimensionality in Vlasov simulations (see Sect. 4), the amounts of memory and computational steps required for global magnetospheric hybrid-Vlasov simulations are extreme. Therefore, the use of HPC resources and parallelisation techniques is essential. Vlasiator uses three levels of parallelisation, of which the first employed is the decomposition of the spatial simulation domain into subdomains handling individual tasks using the Message-Passing Interface (MPI, MPI Forum 2015). The use of the DCCRG grid library (Honkonen et al. 2013) provides most of the interface for MPI communication and management of computational domain interfaces.

Thanks to the sparse velocity space representation, large savings in memory usage and computational demand can be achieved. However, the sparse velocity space induces a further problem, because the computational effort to solve the Vlasov equation is no longer constant for every spatial simulation cell, but varies in direct relation to the complexity of the velocity space at every given point. Due to the large variety of physical processes present in the magnetospheric domain, this leads to large load imbalances throughout the simulation box, making a simple Cartesian subdivision of space over computational tasks infeasible and necessitating a dynamic

<figure><img src="image_10.png" /><figcaption>**Fig. 8** Spatial plots of dynamic load balancing in Vlasiator in a global magnetospheric simulations. **Left** 2D simulation in the noon–midnight meridional plane. The simulation frame is identical to that of Fig. 1. Each computational rank 0–4799 is mapped to a colour, and the corresponding rectangular domain is coloured accordingly. Domain decomposition is performed by run-time updated recursive coordinate bisection using the Zoltan library (Devine et al. 2002; Boman et al. 2012). **Right:** equatorial-plane ($z = 0$) cut through a 3D simulation, decomposed using the recursive inertial bisection strategy</figcaption></figure>
rebalancing of the work distribution. To this end, Vlasiator relies on the Zoltan library (Devine et al. 2002; Boman et al. 2012), which creates an optimised space decomposition from continuously updated run-time metrics, providing a number of different algorithms to do so (usual production runs are using recursive inertial bisection). The left panel of Fig. 8 shows an example of the resulting spatial decomposition in a 2D global magnetospheric simulation run in the Earth’s noon-midnight meridional plane, in which the incoming solar wind plasma with low-complexity Maxwellian velocity distributions on the right-hand side is processed by visibly larger computational domain sizes than the more complex velocity space structures in the magnetosphere. The same effect can be seen in the 3D equatorial run in the right panel. With increasing MPI task counts, time spent in communication may grow prohibitively large, in which case alternative load distribution strategies may be considered. Vlasiator has been extended to support running the field or Vlasov solvers on dedicated sets of tasks, which may or may not overlap. This could allow simpler network communication patterns in particular inside the field solver, but still has the limitation of one solver requiring data from other solvers in order to be able to continue. Thus, to have an efficient solver decomposition, each task should be assigned independent computations to perform whilst waiting for necessary data from the coupled solver.

The second level of parallelisation is carried out within the computational domain of each MPI task. The domain containing the MPI tasks is typically handled by a full supercomputer node (or a fraction of it) with multiple CPU cores and threads, which includes a local parallelisation level based on OpenMP (OpenMP Architecture Review Board 2011). All computationally intensive solver steps have been designed to be run thread-parallel over multiple spatial cells, or in the case of the SLICE-3D position space update (see Sect. 5.3.1) multiple parallel evaluations over the velocity space to make optimum use of the available shared memory parallel computing architecture within one node. Parallelisation efficiency over these first two levels can be challenging, but gathering metrics or counters and interfacing with profiling tools and libraries (see, e.g., Coti et al. 2024) can aid in identifying bottlenecks.

As a third level of parallelisation, all data structures involved in computationally expensive solver steps have been designed to benefit from vector processing of modern CPUs. Specifically, the velocity space representation in Vlasiator is based on 4 × 4 × 4 cell blocks, which are always processed as a whole. This allows multiple velocity cells to be solved at the same time, using single instruction multiple data techniques (Fog 2016).

### 5.6.1 Parallelisation on accelerated hardware

Beyond general-purpose computing partitions, supercomputing clusters are increasingly investing in accelerated partitions where a CPU is coupled with one or more accelerator chips such as GPGPUs or GPUs for short, as introduced in Sect. 4.10.1. An early version of Vlasiator utilizing a finite volume solver was parallelised using the CUDA standard and run on small numbers of GPUs (Sandroos et al. 2013). In more recent years, Vlasiator with its current semi-Lagrangian solvers has been ported to support accelerated hardware through a performance-portable approach wherein
solver code is converted to either the CUDA or HIP languages, depending on the language supported by the system. One major challenge in porting Vlasiator to accelerated hardware is the way in which phase-space sparsity is handled (von Alfthan et al. 2014; Battarbee et al. 2025), in that the local grid in each spatial cell changes several times per time step. Careful use of different types of memory allocations, re-designing algorithms to be more suitable for a parallel implementation, and fusing kernels into larger launches has already resulted in a code architecture which with further optimisations and improved memory access patterns should allow exascale 3D-3V simulations of the terrestrial magnetosphere in the near future.

### 5.6.2 Parallel input / output

A further complication of parallel Vlasov simulations are the associated input/output requirements. Not only does it require a parallel input/output system that scales to the required number of nodes, but the sparse velocity space structure requires an appropriate file format able to represent the sparsity, without relying on fixed data offsets. For Vlasiator's specific use case, the VLSV library and file format have been developed (https://github.com/fmihpc/vlsv). Using parallel MPI-IO (MPI, MPI Forum 2015), it allows high-performance input/output even for simulation restart files which, given the large system size of Earth's magnetosphere, tend to get up to multiple terabytes in size. A plug-in for the scientific visualisation suite VisIt (Childs et al. 2012) is available, as is a python library that allows for quantitative analysis of the output files (Battarbee et al. 2024).

## 5.7 Verification and test problems

As standard verification tests for a hybrid-Vlasov system do not exist, the first verification effort of Vlasiator was presented by Kempf et al. (2013). A simulation of low-$\beta$ plasma waves (where $\beta$ is the ratio of thermal and magnetic pressures) in a one-dimensional case with various angles of propagation with respect to the magnetic field was used to generate dispersion curves and surfaces. These were then compared to analytical solutions from the linearised plasma wave equations given by the Waves in Homogeneous, Anisotropic Multicomponent Plasmas (WHAMP) code (Rönnmark 1982). Excellent agreement between the results obtained from the two approaches was found in the case of parallel, perpendicular and oblique propagation, the only noticeable difference taking place for high frequencies and wave numbers, likely as a result of too coarse a representation of the Hall term in the Vlasiator simulations at that time.

In the work presented by von Alfthan et al. (2014), the study of the ion/ion right-hand resonant beam instability is another effort to verify the hybrid-Vlasov model implemented in Vlasiator against the analytic solution of the dispersion equation for that instability. The obtained instability growth rates were found to behave as predicted by theory in the cool beam regime, although with slightly lower values which can be explained by the finite size of the simulation box used. That paper also compared results from the hybrid-Vlasov approach with those obtained with hybrid-
PIC codes and showed that the distribution functions are comparable albeit smoother and better-resolved with the former approach.

The increasing computational performance of Vlasiator has allowed significant improvements in spatial resolution of 2D-3V runs. It was still 850 km early on von Alfthan et al. (2014); Kempf et al. (2015) but subsequent runs were performed at 300 km and even 227 km resolution (e.g. Palmroth et al. 2015; Hoilijoki et al. 2017). Nevertheless, even at these finer resolutions the typical kinetic scales are still not properly resolved in magnetospheric simulations. This can lead to the a priori concern that under-resolved hybrid-Vlasov simulations would not fare better than their considerably cheaper MHD forerunners and would similarly lack any kinetic plasma phenomena. A systematic study of the effects of the discretisation parameters of Vlasiator on the modelling of a collisionless shock alleviates this concern. Using a one-dimensional shock setup with conditions comparable to the terrestrial bow shock, Pfau-Kempf et al. (2018) show that even at spatial resolutions of 1000 km the results clearly depart from fluid theory and are consistent with a kinetic description. Of course an increased resolution of 200 km leads to a dramatic improvement in the physical detail accessible to the model, even though not yet fully resolving ion kinetic scales. That study also highlights the importance of choosing the velocity resolution and the phase-space density threshold $f_{min}$ carefully as they affect the conservation properties of the model and as a consequence the physical processes it can describe (see Sect. 5.2.2).

An evaluation of how the simulation spatial resolution affects kinetic phenomena such as the mirror and proton cyclotron instabilities, phenomena commonly observed in the magnetosheath, was performed in Dubart et al. (2020). The growth rates of the mirror and proton cyclotron instability were compared against dispersion solver outputs, and it was shown that a resolution of 300 km was sufficient for capturing many aspects of magnetosheath dynamics, but under-resolving these waves caused them to decrease in effect with reduced resolution. At 900 km resolution, EMIC waves were not found in the magnetosheath, leading to reduced pitch-angle diffusion of particle populations and a larger-than-physical temperature anisotropy. As performing simulations with 3D ordinary space and 3D velocity space is computationally very demanding, currently the highest achievable spatial resolution is 1000 km in 3D-3V runs using AMR (see Sect. 5.2.1). Increasing resolution to sufficiently model these instabilities is currently not feasible, and thus a sub-grid model was developed (see Sect. 4.9) to mimic these physical phenomena (Dubart et al. 2022, 2023).

## 5.8 eVlasiator

A specialized version of Vlasiator has been developed for the purpose of investigating possible electron distribution functions generated under the influence of global, ion-scale fields. This eVlasiator approach (Battarbee et al. 2021, and Sect. 3.2.3) resembles electron-kinetic hybrid-Vlasov models in that electron motion is solved as a 2D-3V or 3D-3V advection equation according to the Vlasov equation and ions are considered a stationary background distribution. However, the closure provided through the field solvers is simplified substantially, as magnetic fields are
also considered static on the timeframes of simulated dynamics (typically on the order of a second) and the only electric field term included is that of electrostatic oscillation. This eVlasiator approach was, however, found to provide remarkable agreement with spacecraft observations (Alho et al. 2022). Further, eVlasiator has been used to provide non-Maxwellian electron precipitation estimates from a global ion-kinetic simulation as an input to the WACCM climate model, to improve modelling of upper atmospheric chemistry (Häkkilä et al. 2024).

## 5.9 Physics results

This section presents science results and validation by comparison with experimental ground-based or spacecraft data after 2018. Prior results were presented in Palmroth et al. (2018a). Vlasiator has proved particularly fruitful in modelling the complex dynamics arising from scale coupling – how small scales affect large scale dynamics and vice versa.

### 5.9.1 Terrestrial foreshock

Recent foreshock (see Fig. 1) studies include both science results and the description of new techniques of data analysis that can be applied to *in situ* satellite observations. An example of the latter is found in Dorfman et al. (2023), which uses Vlasiator data to highlight problems that can arise when applying single-spacecraft wave analysis techniques that assume infinite wave fronts to realistic conditions, especially close to the edge of the foreshock. Notably, the global view provided by Vlasiator allows investigating the complex interplay between shock and foreshock dynamics. Johlander et al. (2022) study bow shock reformation due to non-linear foreshock structures, and develop a new multi-spacecraft technique to identify such events. Battarbee et al. (2020b) show that protons are accelerated at a large portion of the bow shock, and that this is influenced by the magnetic field variation caused by the foreshock waves. Foreshock wave properties are themselves strongly affected by solar wind conditions, and particularly the interplanetary magnetic field which modifies their wavelength and their transverse extent (Turc et al. 2019). This in turn results in irregularities of the shock front at shorter scales, which affect ion reflection. Non-linear foreshock structures, cavitons and spontaneous hot flow anomalies (SHFAs), are formed throughout the foreshock and can cause large-scale erosion of the bow shock (Blanco-Cano et al. 2018). The formation and evolution of these structures is investigated in Tarvus et al. (2021), showing that a large fraction of cavitons turns into SHFAs as they approach the bow shock, as conjectured in previous works. Further, the new Vlasiator capabilities describing multiple ion populations were harnessed in Battarbee et al. (2020a) showing that foreshock properties can be also modified by specular reflection of $He^{2+}$ at the bow shock.
## 5.9.2 Terrestrial magnetosheath

Vlasiator also offers unprecedented insight into the physics of the magnetosheath, which is the region of shocked solar wind plasma surrounding the magnetosphere, and is constantly influenced by the upstream foreshock dynamics. Turc et al. (2020) show that the observed large-scale asymmetries in the magnetosheath are reproduced with Vlasiator. This study also highlights the strong variability of magnetosheath properties due to ion-kinetic physics, which are missing from the previous state-of-the-art simulations modelling the magnetosheath with the MHD approximation. These results also emphasise the role of foreshock-bow shock-magnetosheath interaction in determining the magnetosheath characteristics. An important topic that has recently gained much attention is the magnetosheath high-speed jets that originate at the bow shock and travel through the magnetosheath, sometimes impacting the magnetopause (Plaschke et al. 2018). While the origin of magnetosheath jets is still under debate, Vlasiator simulations initially indicated (Palmroth et al. 2018b) and later confirmed (Suni et al. 2021) the presence of a new formation mechanism. These works showed that jets can be directly related to steepened foreshock ULF waves impacting the bow shock like bullets that then travel through the magnetosheath. Finally, Vlasiator has also been used to identify a novel mechanism that modifies magnetosheath dynamics. Using a 2D simulation, Jarvinen et al. (2018) demonstrate that large-scale magnetosheath bow waves discovered by Pfau-Kempf et al. (2016) accelerate ions. These waves originate from the interaction between the magnetosheath and the magnetopause, thereby altering the properties of the magnetosheath through physical processes occurring downstream of it. More recently, Tesema et al. (2024) use a 3D-3V simulation to further reveal that the FTE-driven waves excite Pc2-range waves within the magnetosphere. This findings highlight ways through which magnetopause reconnection influences both the magnetosheath and the magnetosphere.

## 5.9.3 ULF wave transmission

One of the long-standing puzzles in space physics is how waves generated in the foreshock are transmitted into the magnetosphere, where they are observed as so-called Pc3 waves (Turc et al. 2025). The main issue is how these waves can survive through the bow shock, as there is no conclusive evidence for the existence of such waves in the magnetosheath. While previous numerical simulations using a local setup indicated that the waves should convert from fast-mode waves to Alfvén waves through the bow shock, preventing their transmission into the magnetosphere, Vlasiator revealed fast-mode signals propagating from the foreshock into the magnetosheath (Turc et al. 2023). The magnetosheath waves are caused by a modulation of the shock properties by the foreshock waves, which in turn send modulated disturbances through the magnetosheath which preserve the wave period that is thus transmitted all the way to the magnetosphere. The simulation results were confirmed by data from the MMS mission (Turc et al. 2023). Foreshock wave transmission was also investigated by Takahashi et al. (2021) and Turc et al. (2022), revealing that the shorter spatial scales associated with these waves during high IMF

<figure><img src="image_11.png" /></figure>
conditions hinder their propagation into the magnetosphere. Recent developments to Vlasiator have also enabled running simulations with time-varying solar wind conditions, including, for example, ULF density fluctuations, which were shown to launch fast-mode signals through the magnetosheath that alter the magnetosheath wave properties (Zhou et al. 2022).

### 5.9.4 Magnetic reconnection at the terrestrial magnetopause and magnetotail

Magnetic reconnection is a topic of intensive research, as it is the main process through which plasma and energy are transferred from the solar wind into the magnetosphere. Ala-Lahti et al. (2022) demonstrate using Vlasiator that magnetopause reconnection and the resulting FTEs modify the energy transferred through the magnetopause, and especially in the nightside magnetopause they have a significant role in energy injection to the magnetosphere. Akhavan-Tafti et al. (2020) found that the FTEs created by dayside magnetopause reconnection are similar to those measured by the MMS satellite observations. Pfau-Kempf et al. (2020) used a azimuthally restricted 3D-3V setup of Vlasiator and confirmed the earlier results regarding dayside reconnection in 2D-3V, summarised by Palmroth et al. (2018a), indicating that magnetosheath waves and FTEs modify reconnection locations at the magnetopause. Magnetopause reconnection was also investigated using the new eVlasiator model (Battarbee et al. 2021, and Sect. 5.8). Alho et al. (2022) found that the eVlasiator electron distribution functions near magnetopause reconnection sites were remarkably similar to the ones measured by the MMS satellite constellation. This is a somewhat surprising result as the eVlasiator model does not include electron dynamic effects to the electromagnetic field, but is more a test electron representation of the electron physics. Hence, Alho et al. (2022) results indicate that electron physics is largely controlled by ion-scale dynamics and ion-driven magnetic fields even in the vicinity of reconnection sites.

Before moving to full 3D-3V simulations of the magnetosphere, reconnection in the magnetotail was also studied in 2D-3V setups that are worth mentioning. Juusola et al. (2018) report a general agreement between Vlasiator tail reconnection and satellite observations, and that the velocity of the reconnection points is determined by the local characteristics of the reconnection electric field. Runov et al. (2021) investigate ion distribution functions within and near magnetotail reconnection sites, finding overall agreement between the Vlasiator results and spacecraft observations. The study identifies multiple distinct ion-kinetic signatures in the simulated ion distributions, which can be considered as characteristic 'fingerprints' of ion-kinetic dynamics. Specifically, they report the presence of gyrotropic and counterstreaming ion populations, indicative of betatron and Fermi acceleration. Additionally, they observe distributions resulting from Speiser orbits and D-shaped populations, all of which are consistent with observational studies.

### 5.9.5 Results based on global 3D-3V simulations

The new 3D-3V runs enabled unprecedented studies of reconnection in the magnetotail. Palmroth et al. (2023) presented a breakthrough in space physics by
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Evolution of the dynamic magnetotail current sheet in a 3D-3V Vlasiator simulation, featuring in particular magnetic reconnection and the current sheet kink instability. Data are shown on the plasma sheet centre, where colour coding shows current density. Yellow lines indicate positions where plasma flow changes sign (reconnection flow reversal). Fuchsia points show reconnection X-points, and green points show magnetic field O-topologies, which are also further illustrated with light green lines displaying the twisted magnetic field lines of the plasmoids. The image sequence shows how a large scale plasmoid is released in the tail. Figure from Palmroth et al. (2023)</figcaption></figure>

showing for the first time that both reconnection and ion-kinetic instabilities are involved in the large-scale topological reconfiguration of the magnetotail that ejects a plasma structure called a plasmoid (see green structures in Fig. 9). Previous studies did not consider that these processes could operate simultaneously, as spacecraft observations provide only localized measurements, making it difficult to capture the cross-scale sequence of events in the magnetotail. Additionally, earlier state-of-the-art MHD models reproduced reconnection but did not account for ion-kinetic instabilities. Palmroth et al. (2023) show that ion-kinetic instabilities lead to large-scale flapping of the magnetotail, which then disrupts the central plasma sheet
current, while reconnection changes the large-scale topology. In addition, Alho et al. (2024) present an important method to find reconnection locations in 3D. This method is a basis for more in-depth investigations of 3D reconnection in the magnetosphere. For example, Cozzani et al. (2025) reveal the intricate interplay between reconnection and kink instability in the magnetotail by showing that the reconnection rate is significantly influenced by the instability. In particular, the reconnection rate is reduced during the growth phase of the instability. Further, Zaitsev et al. (2025) have shown that demagnetized protons in the thin current layer in the magnetotail drive the development of tearing and kink instabilities. These results show that magnetotail dynamics are far more complex than previously thought based on global MHD models, and more studies and spacecraft missions are required to fully understand the fundamental workings of the magnetosphere.

A recent development in Vlasiator is the computation of the precipitating ion differential number flux at locations magnetically connected to the ionosphere. This operation is applied directly to the velocity distribution and yields the energy-dependent number flux of magnetospheric ions bound to precipitate into the ionosphere. These ions are considered to precipitate into the ionosphere because their pitch angle is too low to enable mirroring at a high enough altitude above the atmosphere. The calculation of the differential number flux of precipitating ions is described in detail by Grandin et al. (2019). Initial studies in 2D-3V were carried out at selected locations in the noon–midnight meridional plane to investigate nightside (Grandin et al. 2019) and dayside (Grandin et al. 2020) proton precipitation. Later on, 3D-3V runs made it possible to examine the precipitating particle fluxes at any location magnetically connected to the Earth. Grandin et al. (2023) investigated the spatial extent and temporal variations of auroral (0.5–50 keV) proton precipitation in the polar cusps and in the nightside auroral oval during southward IMF driving. Comparison of the obtained fluxes with observations from the Defense Meteorological Satellite Program (DMSP) Special Sensor J (SSJ) particle telescope during similar driving conditions revealed good agreement between the simulation and the measured data.

As an important new topic, Vlasiator was used to investigate soft X-ray emissions associated with solar wind charge exchange processes involving interactions between multiply charged heavy ions (such as $O^{7+}$ and $O^{8+}$) and neutral atoms (Sibeck et al. 2018). In view of the upcoming Solar wind Magnetosphere Ionosphere Link Explorer (SMILE) and Lunar Environment heliospheric X-ray Imager (LEXI) satellite missions designed to image near-Earth space in soft X-rays, Grandin et al. (2024) evaluated the detectability of the dayside magnetopause, magnetosheath mirror-mode waves, and flux transfer events in soft X-ray images using outputs from a Vlasiator 3D-3V run. These results contribute to better anticipate the future observations of the SMILE and LEXI missions, and they also give a reference point about how the observations should be understood in terms of global ion-kinetic physics.
## 5.10 Future avenues

Major building blocks for making Vlasiator possible in the past were not only the increase of the computational resources, but also several algorithmic innovations following the upcoming technologies. Examples of these are the sparse grid for distribution functions (see Sect. 5.2.2), and the semi-Lagrangian solver (see Sect. 4.4). Further, the code has been continuously optimised to fit better in different new parallel architectures. With these main steps, the efficiency of Vlasiator has been improved by orders of magnitude relative to the performance at the beginning of the project, allowing to simulate first 2D-3V systems with high resolution (Palmroth et al. 2017), and then transition to 3D-3V systems (Palmroth et al. 2023), and implement the coupling to the ionosphere (Ganse et al. 2025). These achievements were essentially made possible by adaptive mesh refinement (Kotipalo et al. 2024) in ordinary space, implementing higher resolution in regions of interest, and consequently reducing computational cost and data usage allowing for a greater number of simulated and stored time steps (see Sect. 5.2.1).

The development of Vlasiator is closely tied to the awarded grants. The near-term goals include pushing the inner boundary inwards from the current 30,000 km closer to 15,000 km to host the entire region of the ring current, and to increase the run duration towards several hours. This requires improving the computational efficiency of the already highly efficient parallel algorithms. The main avenues to achieve this are to add ‘time classes’ to allow longer simulation time steps in regions of weak magnetic fields, to improve the use of GPU resources, and to harness machine learning methods to increase runtime efficiency. One very promising avenue is to decrease the size of the checkpoint files by compressing and then reconstructing the VDFs using machine learning methods. In the field of weather modelling, Kurth et al. (2023) were the first to demonstrate that neural networks trained on surface and atmospheric variables can be used to produce forecasts with accuracy comparable to current state-of-the-art numerical weather prediction models. These results were soon followed by high-impact papers (Bi et al. 2023; Lam et al. 2023), in which further increases in accuracy were presented. These achievements are possible in the field of space physics as well.

In terms of physics, the most important addition was the coupling to the ionosphere. In the near term, the ionospheric boundary will be further improved. To enable the evaluation of the influence of the ion-kinetic magnetosphere on the ionosphere, the first coupling scheme implemented in Vlasiator is based on the type of boundary used in global MHD simulations, which couples the field-aligned currents, electric potential and effects of electron precipitation between the ionosphere and magnetosphere (e.g., Goodman 1995; Janhunen et al. 2012). However, it was noted (Ganse et al. 2025) that, since in the current scheme electron precipitation is based on the ion-kinetic magnetosphere and hence on the properties of the protons, the ionospheric mapping of the precipitation does not reliably reproduce electron precipitation. Hence, the next step is to improve the ionospheric conductivity to make it more representative of the patterns resulting from electron precipitation in nature. A future goal is to investigate inner magnetospheric physics in terms of solar wind driving, which would close the circle: The problems associated
with reproducing inner magnetospheric physics by global MHD simulations was one of the main motivations for developing Vlasiator in the first place.

Very recently, Vlasiator has also been applied to other environments such as Mars and Mercury, which are smaller objects with less intense (or nearly non-existent) magnetic fields. Hence, they are computationally less demanding to simulate. However, the inner boundary involves considerable differences as in Mercury there is no ionosphere, while in Mars there is mass loading to the exosphere. Cometary environments and comet-solar wind interactions should be interesting in terms of the available heavy ion support, in the context of the Rosetta and Comet Interceptor mission data analysis. Further, the Juice mission visiting the icy moons of Jupiter leads to hypothesise that e.g., the Ganymede-Jupiter interaction may also be one viable option for study in the future.

# 6 Conclusions and outlook

There are several main conclusions that can be made from the progress in the use of Vlasov-based methods in space plasma physics so far. The first one is related to the applicability of the ion-kinetic hybrid-Vlasov approach within the global magnetospheric context. When Vlasiator was first proposed, concerns arose as to whether ions are the dominant species controlling the system dynamics or does one need electrons as well. In particular, a physical representation of the reconnection process may require electrons, while the ion-kinetic Vlasiator approach would still model reconnection similarly as global MHD simulations, namely by enabling it via numerical diffusion. However, even an MHD simulation, treating both ions and electrons as a fluid, is capable of modelling global magnetospheric dynamics (Palmroth et al. 2006), indicating that the reconnection driving global dynamics must be within the right ballpark. Since Vlasiator produces results that are in agreement with *in situ* measurements (e.g., Runov et al. 2021; Grandin et al. 2023), ion-kinetic effects seem to be a major contributor in driving global dynamics. eVlasiator results, on the other hand, show that even in the vicinity of reconnection sites, ion-based electromagnetic fields are able to drive realistic electron velocity distribution functions (Alho et al. 2022). Hence, the leap from MHD to an ion-kinetic magnetospheric simulation allowed to take into account new physical mechanisms arising from ion-kinetic dynamics without compromising the accuracy of electromagnetic fields by simplifying electron physics. This is a major leap in understanding, which was not self-evident at the beginning of the project.

Another major conclusion based on Vlasiator results is the role of grid resolution in global setups. Again, one of the largest concerns in the beginning of Vlasiator development was that the ion gyro-scales could not be resolved within a global simulation volume, raising fears that the outcomes would be MHD-like, even though early hybrid-PIC simulations also used a spatial resolution comparable to the solar wind ion inertial length which was sufficient for ion-kinetic effects to emerge (e.g., Omidi et al. 2005). In this context, the first Vlasiator runs included an element of surprise, as even rather coarse resolution grids produced kinetic phenomena that were in agreement with *in situ* observations (Pokhotelov et al. 2013). The latest results
have clearly indicated that kinetic physics emerges even at coarse spatial resolution (Pfau-Kempf et al. 2018; Palmroth et al. 2023). This is essentially due to two facts: 1) most ion-kinetic phenomena that have implications at the global scale originate from the higher-energy ions that have large gyroradii and are well-represented by a spatial resolution that might not resolve the gyroradii of thermal ions; and 2) because the VDFs are well-described by the choice to model the distribution functions and not by reconstruction from macroparticles. It should be emphasised that this result could not have been foreseen before developing the model. Further, it indicates that also some electron physics could be evaluated without resolving the actual electron spatial scales, as demonstrated with the eVlasiator approach (Battarbee et al. 2021; Alho et al. 2022). One can hence conclude that others attempting to develop a (hybrid-) Vlasov simulation may face less concerns due to the grid resolution, even in setups with major computational challenges like, for example, portions of the Sun.

The most common physical conclusion based on Vlasiator simulations is that ‘everything affects everything’, indicating that scale coupling is important in global magnetospheric dynamics. A general avenue of development for global MHD simulations has been code coupling, where for example problem-specific codes have been coupled into the global context (e.g., Huang et al. 2006; Poedts et al. 2020; Gombosi et al. 2021), or, for example, full-PIC simulations have been embedded within an MHD domain (Tóth et al. 2016; Gombosi et al. 2023). While these approaches are interesting and advance physical understanding, they cannot approach scale coupling as specific kinetic phenomena are only addressed within their respective simulation volumes. A prime example of scale coupling is instability and reconnection acting in different parts of the magnetotail to cause the plasmoid eruption (Palmroth et al. 2023). These results could not have been achieved without a simulation that resolves both small and large scales simultaneously in the same simulation volume.

Lastly, an almost completely unexplored regime awaits in the application of relativistic Vlasov methods to plethora of astrophysical applications. Vlasov-based methods have not yet been widely adopted in the fields of astrophysics and space physics for modelling large-scale systems beyond the few examples cited in Tables 2 and 3, mainly due to the truly astronomical computational cost such simulations can incur. Some interesting systems to explore include, for example, the plasma dynamics of neutron-star magnetospheres, where high momentum-space accuracy is needed to resolve the existence of some speculated beam instabilities (e.g., Lyubarsky and Eichler 2007). Another avenue is the combination of radiation mechanisms and quantum-electrodynamic (QED) processes to simulate plasmas embedded in high-intensity radiation fields. Recently, such systems—in the context of black-hole surroundings—have, for example, been simulated with hybrid methods, where the radiation field is captured with the Vlasov technique for some QED processes, but the plasma dynamics are still solved with traditional particle-in-cell methods (Nättilä 2024). In summary, the Vlasov technique is also starting to take a foothold in astrophysics and is actively developed, for example, within the Runko framework (Nättilä 2022).

We emphasise that a critical success factor in the development of Vlasiator has been the close involvement with technological advances in the field of high-
performance computing. European research infrastructures for supercomputing have been developed almost hand-in-hand with Vlasiator, providing an opportunity to always target the newest platforms, thus feeding directly into code development. Should similar computationally intensive codes be designed and implemented elsewhere, it is recommended to keep a keen eye on the technological development of the supercomputing platforms. The experience with Vlasiator demonstrates that Vlasov-based modelling is complementary to other methods and provides unprecedented insights that are well worth the implementation effort.

**Acknowledgements** We acknowledge The European Research Council for Starting grant 200141-QuESpace, with which Vlasiator (https://www.helsinki.fi/en/researchgroups/vlasiator) was developed, and Consolidator grant 682068-PRESTISSIMO awarded to further develop Vlasiator and use it for scientific investigations. We gratefully also acknowledge the Research Council (formerly Academy) of Finland (grant numbers 361901, 339327, 335554, 347795, 345701). The Finnish Centre of Excellence in Research of Sustainable Space, funded through the Research Council of Finland with grant number 352846, supports Vlasiator development and science as well. Vlasiator parallel architecture is developed in the European Centre of Excellence in Code Development PLASMA-PEPSC, grant number 4100455, and we acknowledge the Research Council of Finland matching funding (grant number 359806). AI methods are developed in a EuroHPC Inno4Scale project ASTERIX, grant number 4100534. We acknowledge all computational grants we have received: PRACE/Tier-0 2012061111 on Hermit/HLRS, PRACE/Tier-1 on Abel/UiO-NOTUR, PRACE/Tier-0 2014112573 on HazelHel/HLRS, PRACE/Tier-0 2016153521 on Marconi/CINECA, PRACE/Tier-0 2019204998 on Hawk/HLRS; EHPC-REG-2022R02-238/LUMI-C and EHPC-EXT-2023E02-060/LUMI-C. CSC - IT Center of Science Grand Challenge grants in 2015 and 2016, 2017, 2018, 2020, 2021 and 2023, the pilot use in summer as well as the special Christmas present pilot use of sisu.csc.fi in 2014 are gratefully acknowledged. The authors thank the Finnish Computing Competence Infrastructure (FCCI), the Finnish Grid and Cloud Infrastructure (FGCI) and the University of Helsinki IT4SCI team for supporting this project with computational and data storage resources.

YPK acknowledges Research Council of Finland grant 339756 (KIMCHI). MG acknowledges Research Council of Finland grants 338629 (AERGELC'H) and 360433 (ANAON). The work of LT is funded by the European Union (ERC grant WAVESTORMS - 101124500). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. LT also acknowledges support from the Research Council of Finland (grant number 322544). GC is supported by the Integration Fellowship of Le Studium Loire Valley Institute for Advanced Studies.

**Author contributions** All authors contributed in generating the structure of the manuscript. MP contributed to all sections and consolidated the inputs from other authors to get the final manuscript. UG, YPK, MB, MA, KP, LK contributed to Sects. 4 and 5. MA contributed to Sects. 2, 4 and 5. IZ contributed to Sects. 3. GC, LT, SH, MG contributed to Sects. 2 and 5. HZ contributed to Sects. 5. LP contributed to Sects. 4. AS and SVA read through the manuscript and commented. The authors acknowledge the work of Thiago Brito in the 2018 edition of this review.

**Data availability** No datasets were generated or analysed during the current study.

## Declarations

**Conflict of interest** The authors declare no conflict of interest.
**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

## References

Afanasiev A, Battarbee M, Vainio R (2015) Self-consistent Monte Carlo simulations of proton acceleration in coronal shocks: effect of anisotropic pitch-angle scattering of particles. Astron Astrophys 584:A81. https://doi.org/10.1051/0004-6361/201526750. arXiv:1603.08857 [physics.space-ph]

Afanasiev A, Vainio R, Rouillard AP et al (2018) Modelling of proton acceleration in application to a ground level enhancement. Astron Astrophys. https://doi.org/10.1051/0004-6361/201731343

Aguilar M, Aisa D, Alpat B et al (2015) Precision measurement of the proton flux in primary cosmic rays from rigidity 1 GV to 1.8 TV with the Alpha Magnetic Spectrometer on the International Space Station. Phys Rev Lett 114(17):171103. https://doi.org/10.1103/PhysRevLett.114.171103

Akhavan-Tafti M, Palmroth M, Slavin JA et al (2020) Comparative analysis of the Vlasiator simulations and MMS observations of multiple X-line reconnection and flux transfer events. J Geophys Res Space Phys 125(7):e2019JA027410. https://doi.org/10.1029/2019JA027410

Ala-Lahti M, Pulkkinen TI, Pfau-Kempf Y et al (2022) Energy flux through the magnetopause during flux transfer events in hybrid-Vlasov 2D simulations. Geophys Res Lett 49(19):e2022GL100079. https://doi.org/10.1029/2022GL100079

Alho M, Battarbee M, Pfau-Kempf Y et al (2022) Electron signatures of reconnection in a global eVlasiator simulation. Geophys Res Lett 49(14):e2022GL098329. https://doi.org/10.1029/2022GL098329

Alho M, Cozzani G, Zaitsev I et al (2024) Finding reconnection lines and flux rope axes via local coordinates in global ion-kinetic magnetospheric simulations. Ann Geophys 42(1):145–161. https://doi.org/10.5194/angeo-42-145-2024

Allmann-Rahn F, Lautenbach S, Grauer R (2022) An energy conserving Vlasov solver that tolerates coarse velocity space resolutions: simulation of MMS reconnection events. J Geophys Res Space Phys 127 (2):e2021JA029976. https://doi.org/10.1029/2021JA029976

Allmann-Rahn F, Lautenbach S, Deisenhofer M et al (2024) The muphyII code: multiphysics plasma simulation on large HPC systems. Comput Phys Commun 296:109064. https://doi.org/10.1016/j.cpc.2023.109064

André M, Vaivads A, Khotyaintsev YV et al (2010) Magnetic reconnection and cold plasma at the magnetopause. Geophys Res Lett 37:L22108. https://doi.org/10.1029/2010GL044611

Anekallu CR, Palmroth M, Pulkkinen TI et al (2011) Energy conversion at the Earth's magnetopause using single and multispacecraft methods. J Geophys Res Space Phys 116(A15):A11204. https://doi.org/10.1029/2011JA016783

Angelopoulos V (2008) The THEMIS mission. Space Sci Rev 141:5–34. https://doi.org/10.1007/s11214-008-9336-1

Arber T, Vann R (2002) A critical comparison of Eulerian-grid-based Vlasov solvers. J Comput Phys 180 (1):339–357. https://doi.org/10.1006/jcp.2002.7098

Arrò G, Califano F, Pucci F et al (2024) Large-scale linear magnetic holes with magnetic mirror properties in hybrid simulations of solar wind turbulence. Astrophys J Lett 970(1):L6. https://doi.org/10.3847/2041-8213/ad61da

Arrò G, Califano F, Lapenta G (2020) Statistical properties of turbulent fluctuations associated with electron-only magnetic reconnection. Astron Astrophys 642:A45. https://doi.org/10.1051/0004-6361/202038696

Asahi Y, Latu G, Grandgirard V et al (2020) Performance portable implementation of a kinetic plasma simulation mini-app. In: Wienke S, Bhalachandra S (eds) Accelerator programming using directives
(WACCPD 2019). Lecture Notes in Computer Science, vol 12017. Springer, Cham, pp 117–139.
https://doi.org/10.1007/978-3-030-49943-6_5

Aunai N, Hesse M, Black C et al (2013) Influence of the dissipation mechanism on collisionless magnetic reconnection in symmetric and asymmetric current layers. Phys Plasmas 20(4):042901. https://doi.org/10.1063/1.4795727

Aunai N, Smets R, Ciardi A et al (2024) PHARE: parallel hybrid particle-in-cell code with patch-based adaptive mesh refinement. Comput Phys Commun 295:108966. https://doi.org/10.1016/j.cpc.2023.108966

Axford WI, Leer E, Skadron G (1977) The Acceleration of Cosmic Rays by Shock Waves. In: international cosmic ray conference, pp 132–137

Baker DN (1995) The inner magnetosphere: a review. Surv Geophys 16(3):331–362. https://doi.org/10.1007/BF01044572

Balbus SA, Hawley JF (1998) Instability, turbulence, and enhanced transport in accretion disks. Rev Mod Phys 70(1):1–53. https://doi.org/10.1103/RevModPhys.70.1

Ballerini G, Lavorenti F, Califano F et al (2024) Whistler-mode waves in the tail of Mercury's magnetosphere: a numerical study. Astron Astrophys 687:A204. https://doi.org/10.1051/0004-6361/202349093

Balogh A, Treumann RA (2013) Physics of collisionless shocks: space plasma shock waves. ISSI Scientific Report Series, vol 12. Springer, New York. https://doi.org/10.1007/978-1-4614-6099-2

Balsara DS (2009) Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics. J Comput Phys 228(14):5040–5056. https://doi.org/10.1016/j.jcp.2009.03.038

Balsara DS (2017) Higher-order accurate space-time schemes for computational astrophysics–Part I: finite volume methods. Living Rev Comput Astrophys 3(1):2. https://doi.org/10.1007/s41115-017-0002-8

Balsara DS, Kim J (2004) A comparison between divergence-cleaning and staggered-mesh formulations for numerical magnetohydrodynamics. Astrophys J 602(2):1079. https://doi.org/10.1086/381051

Banks JW, Hittinger JAF (2010) A new class of nonlinear finite-volume methods for Vlasov simulation. IEEE Trans Plasma Sci 38(9):2198–2207. https://doi.org/10.1109/TPS.2010.2056937

Battarbee M, Blanco-Cano X, Turc L et al (2020a) Helium in the Earth's foreshock: a global Vlasiator survey. Ann Geophys 38(5):1081–1099. https://doi.org/10.5194/angeo-38-1081-2020

Battarbee M, Ganse U, Pfau-Kempf Y et al (2020b) Non-locality of Earth's quasi-parallel bow shock: injection of thermal protons in a hybrid-Vlasov simulation. Ann Geophys 38(3):625–643. https://doi.org/10.5194/angeo-38-625-2020

Battarbee M, Brito T, Alho M et al (2021) Vlasov simulation of electrons in the context of hybrid global models: an eVlasiator approach. Ann Geophys 39(1):85–103. https://doi.org/10.5194/angeo-39-85-2021

Battarbee M, Alho M, Pfau-Kempf Y et al (2024) Analysator. https://doi.org/10.5281/zenodo.4462514, https://github.com/fmihpc/analysator

Battarbee M, Papadakis K, Ganse U et al (2025) Porting the grid-based 3D+3V hybrid-Vlasov kinetic plasma simulation Vlasiator to heterogeneous GPU architectures. J Phys: Conf Ser 2997(1):012010. https://doi.org/10.1088/1742-6596/2997/1/012010, arXiv:2406.02201

Bauer S, Kunze M (2005) The darwin approximation of the relativistic Vlasov-Maxwell system. Ann Henri Poincaré 6(2):283–308. https://doi.org/10.1007/s00023-005-0207-y

Becerra-Sagredo J, Málaga C, Mandujano F (2016) Moments preserving and high-resolution semi-Lagrangian advection scheme. SIAM J Sci Comput 38(4):A2141–A2161. https://doi.org/10.1137/140990619

Beckingsale DA, Burmark J, Hornung R, et al (2019) RAJA: Portable performance for large-scale scientific applications. In: IEEE/ACM international workshop on performance, portability and productivity in HPC (P3HPC), pp 71–81. https://doi.org/10.1109/P3HPC49587.2019.00012

Behar E, Fatemi S, Henri P et al (2022) Menura: a code for simulating the interaction between a turbulent solar wind and solar system bodies. Ann Geophys 40(3):281–297. https://doi.org/10.5194/angeo-40-281-2022

Bell AR (1978) The acceleration of cosmic rays in shock fronts – I. MNRAS 182:147–156. https://doi.org/10.1093/mnras/182.2.147

Berenger JP (1994) A perfectly matched layer for the absorption of electromagnetic waves. J Comput Phys 114(2):185–200. https://doi.org/10.1006/jcph.1994.1159