Simulations of the coupled ionosphere–magnetosphere system are a key tool to understand geospace and its response to space weather. For the most part, they are based on fluid descriptions of plasma (magnetohydrodynamics, MHD) formalism, coupled to an electrostatic ionosphere. Kinetic approaches to modeling the global magnetosphere with a coupled ionosphere system are still a rarity. We present an ionospheric boundary model for the global near-Earth plasma simulation system Vlasiator. It complements the magnetospheric hybrid-Vlasov simulations with an inner boundary condition that solves the ionospheric potential based on field-aligned current and plasma quantities from the magnetospheric domain. This new ionospheric module solves the ionospheric potential in a height-integrated approach on an unstructured grid and couples back to the hybrid-kinetic simulation by mapping the resulting electric field to the magnetosphere's inner boundary. The solver is benchmarked against a set of well-established analytic reference cases, and we discuss the benefits of a spherical Fibonacci mesh for use in ionospheric modeling. Preliminary results from coupled global magnetospheric–ionospheric simulations are presented, showing formation of both Region 1 and Region 2 current systems.

Bursty Bulk Flows (BBFs) are transient plasma flows in the Earth's magnetotail plasma sheet. These short-lived, high-speed flows play a key role in the magnetosphere-ionosphere coupling. Currently, most insights into the ionospheric signatures of BBFs come from individual case studies that include conjugate observations of BBFs in the magnetotail and field-aligned currents (FACs) in the nightside ionosphere. In this study, we utilise the 6D hybrid-Vlasov simulations to study the ionospheric signatures of BBFs in the near-Earth magnetotail. We show that a BBF with Vx≥400 km s−1 emerges shortly after magnetic reconnection occurs on the duskside at a radial distance between 11 and 14 RE (where RE=6371 km is the radius of the Earth) in the current sheet. As the BBF moves Earthward, clockwise (counterclockwise) flow vortices are induced on its dawn (dusk) sides. These vortical flows generate FACs flowing upward (out of the current sheet) on the dawnside and downward (into the current sheet) on the duskside flank, respectively. The mapping of BBF structures onto the ionosphere shows that the structure is primarily aligned in the east-west direction, with its ionospheric signatures appearing as enhancements in FACs, ionospheric conductances, horizontal ionospheric currents, energies of precipitating electrons and protons, and the formation of localised plasma flow channels. The upward and downward FACs associated with BBFs in the magnetotail consistently map to enhanced Region 2 (R2) and Region 1 (R1) FAC structures at ionospheric altitude, which are then closed in the ionosphere by north-west flowing Pedersen currents. The ionospheric counterpart of the Earthward plasma flow of the BBF is a channel of equatorward plasma flow, while the westward drift of these enhanced structures corresponds to the duskward motion of the BBF in the magnetotail.

Solar wind directional discontinuities can generate transient mesoscale structures such as foreshock bubbles and hot flow anomalies (HFAs) upstream of Earth's bow shock. These structures can have a global impact on the near-Earth environment, and understanding their formation conditions is crucial to evaluate their contribution to solar wind-magnetosphere coupling. Here we present the results of a global 2D hybrid-Vlasov simulation (with 3D electromagnetic fields) of the interaction of a rotational discontinuity with near-Earth space, performed with the Vlasiator model. The magnetic field rotates by 90 degrees from ortho-Parker spiral to Parker spiral orientation across the discontinuity. As the discontinuity enters the simulation domain, a foreshock bubble forms duskward of the Sun-Earth line, where the foreshock is initially located. Shortly after the discontinuity makes first contact with the bow shock at the subsolar point, we find that a structure with enhanced temperature and strongly deflected flows develops at the intersection of the discontinuity with the bow shock. This structure displays typical features of an HFA. However, HFA formation requires electric fields pointing towards the discontinuity on at least one side, a condition which is not initially met in our simulation. We demonstrate that the prior generation of the foreshock bubble provides the necessary conditions for HFA formation. We then investigate the evolution of both structures as the discontinuity travels antisunward, showing that the foreshock bubble signatures tend to weaken while the HFA grows. We also report a large-scale bow shock deformation, with the bow shock expanding several Earth radii outward of its initial position within the compressed edge of the foreshock bubble. Our results provide new clues regarding the formation and evolution of large-scale foreshock transients and their impact on the shock.

We explore the three-dimensional structure of ion-kinetic instabilities in a thin current layer using a hybrid-Vlasov simulation of the Earth's magnetosphere. The simulation shows the simultaneous growth of tearing and kinking instabilities, which develop in the Sun-Earth and dawn-dusk directions, respectively, within the magnetotail current sheet. The formation of flux ropes indicates the development of the tearing instability, while flapping-type cross-tail oscillations arise from the kink instability. We consider both instabilities as independent polarizations, albeit sharing a common source: demagnetized ions forming crescent-shape distributions at the center of the current layer. These oscillations exhibit spatiotemporal characteristics within the proton-scale range, featuring a growth time on the order of 40-80 proton gyroperiods and a wavelength of approximately 15-30 proton skin depths.

Parallelization is a necessity for large-scale simulations due to the amount of data processed. In this article we investigate different load balancing methods using Vlasiator, a global magnetospheric simulation as our case study. The theoretical basis for load balancing is the (hyper)graph partitioning problem, modeling simulation units as vertices and their data dependencies as edges. As it is an NP-hard problem, heuristics are necessary for dynamic runtime balancing. We consider first hypergraph partitioning via an algorithm called parallel hypergraph partitioner (PHG); this is done by partitioning a simplified grid and then attempting to optimize the solution on the finer grid. The second and third are the geometric methods of recursive coordinate bisection (RCB) and recursive inertial bisection (RIB). Finally we consider the method of Hilbert space filling curves (HSFC). The algorithm projects simulation cells along a Hilbert curve and makes cuts along the curve. This works well due to the excellent locality of Hilbert curves, and can be optimized further by choice of curve. We introduce and investigate six three-dimensional Hilbert curves in total. Our findings on runs of two different scales indicate the HSFC method provides optimal load balance, followed by RIB and PHG methods and finally by RCB. Of the Hilbert curves evaluated, the Beta curve outperformed the most commonly used curve by a few percent.

In this study we investigate the formation of magnetosheath jets before, during, and after the interaction between Earth's bow shock and a solar wind rotational discontinuity in a 2D ecliptic simulation run of the global magnetospheric hybrid-Vlasov model Vlasiator. Magnetosheath jets are transient enhancements of dynamic pressure downstream of collisionless shocks, and they have been observed in Earth's magnetosheath, the magnetosheaths of other planets, as well as the sheaths of interplanetary shocks. Rotational discontinuities (RD) are boundaries where the components of the magnetic field and velocity tangential to the boundary change abruptly, and they have been observed by spacecraft in the solar wind and in Earth's magnetosheath. Both spacecraft observations and previous simulation studies have shown that RDs interacting with the bow shock can generate dynamic pressure pulses in the magnetosheath.Studying magnetosheath jets is important because they have been shown to potentially have magnetospheric effects if impacting the magnetopause, and while travelling through the magnetosheath they can modify its properties. Statistical studies of simulations and spacecraft observations have shown that jets tend to form mainly at Earth's quasi-parallel bow shock, that is where the interplanetary magnetic field (IMF) direction is nearly parallel to the shock normal, but they have also been observed downstream of the quasi-perpendicular shock. By studying the formation and properties of jets at the quasi-parallel and quasi-perpendicular shock at different times in the simulation, we aim to shed light on the differences between jets forming at different parts of the shock, and during different stages of interaction between an RD and the bow shock.

Modelling the distribution of odd nitrogen (NOx) in the polar middle and upper atmosphere has proven to be a complex task. Firstly, its production by energetic electron precipitation is highly variable across a range of temporal scales from seconds to decades. Secondly, there are uncertainties in the measurement-based but simplified electron flux datasets that are currently used in atmosphere and climate models. The altitude distribution of NOx is also strongly affected by atmospheric dynamics on monthly timescales, particularly in the polar winter periods when the isolated air inside the polar vortex descends from the lower thermosphere to mesosphere and stratosphere. Recent comparisons between measurements and simulations have revealed strong differences in the NOx distribution, with questions remaining about the representation of both production and transport in models. Here we present for the first time a novel approach, where the electron atmospheric forcing in the auroral energy range (50 eV–50 keV) is derived from a magnetospheric hybrid-kinetic simulation with a detailed description of energy range and resolution, as well as spatial and diurnal distribution. These electron data are used as input in a global whole-atmosphere model to study the impact on polar NOx and ozone. We show that the magnetospheric electron data provide a realistic representation of the forcing, which leads to considerable impact in the lower thermosphere, mesosphere, and stratosphere. We find that during the polar winter the simulated auroral electron precipitation increases polar NOx concentrations up to 215 %, 59 %, and 7.8 % in the lower thermosphere, mesosphere, and upper stratosphere, respectively, when compared to no auroral electron forcing in the atmospheric model. These results demonstrate the potential of combining magnetospheric and atmospheric simulations for detailed studies of solar wind–atmosphere coupling.

When plasma and magnetic field discontinuities in the solar wind interact with Earth’s magnetic field, they can significantly alter the properties and dynamics of Earth’s bow shock, magnetosheath, and magnetopause. In this study, we investigate magnetosheath dynamic pressure enhancements generated by the interaction between a solar wind rotational discontinuity with Earth’s bow shock in a 2D simulation run of the global magnetospheric hybrid-Vlasov model Vlasiator. We find that as the discontinuity is advected into the bow shock, several fast-mode pulses associated with enhanced dynamic pressure are launched toward the Earth. In addition, the interaction between discontinuity and shock generates transient enhancements of dynamic pressure at the bow shock that move Earthward together with the discontinuity. We find that the fast-mode pulses are able to traverse the magnetosheath and disturb the magnetopause. This finding differs from the results of previous studies using 2D and 3D MHD simulations as well as spacecraft measurements, which concluded that magnetopause disturbances should be caused by the rotational discontinuity itself and the dynamic pressure enhancement associated with it.

Vlasiator is a space plasma simulation code which models near-Earth ion-kinetic dynamics in three spatial and three velocity dimensions. It is highly parallelized, modeling the Vlasov equation directly through the distribution function, discretized on a Cartesian grid, instead of the more common particle-in-cell approach. Modeling near-Earth space, plasma properties span several orders of magnitude in temperature, density, and magnetic field strength. In order to fit the required six-dimensional grids in memory, Vlasiator utilizes a sparse block-based velocity mesh, where chunks of velocity space are added or deleted based on the advection requirements of the Vlasov solver. In addition, the spatial mesh is adaptively refined through cell-based octree refinement. In this paper, we describe the design choices of porting Vlasiator to heterogeneous CPU/GPU architectures. We detail the memory management, algorithmic changes, and kernel construction as well as our unified codebase approach, resulting in portability to both NVIDIA and AMD hardware (CUDA and HIP languages, respectively). In particular, we showcase a highly parallel block adjustment approach allowing efficient re-ordering of a sparse velocity mesh. We detail pitfalls we have overcome and lay out a plan for optimization to facilitate future exascale simulations using multi-node GPU supercomputing.

Magnetic flux ropes are helical structures of magnetic field which form in a variety of magnetized plasmas. In near-Earth space, flux ropes are a manifestation of energy transfer at the magnetopause and in the magnetotail current sheet. We present a new method to detect magnetic flux ropes in large-scale simulations using only magnetic field line tracing. The method does not require prior identification of structures of interest such as current sheets or null lines and thus allows one to identify flux ropes of any size and orientation anywhere in the simulation domain. In this work, the new method is implemented in the hybrid-Vlasov model Vlasiator and demonstrated in global simulations of the terrestrial magnetosphere.We study the evolution of flux ropes forming during flux transfer events on the dayside magnetopause under a southward interplanetary magnetic field. It is found that flux ropes with an axial orientation along the dawn-dusk direction and propagating beyond the cusps will rapidly reconnect with the lobe magnetic field and vanish. In contrast, the flux ropes remaining near the equatorial plane and with an axial orientation along the flow direction - that is, tangential to the magnetopause - can maintain their structure and propagate tens of Earth radii down the tail in the absence of a reconnecting shear magnetic field component. These results are a step forward in the global characterization of flux ropes in and around the magnetosphere and may help in guiding the search for elusive far-tail flux ropes in satellite measurements.

AbstractMagnetic reconnection and current sheet kink instability often develop concurrently in current sheets, yet their dynamic interplay remains unclear. We investigate their interaction in the magnetotail of a 3D global magnetospheric hybrid‐Vlasov simulation. We identify the instability growth and saturation phase and estimate the evolution of the reconnection rate during the same interval. Our findings indicate that the reconnection rate decreases during the instability growth phase, especially at locations where the current sheet undergoes significant perturbations. These results highlight the intricate three‐dimensional relationship between reconnection and kink instabilities, suggesting that the kink instability plays a significant role in modulating the reconnection rate.

The Dungey cycle is a fundamental process governing large-scale plasma dynamics in the near-Earth space, traditionally examined through Magnetohydrodynamic (MHD) simulations and ionospheric observations. However, MHD models often oversimplify the complexities of driving dynamics and kinetic processes, while observational data tend to lack sufficient coverage. In this study, we utilize a hybrid-Vlasov simulation to investigate the Dungey cycle, and introduce a novel method for quantifying reconnection voltages in different Magnetic Local Time (MLT) sectors. This method is validated by comparing it with the ionospheric open flux change rate in the simulation. Our analysis identifies discrete azimuthal convection channels of closed field lines, clearly initiated by dayside reconnection and propagating to the nightside. These channels are prominent even during intervals of intense nightside reconnection. Notably, we observe that the effective length of dayside reconnection fluctuates, even under steady solar wind conditions. Our results reveal significant deviations from MHD theory, which predicts that plasma flows within the magnetosphere should follow flux tube entropy isocontours. Instead, we demonstrate that plasma flows near reconnection sites and at the terminators deviate from isentropic behavior, suggesting the presence of non-adiabatic processes in these regions. This study validates the representation of the Dungey cycle in the Vlasiator 3D simulation and enhances our understanding of global plasma convection. Future work should focus on identifying the kinetic processes that explain the deviations in the plasma convection with flux tube entropy isocontours between MHD theory and kinetic approach.

Magnetosheath jets are transient enhancements of dynamic pressure downstream of collisionless shocks. In Earth's magnetosheath they are mostly found downstream of the quasi‐parallel bow shock during steady solar wind and low interplanetary magnetic field (IMF) cone angle conditions, but they have also been observed in the quasi‐perpendicular magnetosheath and during different solar wind conditions. In this study we use a 2D simulation run of the global hybrid‐Vlasov model Vlasiator to investigate how the interaction between the bow shock and a solar wind rotational discontinuity influences the formation of magnetosheath jets. Separating the jets identified in the simulation based on formation site and time relative to the interaction between the discontinuity and the shock, we conduct a statistical study to find the characteristic properties of the different jet types. We find that jets forming at the quasi‐parallel shock are similar to each other regardless of the stage of the shock‐discontinuity interaction. Jets forming at the quasi‐perpendicular shock after the discontinuity has passed are small and short‐lived. The jets forming at the quasi‐perpendicular shock as the discontinuity impacts it, on the other hand, merge with each other into a large and long‐lived transient density enhancement that propagates deep into the magnetosheath together with the discontinuity, giving it the potential to be more geoeffective than the other types. This study sheds light on the properties of jets and jet‐like structures that form during non‐steady solar wind and IMF conditions, and the results can be of use when classifying similar events from spacecraft observations.

Magnetic reconnection in Earth's magnetotail is thought to create bursty bulk flows (BBFs), short-lived plasma bulk velocity enhancements in the magnetotail's central plasma sheet (CPS) region. Closely related to BBFs are dipolarization fronts (DFs), sudden increases in Bz, the magnetic field component aligned with Earth's magnetic dipole axis. Both phenomena affect energy distribution and flux transport in the magnetotail.We demonstrate novel methods of identifying BBFs and DFs in a 3D global magnetospheric simulation and present results for multiple case studies. BBFs and DFs are searched for in a simulation conducted using Vlasiator, a global magnetospheric hybrid-Vlasov code where ions are modeled using distribution functions and electrons are treated as a charge-neutralizing fluid. DFs are identified using a magnetic field time derivative threshold dBz /dt > 0.35 nT/s. BBFs are defined based on a velocity threshold, and they are studied on a case-by-case basis. Tailward DFs (anti-dipolarization fronts) are found at magnetic islands, while earthward DFs are mostly seen in finger-like structures of high earthward bulk velocity alongside BBFs. Signatures registered as BBFs in spacecraft view also originate due to moving reconnection locations and movement of the current sheet within the CPS while the reconnection outflow stays roughly constant. The results show that rapid Bz variations in the simulation have multiple sources, and similar satellite measurements of BBFs can arise from different physical phenomena. The findings may help with interpreting satellite observations in the magnetotail.

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.

High-performance computing is used for diverse simulations, some of which parallelize over the Message Passing Interface (MPI) with ease, whilst others may have challenges related to uniform balancing of computational load and communication between simulation domains. We introduce an alternative approach to solving advection equations, specifically in an application to solving the six-dimensional Vlasov equation for modelling space plasmas. Communicating larger ghost domains around the partition assigned to each MPI task and computing on these ghost cells allows for coalescing several discrete communication calls into one. This approach needs more overall data communication and computation, but provides interesting new avenues for the balancing of computational load between MPI tasks. We discuss this trade-off, how it may assist in developing other algorithmic improvements, and how the transition to heterogeneous CPU-GPU architectures may impact its usefulness.

Solar wind directional discontinuities can generate transient mesoscale structures such as foreshock bubbles and hot flow anomalies (HFAs) upstream of Earth's bow shock. These structures can have a global impact on near‐Earth space, so understanding their formation conditions is essential. We investigate foreshock transient generation at a rotational discontinuity using a global 2D hybrid‐Vlasov simulation. As expected, a foreshock bubble forms on the sunward side of the discontinuity. Later, when the discontinuity reaches the shock, new structures identified as HFAs develop, despite the initial discontinuity not being favorable to HFA formation. We demonstrate that the foreshock bubble provides the necessary conditions for their generation. We then investigate the evolution of the transient structures and the large‐scale bow shock deformation they induce. Our results provide new insights on the formation and evolution of foreshock transients and their impact on the shock.

Magnetic reconnection in Earth's magnetotail is thought to create bursty bulk flows (BBFs), short‐lived plasma bulk velocity enhancements in the magnetotail's central plasma sheet (CPS) region. Closely related to BBFs are dipolarization fronts (DFs), sudden increases in , the magnetic field component aligned with Earth's magnetic dipole axis. BBFs are involved in, for example, transport of energy and mass, and DFs contribute to for instance energy conversion processes and plasma acceleration. Both phenomena increase magnetic flux transport in the magnetotail. We demonstrate novel methods of identifying BBFs and DFs in 3D global magnetospheric simulations and present results for multiple case studies. We search for BBFs and DFs in a simulation conducted using the 3D global magnetospheric hybrid‐Vlasov code Vlasiator. DFs are identified using a time derivative threshold, whereas BBFs are defined based on a velocity threshold. The tailward DF (anti‐DF) identification criteria are fulfilled by tailward‐propagating flux ropes, while earthward‐propagating DFs are mostly seen in finger‐like structures of high earthward bulk velocity alongside BBFs. We show that detections of fast flows meeting the BBF criteria in virtual spacecraft time series also originate due to moving reconnection locations and movement of the current sheet within the CPS region, while the reconnection outflow stays roughly constant. The results show that rapid variations in the simulation have multiple sources, and similar satellite measurements of BBFs can arise from different physical phenomena. Our findings may help with interpreting satellite observations in the magnetotail.

Magnetic reconnection is a crucially important process for energy conversion in plasma physics, with the substorm cycle of Earth's magnetosphere and solar flares being prime examples. While 2D models have been widely applied to study reconnection, investigating reconnection in 3D is still, in many aspects, an open problem. Finding sites of magnetic reconnection in a 3D setting is not a trivial task, with several approaches, from topological skeletons to Lorentz transformations, having been proposed to tackle the issue. This work presents a complementary method for quasi-2D structures in 3D settings by noting that the magnetic field structures near reconnection lines exhibit 2D features that can be identified in a suitably chosen local coordinate system. We present applications of this method to a hybrid-Vlasov Vlasiator simulation of Earth's magnetosphere, showing the complex magnetic topologies created by reconnection for simulations dominated by quasi-2D reconnection. We also quantify the dimensionalities of magnetic field structures in the simulation to justify the use of such coordinate systems.