# Simulation Models for Exploring Magnetic Reconnection

Michael Shay, Subash Adhikari, Naoki Beesho, Joachim Birn, Jörg Büchner, Paul Cassak, Li-Jen Chen, Yuxi Chen, Giulia Cozzani, James Drake, et al.

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

## To cite this version:

Michael Shay, Subash Adhikari, Naoki Beesho, Joachim Birn, Jörg Büchner, et al. Simulation Models for Exploring Magnetic Reconnection. Space Science Reviews, 2025, 221, 10.1007/s11214-025-01210-5. insu-05281020

HAL Id: insu-05281020

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

Submitted on 24 Sep 2025

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_2.png" /></figure>
# Simulation Models for Exploring Magnetic Reconnection

Michael Shay¹ · Subash Adhikari²,¹ · Naoki Beesho³ · Joachim Birn⁴ · Jörg Büchner⁵ ·

Paul Cassak² · Li-Jen Chen⁶ · Yuxi Chen⁷ · Giulia Cozzani⁸,⁹ · James Drake¹⁰,¹¹ ·

Fan Guo¹² · Michael Hesse¹³ · Neeraj Jain¹⁴ · Yann Pfau-Kempf⁸,¹⁵ · Yu Lin¹⁶ ·

Yi-Hsin Liu¹⁷ · Mitsuo Oka¹⁸ · Yuri Omelchenko⁴,¹⁹ · Minna Palmroth⁸ · Oreste Pezzi²⁰ ·

Patricia H. Reiff²¹ · Marc Swisdak¹¹ · Frank Toffoletto²¹ · Gabor Toth⁷ · Richard A. Wolf²¹

Received: 11 July 2024 / Accepted: 6 August 2025 / Published online: 9 September 2025
© The Author(s) 2025

## Abstract

Simulations have played a critical role in the advancement of our knowledge of magnetic reconnection. However, due to the inherently multiscale nature of reconnection, it is impossible to simulate all physics at all scales. For this reason, a wide range of simulation methods have been crafted to study particular aspects and consequences of magnetic reconnection. This article reviews many of these methods, laying out critical assumptions, numerical techniques, and giving examples of scientific results. Plasma models described include magnetohydrodynamics (MHD), Hall MHD, Hybrid, kinetic particle-in-cell (PIC), kinetic Vlasov, Fluid models with embedded PIC, Fluid models with direct feedback from energetic populations, and the Rice Convection Model (RCM).

**Keywords** Plasma simulation · Magnetic reconnection · Plasma physics · Magnetosphere · Solar corona · Turbulence · Numerical methods

## 1 Introduction

Numerical computation has always played an important role in science. The term “computer” was used during the Renaissance to describe a person who performed mathematical calculations, and such computers were used extensively to calculate the positions of the planets. However, with the advent of digital computers last century, the role of such computation has exploded and revolutionized science in general. The study of magnetic reconnection has seen such a revolution in the last several decades as both numerical power has increased and numerical techniques have become more sophisticated.

Magnetic reconnection is considered a multiscale process because it allows physics that emerges at very small length and time scales to have global consequences in the system. A straightforward example of this large separation of scales is magnetic reconnection on the sun side of Earth’s magnetosphere. In this region, magnetic field lines are finally broken on a length scale of the order of 5 km which is the electron inertial length $d_e \equiv c/\omega_{pe}$. However, the dynamical effects of this breaking of field lines include driving global convection of the magnetosphere, a system spanning 100s of Earth radii ($R_E$) which is hundreds of thousands of $d_e$. A grid scale of about a $d_e$ over 100 Earth Radii requires about 100,000 spatial grid
points in only 1 dimension. Clearly, accurately resolving the physics breaking the frozen-in constraint while simulating global scales is impossible.

The impossibility of globally simulating the whole 3D system and resolving all scales has led to the generation of a wide range of simulation models, each of which has its own strengths and weaknesses. Through many decades of research, scientists have carefully crafted these models for the particular application or applications they are studying. Typically, the more realistic physics that is included in the simulation, the more computationally expensive it is. Studies of the basic physics of magnetic reconnection (Biskamp 1996) have very often used kinetic PIC simulations which include all relevant physics, but require a simplified geometry and boundary conditions. Global magnetospheric simulations include the complex boundaries associated with the solar wind and the ionosphere, but until recently were required to be fluid models due to the cost of including kinetic effects.

In this paper we will provide an overview of the primary simulation models that are currently being used to study magnetic reconnection. Please note, however, that the topic of plasma simulation is extremely complex and detailed and cannot be fully covered in a single book, much less a single article. To assist the reader, we have included a table of information for representative simulations codes associated with the types of models described in this article (See Table 1 at the end of this section). If the reader wishes to dive even deeper into a particular model, there are many references available, many of which are cited in the individual sections of this paper and in the table. There are also excellent books devoted to the subject (e.g., Büchner et al. 2003; Büchner 2023).

In the field of magnetic reconnection research, more than one system of units is used. As of this writing, one can generally say that scientists specializing in theory/simulation primarily use cgs units and scientist specializing in observational analysis use SI units. We have chosen as much as possible to use cgs units in this paper, although the section on the Rice Convection Model has been left in SI units. For an excellent description on how to convert units between cgs and SI, please see the *NRL Plasma Formulary* (Huba 1998).

In magnetic reconnection the diffusion region occurs in thin boundary layers where the physics changes, ultimately allowing magnetic topology to change (Liu et al. 2025a). Although not exhaustive, the new physics which emerges in the diffusion region can be characterized by examining Ohm's law, which comes from the fluid electron momentum equation.

$$E = -\frac{\mathbf{u} \times \mathbf{B}}{c} + \eta \mathbf{J} + \frac{\mathbf{J} \times \mathbf{B}}{n_e e c} + \frac{m_e}{e^2} \frac{\partial}{\partial t} \left( \frac{\mathbf{J}}{n_e} \right) - \frac{m_e}{e^2} \left( \frac{\mathbf{J}}{n_e e} \right) \cdot \nabla \left( \frac{\mathbf{J}}{n_e} \right) - \frac{1}{n_e e} \nabla \cdot \mathbf{p}_e, \quad (1)$$

where **E** and **B** are the electric and magnetic fields, **J** is the current density, **u** is the single fluid bulk flow velocity, **η** is the resistivity, **p**<sub>e</sub> is the electron pressure tensor, $n_e$ is the number density, **e** is the proton charge, $m_e$ is the electron mass and **c** is the speed of light. The first term on the right hand side of Eq. (1) is the ideal term, the second term is the resistive term, the third term is the Hall term, the fourth and fifth terms collectively represent electron inertia and the final term is due to the electron pressure tensor.

For the organization of the paper, we choose to move generally from fluid models to kinetic models as exemplified by terms on the right hand side of Eq. (1). We start with magnetohydrodynamics (MHD) and gradually increase in physical complexity until reaching fully kinetic simulations. We end with the Rice Convection Model (RCM), a widely used model for the inner magnetosphere, which acts as an inner boundary for magnetic field lines which are reconnecting in the magnetosphere. Section 2 describes MHD – first two terms. Section 3 describes Hall MHD including electron inertia – third and fourth terms. Section 4 describes Hybrid Simulations – generally also third and fourth terms. Section 5 describes
<table><caption>**Table 1** List of representative simulation codes discussed in this manuscript</caption><thead><tr><th>Code</th><th>Language</th><th>Parallelization</th><th>Reference</th><th>Public</th></tr><tr><td></td><td colspan="4">URL</td></tr></thead><tbody><tr><td>Regional MHD</td><td>Fortran 77</td><td>None</td><td>Birn et al. (2025)<br>https://doi.org/10.1029/2024JA033648</td><td>No</td></tr><tr><td>Global MHD (BATSRUS)</td><td>Fortran 90</td><td>MPI, OpenMP,<br>OpenACC</td><td>Tóth et al. (2012)<br>https://github.com/SWMFsoftware/BATSRUS</td><td>Yes</td></tr><tr><td>Test particle in MHD</td><td>Fortran 77</td><td>None</td><td>Birn et al. (2022)<br>https://doi.org/10.3389/fspas.2022.908730</td><td>No</td></tr><tr><td>Hall MHD/EMHD (F3D)</td><td>Fortran 90</td><td>MPI</td><td>Shay et al. (2004)<br>https://doi.org/10.1063/1.1705650</td><td>No</td></tr><tr><td>Hybrid PIC/EMHD (CHIEF)</td><td>Fortran 90,<br>C++</td><td>MPI</td><td>Muñoz et al. (2018)</td><td>No</td></tr><tr><td>Global Hybrid PIC (HYPERS)</td><td>Fortran 77,<br>C++</td><td>MPI</td><td>Omelchenko et al.<br>(2021a)<br>https://ccmc.gsfc.nasa.gov/models/HYPERS-Global~2021/</td><td>Yes<sup>1</sup></td></tr><tr><td>Kinetic PIC (P3D)</td><td>Fortran 90</td><td>MPI</td><td>Zeiler et al. (2002)<br>https://github.com/spudam/P3D-PLASMA-PIC</td><td>Yes</td></tr><tr><td>Kinetic PIC (VPIC)</td><td>C++</td><td>MPI / OpenMP /<br>Kokkos</td><td>Bowers et al. (2008),<br>Bird et al. (2021)<br>https://github.com/lanl/vpic / https://github.com/lanl/vpic-kokkos</td><td>Yes</td></tr><tr><td>Kinetic PIC (ExPIC)</td><td>Fortran 77</td><td>MPI</td><td>Bessho et al. (2019)<br>https://doi.org/10.1029/2019GL083397</td><td>No</td></tr><tr><td>Embedded Kinetic PIC (FLEKS)</td><td>C++</td><td>MPI</td><td>Chen et al. (2023)<br>https://github.com/SWMFsoftware/FLEKS</td><td>Yes</td></tr><tr><td>Kglobal</td><td>Fortran 90</td><td>MPI</td><td>Drake et al. (2019)</td><td>No</td></tr><tr><td>Hybrid Vlasov (HVM)</td><td>Fortran 90,<br>Fortran 77</td><td>MPI, OpenMP,<br>OpenACC, CUDA</td><td>Valentini et al. (2007)</td><td>No</td></tr><tr><td>Global Hybrid Vlasov (Vlasiator)</td><td>C++</td><td>MPI, OpenMP</td><td>Palmroth et al. (2018a),<br>Ganse et al. (2023)<br>https://github.com/fmihpc/vlasiator</td><td>Yes</td></tr><tr><td>Inner Magnetosphere (RCM)</td><td>Fortran 95</td><td>MPI</td><td>Toffoletto et al. (2003)</td><td>Yes<sup>2</sup></td></tr></tbody></table>

<footnote n="1">¹Available at NASA CCMC</footnote>

<footnote n="2">²By request</footnote>
kinetic particle-in-cell simulations – the physics of the fifth term is added, although it arises from the collective effects of many individual electrons. Section 6 describes embedding PIC codes into fluid models like MHD. Section 7 describes Kglobal, an MHD model which self consistently evolves energetic particles. Section 8 describes kinetic Vlasov models. And Section 9 describes the Rice Convection Model.

# 2 MHD

## 2.1 Equations of MHD

Magnetohydrodynamics (MHD) is the simplest fluid model used to study large scale plasma dynamics. MHD is based on the assumption that the characteristic length and time scales of the system under study are much larger than the length and time scales of the plasma species, usually Debye length ($\lambda_D$) or gyroradius and gyroperiod. Therefore, MHD represents the slow evolution of plasmas, often electrons and ions as a single fluid. The macroscopic behavior of the fluid in presence of a magnetic field is described by MHD using hydrodynamics and Maxwell's equations.

Let us consider a fluid (in this case a plasma with ions and electrons), moving with a flow velocity $\mathbf{u}$, characterized by a mass density $\rho = m_i n_i$ where $m_i$ is the mass of protons ($m_i \gg m_e$) and $n_i$ is the number density of protons (with quasi-neutrality $n_i = n_e = n$), thermal pressure $p$, and a magnetic field $\mathbf{B}$. The evolution of these fields in space and time are governed by the MHD equations given by

$$ \frac{\partial \rho}{\partial t} = -\nabla \cdot (\rho \mathbf{u}), \qquad (2) $$

$$ \rho \frac{\partial \mathbf{u}}{\partial t} + \rho(\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \frac{1}{4\pi} (\nabla \times \mathbf{B}) \times \mathbf{B} + \nu \nabla^2 \mathbf{u}, \qquad (3) $$

$$ \frac{d}{dt} \left( \frac{p}{\rho^\gamma} \right) = 0, \qquad (4) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \frac{\eta c^2}{4\pi} \nabla^2 \mathbf{B}, \qquad (5) $$

where $\gamma$ is the adiabatic index (usually 5/3), $\nu$ is the dynamic viscosity, $\eta$ is the resistivity and $\eta c^2/4\pi$ collectively is known as the magnetic diffusivity. Here Eq. (2) is the continuity equation representing conservation of mass density, Eq. (3) is the momentum conservation equation, Eq. (4) is the simple adiabatic gas equation representing the conservation of energy and Eq. (5) is the induction equation, where the first term on the right is the advection term and the second one represents diffusion. Note, as a result of MHD approximations, the displacement current term is omitted in the induction equation. In an ideal situation, there are no dissipative processes and therefore $\nu$ and $\eta = 0$ gives ideal MHD equations. Studies have shown that ideal MHD description is a very good approximation to study dynamical properties of strongly magnetized plasmas.
## 2.2 Regional MHD Simulations

Several large-scale MHD approaches do not model the entire magnetosphere but only sections of it, such as certain magnetopause regions (dealt with elsewhere in this volume) or the magnetotail. Here we focus particularly on the magnetotail. The basic numerical approach used in regional MHD simulations is essentially similar to that used in (some) global simulations. It is typically based on explicit finite difference methods to solve the MHD equations. Minor differences might exist in adding resistive terms, which are usually necessary in local MHD to initiate reconnection. On the other hand, many global models use finite volume approaches in part due to the complex geometries involved (e.g., Powell et al. 1999).

The main difference, however, consists of the setup or initialization. Whereas global simulations typically involve of period of interaction with the solar wind to create a realistic magnetotail, local MHD simulations generally start from some equilibrium or near-equilibrium that models the stretched magnetotail (e.g., Schindler 1972). This approach provides more flexibility in treating different scenarios, for instance, varying the tail flaring between $y$ and $z$ (Birn and Hesse 2000) or including a local $B_z$ hump (Merkin and Sitnov 2016; Birn et al. 2018). This flexibility has also proven useful in PIC simulations that go beyond the commonly used initial 1D Harris sheet, most notably in addressing the holy grail of reconnection onset (e.g., Liu et al. 2014; see also Liu et al. 2025a).

On the other side, interactions with the ionosphere or the solar wind are incorporated only in some ad hoc fashion, if at all. Regional magnetotail MHD simulations therefore have been most successful in treating dynamic tail phenomena on relatively short time scales that are typically substorm related. The successes include

1. The demonstration that X-line formation and plasmoid ejection can be part of a 2D or 3D tearing-type instability of the tail (e.g., Birn and Hones 1981).

2. The demonstration that the build-up of the substorm current wedge (SCW), involving dipolarization and Region-1-type field-aligned currents (McPherron et al. 1973), can be due to the braking and azimuthal diversion of earthward flow from a near-tail reconnection site (Birn and Hesse 1991; Scholer and Otto 1991). This basic picture has been modified more recently, most notably by the addition of a Region-2 current system connecting to the ionosphere at lower latitude (Birn and Hesse 2014; Kepko et al. 2015), in agreement with observations (Sergeev et al. 2014). While the buildup of the SCW in the simulations is based on the shear and vorticity of the earthward flow, the persistence of the currents relies on the changes of the magnetic flux and pressure patterns brought about by the severance of a plasmoid and the resulting loss of entropy and redistribution of the pressure. It is noteworthy that these features can be, and have been, found also in global simulations. In the regional simulations, however, they arise as consequences of an instability without involvement of external driving or feedback from the ionosphere. This would be harder to extract from the global simulations.

3. Regional MHD simulations have also been used to address the evolution prior to the onset of reconnection in the tail, demonstrating, specifically, the formation of a thin concentrated current sheet embedded in the near-tail plasma sheet. These approaches have included interaction with the solar wind in two complementary, ad-hoc ways. In one approach magnetic flux is added to the tail lobes (Birn and Schindler 2002), the other is based on low-latitude magnetic flux reduction from convection around the Earth toward the dayside (Hsieh and Otto 2014, 2015). Both mechanisms are expected from solar wind interaction. For more details, see the review by Sitnov et al. (2019).

4. Recently, regional MHD simulations have also demonstrated that the magnetotail may become unstable even under ideal 2D MHD constraints, when it includes a region of
inverse (i.e. tailward) gradient of the normal magnetic field $B_z$ (denoted ‘$B_z$ hump’ instability; (Merkin and Sitnov 2016; Birn et al. 2018)).

## 2.3 Global MHD

Global MHD models representing the (outer) magnetosphere of Earth typically extend around 100–200 $R_E$ in the flank and tail directions, and around 30 $R_E$ (beyond the bow shock) towards the Sun, where the solar wind is coming from. When, occasionally, the solar wind becomes sub-Alfvénic, the bow shock disappears and Alfvén wings form. In this case the upstream boundary has to be moved much further to minimize the boundary effects. It is computationally very demanding to obtain an accurate solution in such a large domain while resolving various structures such as the current sheets at the dayside magnetopause and in the tail. There are various approaches to overcome this difficulty, including the block-adaptive grid of the BATSRUS code (Powell et al. 1999; Tóth et al. 2012), the stretched Cartesian grid of OpenGGCM (Raeder et al. 1997) or the stretched spherical grid of the Lyon-Fedder-Mobarry (LFM) code (Lyon et al. 2004).

A general issue associated with all MHD simulations is preserving $\nabla \cdot \mathbf{B} = 0$ (e.g., Tóth et al. 2012). This issue also creates issues associated with the upstream boundary condition. Typically, we have observations at a single point near L1, a Lagrange point in between the Sun and the Earth, and assume that the solar wind and interplanetary magnetic field (IMF) have no variation in the transverse direction. If $B_x$ varies, and it certainly does, these assumptions lead to a finite $\nabla \cdot \mathbf{B} = \partial B_x / \partial x$ propagating into the domain. There are various approaches to handle this situation. One is to ignore the problem and propagate the finite $\nabla \cdot \mathbf{B}$ with the flow using some variation of the 8-wave scheme (Powell 1994). Another common approach is to set $B_x$ to a constant value, for example 0. Finally, one can relax the condition that the transverse gradients in the y and z directions are zero, and guess those gradients from the temporal evolution of $B_x$. Unfortunately this is an underspecified problem, so there is no unique solution. A typical approach is to smooth $B_x$ in time, and apply some minimum variance constraint. While theoretically nice, in practice the minimum variance approach does not work great. The likely reason is that the magnetic field is turbulent, so local changes in $B_x$ are not representative of the large scale tilt of the propagation plain.

The inner boundary conditions are usually applied at a sphere of radius 1.5 $R_E$ to 3 $R_E$ surrounding the Earth. One can use a semi-empirical electrodynamic solver (Ridley et al. 2004), or a fully empirical model (Weimer 1996, 2001) to calculate the $\mathbf{E} \times \mathbf{B}$ drift velocity at the inner boundary. Another important use of an electrodynamic solver is that it can also provide the $\mathbf{E} \times \mathbf{B}$ drift to an inner magnetosphere model (Wolf et al. 1982; Toffoletto et al. 2003; Buzulukova et al. 2010; Liemohn et al. 2001; Jordanova et al. 1994; Zaharia et al. 2006), which can calculate realistic ring current and associated pressure (and density) in the closed field line region during geomagnetic storms (Liemohn et al. 2018). The global MHD model can then relax its pressure (and density) towards the values supplied by the inner magnetosphere model. In return, the global MHD model can supply the plasma boundary conditions for the inner magnetosphere model at the edge of the closed field region as well as the magnetic field configuration (De Zeeuw et al. 2004; Meng et al. 2013).

Global models can properly represent the overall dynamics of the interaction of the solar wind with the magnetosphere, including the formation of the bow shock and the magnetopause, as well as the main current sheet in the magnetotail. Magnetic reconnection will happen on the dayside magnetopause and in the magnetotail in agreement with the theory of the Dungey cycle (Dungey 1961). The magnetic reconnection in the MHD simulation is not represented by the actual kinetic physics but it is approximated by numerical diffusion or
artificial resistivity. Despite these caveats, global MHD models generate reconnection sites where the magnetic field changes sign and the reconnection rate is approximately correct (see Appendix A of Wang et al. (2022a)).

When the grid is relatively coarse, the numerical diffusion will easily adjust to reconnect the incoming magnetic flux carried by the solar wind. For a constant solar wind and IMF driving the simulation will settle to a steady state solution. Using fine grids in combination with low dissipation numerical methods can lead to a more dynamic reconnection process in the model. On the dayside, simulations can produce Flux Transfer Events (Raeder 2006) and in the tail flux ropes can be produced even by ideal MHD simulations. If the flux ropes in the tail are triggered by sign changes of the IMF $B_z$, the MHD simulations can match observations very well. A more challenging problem is reproducing substorms and sawtooth events. MHD models cannot do this well (Haiducek et al. 2020), and one needs to add either ionospheric outflow to regulate the reconnection rate (Brambles et al. 2013; Zhang et al. 2020), or kinetic reconnection physics (Wang et al. 2022a) to produce the typical spatial and temporal scales of sawtooth oscillations.

## 2.4 Test Particles in MHD Simulations

Test particle approaches consist of tracing charged particle orbits in electromagnetic fields that are either prescribed in some plausible fashion or obtained from a simulation that typically does not contain individual particle information, most commonly based on MHD. The approach bridges the gap between large-scale MHD and small-scale particle simulations. In contrast to the latter, it can treat together realistic 3D space and large evolution time scales and realistic electron mass. However, it is not self-consistent and relies on whether the MHD model or the postulated **E**, **B** fields capture the main physics. But that may also be considered an advantage as it permits studying the effects of large-scale fields in isolation.

In the magnetospheric context ions are usually treated by integration of the full orbit

$$ \frac{Du}{Dt} = \frac{e}{m} \left( \mathbf{E} + \frac{1}{\gamma c} \mathbf{u} \times \mathbf{B} \right). \qquad (6) $$

Here $\gamma$ is the relativistic factor, which may be more relevant for electrons, $\mathbf{u} = \gamma \mathbf{w}$, where $\mathbf{w}$ is the particle velocity, $c$ is the speed of light and $D/Dt = \partial/\partial t + \mathbf{w} \cdot \nabla$ denotes the derivative along the full orbit. This equation is typically evolved using a high order Runge-Kutta (Press et al. 1992) or the Boris method (Birdsall and Langdon 1991). Full integration of Equation (6) over extended orbits is not practical for electrons, as it is more time consuming and might accumulate too large errors. Also, the adiabatic drift approximation, based on conservation of the magnetic moment $\mu$, is valid over larger areas in the magnetosphere and can be adequate for identifying typical acceleration mechanisms (e.g., Delcourt and Sauvaud 1994; Li et al. 1998; Zaharia et al. 2000; Gabrielse et al. 2012).

However, when the full history of electron orbits is considered; this may involve encounters of the reconnection site and low magnetic field, or high curvature regions, where the conservation of $\mu$ breaks down, and full orbit integration is required. Consequently, several codes have been developed that involve a transition between full orbits, integrated by Eq. (6), and drift orbits (e.g., Birn et al. 2004; Schriver et al. 2005; Ashour-Abdalla et al. 2011; Sorathia et al. 2017). The drift is described by the guiding center drift velocity (e.g., Birn et al. 2004)

$$ \mathbf{v}_d = \mathbf{v}_E - \frac{\mu c}{\gamma e} \frac{\mathbf{B} \times \nabla B}{B^2} - \frac{\gamma m_e c v_\parallel}{e} \frac{\mathbf{B}}{B^2} \times \frac{d\mathbf{b}}{dt} - \frac{m_e c}{e} \frac{\mathbf{B}}{B^2} \times \frac{d(\gamma \mathbf{v}_E)}{dt} \qquad (7) $$
where $\mu$ is the (relativistic) magnetic moment, $\mathbf{v}_E = \mathbf{E} \times \mathbf{B}/B^2$ and $\mathbf{b} = \mathbf{B}/B$. In addition, the field-aligned velocity is advanced by

$$ \frac{du_{\parallel}}{dt} = -\frac{e}{m_e} E_{\parallel} - \frac{\mu}{\gamma m_e} \frac{\partial B}{\partial s} - (\mathbf{u}_E + \mathbf{u}_{\nabla B}) \cdot \frac{d\mathbf{b}}{dt} \quad (8) $$

where now $\mathbf{u} = \gamma \mathbf{v}$ and $\mathbf{v} = \mathbf{v}_E + \mathbf{v}_d + \mathbf{v}_{\parallel}$ describes the guiding center velocity, $\mathbf{v}_{\nabla B}$ is the grad B drift, given by the second term on the right side of Eq. (7), and $d/dt$ is the derivative along the guiding center path. The transition between full orbit and drift orbit is typically determined from an adiabaticity criterion that is based on the ratio between the local field line curvature radius and the gyro radius based on the local magnetic field strength (e.g., Buchner and Zelenyi 1989). On the switch from drift to full orbit a phase has to be generated, which is typically chosen randomly. Although this can alter individual orbits and make them not reversible, it was found to have no significant effect on general conclusions about sources and properties of distributions (e.g., Birn et al. 2004).

Two different techniques are used, tracing particle motion either forward or backward in time. Forward tracing requires larger numbers of particles, sometimes comparable to those in PIC simulations, to obtain sufficient numbers at the points of interest (e.g., Scholer and Jamitzky 1987; Sachsenweger et al. 1989; Peroomian and El-Alaoui 2008; Ukhorskiy et al. 2017, 2018). However, since the particles are not interacting this approach is even more suitable for parallel processing than full particle simulations. In principle, this approach can also include wave scattering and collisions (albeit in an ad-hoc non self-consistent manner) to add to the simple collisionless advance.

Backward tracing is generally based on Liouville's theorem of the conservation of phase space density $f$ to map $f$ from source locations to the final location of interest (e.g., Curran and Goertz 1989; Birn and Hesse 1994; Birn et al. 2004). It requires fewer orbits to identify properties at selected final locations, but relies on the validity of Liouville's theorem, i.e. the absence of collisions. Backward tracing permits an easier identification of different sources contributing to the final population. Thus, sometimes a combination of both techniques is employed (e.g., Ashour-Abdalla et al. 2011).

Further complications are related to the use of MHD simulation results. Since the fields are given only on a finite grid, they have to be interpolated in space and time.

1. The advance of the drift equations (7) and (8) requires a third order spatial interpolation in B for continuous transition between grid cells, which could lead to spurious maxima or minima. This can be avoided, however, by employing a monotonicity algorithm (Hyman 1983; Birn et al. 2004).

2. Simple interpolation of the electric field could also yield spurious parallel components. This can be avoided, however, by various techniques, for instance, by interpolating $E_{\parallel}$ and $E_{\perp}$ separately (e.g., Birn et al. 2022).

# 3 Hall MHD

## 3.1 Introduction

The ideal-MHD model, as discussed in Sect. 2.1, is well-suited for magnetized plasmas when the dynamics is slow compared to the gyration time of charged particles around the magnetic fields and the length scales over which quantities vary is much larger than the
gyroradius of the charged particles. However, going back many years in the study of neutral fluids, fluid models can lead to incorrect and paradoxical results at boundaries layers [e.g., d'Alembert's paradox (Sect. 4.7 of Choudhuri 1998)]. In a magnetized plasma, these problems can occur where plasmas of two different origins abut against each other (such as at Earth's magnetopause), at shocks and discontinuities such as Earth's bow shock, and at localized regions where the magnetic field goes to zero, such as in the solar corona near sunspots.

Magnetic reconnection, in particular, occurs at a boundary layer at a region where at least two components of the magnetic field go to zero, so it is a key example of a physical process that cannot be faithfully modeled by ideal-MHD. Often in numerical simulations, ideal-MHD is used anyway, with numerical dissipation allowing reconnection to occur with the hope that it mimics the actual process. Another approach employs resistivity to model the effect of collisions; this is a useful approach in systems for which collisions are dynamically relevant, but many settings where reconnection occurs – especially in space and the solar corona – are weakly collisional or effectively collisionless (Priest and Forbes 2000; Cassak and Shay 2012). There are examples where either approach can be good enough for the questions being asked. For other questions that rely on a faithful representation of the physics in the regions where ideal- and resistive-MHD break down, a new model is necessary. In this section, we discuss a number of approaches within the fluid description that are used to go beyond ideal- and resistive-MHD simulations. In later sections, simulation techniques using the kinetic theory of gases are discussed. There are previous review papers discussing Hall-MHD and numerical approaches (Vasyliunas 1975; Huba 1995, 2003; Gómez 2006).

## 3.2 The Hall-MHD Model

The equations of Hall-MHD are similar to those of MHD with one key difference. In resistive-MHD, Ohm's law is given (in cgs units) by

$$E + \frac{u \times B}{c} = \eta J, \qquad (9)$$

where **E** and **B** are the electric and magnetic fields, **u** is the single fluid bulk flow velocity, $\eta$ is the resistivity, and **J** = $(c/4\pi)\nabla \times **B**$ is the current density; in ideal-MHD, $\eta$ is set to zero. To go beyond this model, we revisit where Ohm's law comes from.

The equation of motion of an electron fluid (i.e., Newton's 2nd law) in a fully ionized plasma (in cgs units) is (Braginskii 1965)

$$m_e \frac{d\mathbf{u}_e}{dt} = -e \left( \mathbf{E} + \frac{\mathbf{u}_e \times \mathbf{B}}{c} \right) - \frac{1}{n_e} \nabla \cdot \mathbf{p}_e + \mathbf{R}_e, \qquad (10)$$

where $m_e$ is the electron mass, $\mathbf{u}_e$ is the electron bulk flow velocity, $-e$ is the electron charge, $n_e$ is the electron density, $\mathbf{p}_e$ is the electron pressure which we write more generally as a tensor for now, and $\mathbf{R}_e$ represents the rate of change of momentum resulting from collisions between electrons and other electrons or other charged or neutral particles in the plasma. There are rigorous ways to determine the role of collisions (e.g., Braginskii 1965) that we do not employ here. Instead, we use the often used simpler approach that assumes $\mathbf{R}_e = m_e v_{ei} (\mathbf{u}_i - \mathbf{u}_e)$, where $v_{ei}$ is the electron-ion collision frequency, and $\mathbf{u}_i$ is the ion bulk flow velocity, and for simplicity we assume the plasma has only electrons and ions (it is fully ionized).
To recover the resistive-MHD Ohm's law from this equation, first the “electron inertia term” $m_e d\mathbf{u}_e/dt$ and the “electron pressure gradient term” $-(\nabla \cdot \mathbf{p}_e)/n_e$ are ignored for reasons we return to in Sect. 3.6. Second, the single fluid bulk flow velocity used in MHD is $\mathbf{u} = (m_i n_i \mathbf{u}_i + m_e n_e \mathbf{u}_e)/(m_i n_i + m_e n_e)$ and the current density $\mathbf{J} = n_i q_i \mathbf{u}_i - n_e q_e \mathbf{u}_e \simeq n_e q_e (\mathbf{u}_i - \mathbf{u}_e)$, where the latter form uses the assumption of quasi-neutrality $n_i q_i - n_e q_e \simeq 0$ with $q_i$ and $q_e$ being the charge of ions and electrons respectively. Using these expressions to write $\mathbf{u}_e$ in terms of $\mathbf{u}$ and $\mathbf{J}$ gives $\mathbf{u}_e = \mathbf{u} - (\mathbf{J}/ne)[m_i n_i/(m_i n_i + m_e n_e)] \simeq \mathbf{u} - (\mathbf{J}/ne)$, where in the latter form, we use the approximation that $m_i \gg m_e$, since it is at least 1836 in an electron-ion plasma. Using these approximations in Eq. (10) and dividing by $e$ gives

$$ \mathbf{E} + \frac{\mathbf{u} \times \mathbf{B}}{c} = \frac{\mathbf{J} \times \mathbf{B}}{n_e e c} + \eta \mathbf{J}, \quad (11) $$

where the resistivity $\eta$ is defined as $m_e v_{ei} / n_e e^2$. If one additionally ignores the $\mathbf{J} \times \mathbf{B} / n_e e c$ term, what remains is the resistive-MHD Ohm's law in Eq. (9). If instead, one ignores the resistive term $\eta \mathbf{J}$, the result is

$$ \mathbf{E} + \frac{\mathbf{u} \times \mathbf{B}}{c} = \frac{\mathbf{J} \times \mathbf{B}}{n_e e c}. \quad (12) $$

The term on the right is called the “Hall electric field” $\mathbf{E}_H$ (or simply the “Hall term”), and Eq. (12) is called the “Hall-MHD Ohm’s law”. Simply coupling this equation to the rest of the ideal-MHD equations gives the Hall-MHD model:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \quad (13) $$

$$ \rho \left[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right] = -\nabla p + \frac{\mathbf{J} \times \mathbf{B}}{c}, \quad (14) $$

$$ \frac{\partial p}{\partial t} + (\mathbf{u} \cdot \nabla) p = -\gamma p (\nabla \cdot \mathbf{u}) \quad (15) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}, \quad (16) $$

$$ \mathbf{E} + \frac{\mathbf{u} \times \mathbf{B}}{c} = \frac{\mathbf{J} \times \mathbf{B}}{n_e e c} \quad (17) $$

$$ \nabla \times \mathbf{B} = \frac{4\pi \mathbf{J}}{c} \quad (18) $$

with the auxiliary equation $\nabla \cdot \mathbf{B} = 0$, and where $\gamma$ is the single fluid ratio of specific heats, typically taken to be 5/3. Note $n_e$ in Eq. (17) is related to the MHD mass density via $\rho = m_i n_i + m_e n_e$, so $n_e = \rho/(m_i e/q_i + m_e) \simeq Z\rho/m_i$, where we again assume $m_i \gg m_e$, we use quasi-neutrality to write $q_i n_i \simeq e n_e$, and define $Z = q_i/e$ as the degree of ionization. With these assumptions, Eqs. (13) - (18) form a closed set of equations, and therefore can be used to model physical systems. Technically, these equations actually give the “ideal Hall-MHD model” since resistivity is not retained here.

It is important to note a confusing aspect of these equations. The Hall term is proportional to $\mathbf{J} \times \mathbf{B}$, and Eq. (14) also contains a term including $\mathbf{J} \times \mathbf{B}$. It is tempting to draw relations between the two terms because of this outward similarity, but this should not be done. The two terms have different dimensions: $\mathbf{J} \times \mathbf{B}/c$ is a force density and $\mathbf{J} \times \mathbf{B}/n_e e c$ is an electric field. They have completely different manifestations and impacts on the physics, and therefore actually are not related despite their similar forms.
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Physics of the Hall effect in reversing magnetic fields. A reversing anti-parallel magnetic field in the $\pm\hat{x}$ direction is drawn using black arrows, with the neutral line as the dashed line. A uniform electric field in the $\hat{z}$ direction is drawn in green. The smaller blue trajectory is that of an electron $E \times B$ drifting towards the neutral line with bulk velocity $u_e$. An ion $E \times B$ drifting far from the neutral line has bulk ion velocity $u_i$ identical to $u_e$, so the electric field in this region is given by $E = -u \times B/c$, where $u$ is the single fluid (MHD) bulk flow velocity. Within an ion gyroradius $r_{Li}$ of the neutral line, the electron continues with bulk flow $u_e$, while the ion demagnetizes (in the red trajectory) and $u_i$ becomes small. In this region, there is a non-zero current density $J$ (orange), and the electric field in this region is predominantly given by the Hall electric field $E = J \times B/n_e e c$</figcaption></figure>

## 3.3 Hall-MHD Physics

The only difference between the ideal-MHD and Hall-MHD models is the Hall electric field, and here we investigate the physics introduced by this term. Let $\tilde{L}$ be a characteristic length scale over which the plasma properties vary, and let $\tilde{r}_{Li}$ be the characteristic (Larmor) radius of the ions as they gyrate around a magnetic field of characteristic strength $\tilde{B}$. The Hall electric field is small and can be neglected if $\tilde{L} \gg \tilde{r}_{Li}$, and doing so brings us back to ideal-MHD. It is important to retain the Hall electric field, and therefore use Hall-MHD, when $\tilde{L} \lesssim \tilde{r}_{Li}$.

To see why structure below the ion gyroscale gives rise to the Hall effect, consider a magnetic field **B** that reverses direction over a scale $\tilde{r}_{Li}$ or less, but on larger scales than the characteristic electron gyroradius $\tilde{r}_{Le}$. For specificity, consider a magnetic field pointing in the $\pm\hat{x}$ direction, as sketched as the black arrows in Fig. 1. The neutral line, where the magnetic field strength vanishes, is the dashed black line. Suppose also there is a uniform electric field **E** pointing everywhere in the $\hat{z}$ direction, as sketched as the green arrows. (We use a reversed magnetic field and uniform electric field for illustrative purposes due to its relation to the reconnection process, but the Hall effect is important for any magnetic field configuration that varies on $\tilde{r}_{Li}$ scales.)

At distances from the magnetic field reversal that greatly exceed $\tilde{r}_{Li}$, ions and electrons undergo the $E \times B$ drift that gives rise to a bulk flow towards the neutral line; the bulk flow velocity of ions and electrons, $u_i$ and $u_e$, respectively, are identical. The $E \times B$ drift of the electrons above the neutral line is sketched as the blue curve. In the region farther from the neutral line than $\tilde{r}_{Li}$, the electric field is given by $E = -u \times B/c$, where $u = u_i = u_e$ is the MHD bulk flow velocity.

As the ions reach a distance from the neutral line that is equal to its gyroradius, the ions cross the neutral line and are immersed in a magnetic field pointing in the opposite direction.
Their gyromotion changes direction, and they make figure 8 orbits, sketched as the red curve. (They also accelerate in the $\hat{z}$ direction due to the electric field, but this is omitted from the sketch.) Consequently, their bulk velocity becomes small, $\mathbf{u}_i \simeq 0$.

Since we assumed $\tilde{L} > \tilde{r}_{Le}$, the electrons have a smaller gyroradius and do not see the magnetic field reversal, so they continue to undergo the $E \times B$ drift towards the neutral line. The key is that the ions and electrons are undergoing different dynamics between distances from the neutral line of $\tilde{r}_{Li}$ and $\tilde{r}_{Le}$!

In this region, the difference in the bulk motion between ions and electrons implies there is a net current density $\mathbf{J}$, sketched as the orange arrow, and called the Hall current. The current density is perpendicular to the magnetic field $\mathbf{B}$. Then, between $\tilde{r}_{Li}$ and $\tilde{r}_{Le}$ from the neutral line, the electric field is given by the Hall electric field $\mathbf{E}_H = \mathbf{J} \times \mathbf{B}/n_e e c$. This exemplifies why the Hall electric field is important between ion and electron gyroscales.

This situation in a plasma is analogous to the Hall effect in condensed matter physics, where it was originally discovered by Edwin Hall (a graduate student) in 1879 (Hall et al. 1879). It has extensive applications in that field of study. The derivation generalizing shear Alfvén waves to include the Hall effect happened as early as 1954 by Jim Dungey (Dungey 1954), the same person who first understood magnetic reconnection and gave the process its name.

## 3.4 Linear Waves in Ideal Hall-MHD

It would take us too far afield to elucidate how the Hall electric field modifies all the physics of ideal-MHD. Rather, we highlight one important example – linear waves. In ideal-MHD, there are three propagating linear waves available to a uniform plasma: the shear Alfvén wave and the fast and slow magnetosonic waves. The Hall electric field introduces two wave modes that become important between ion and electron gyroscales – the whistler wave and the kinetic Alfvén wave.

### 3.4.1 The Whistler Wave

First, we consider transverse waves propagating along a uniform background magnetic field of strength $B_0$, which without loss of generality we take to be in the $y$ direction, in a plasma of equilibrium mass density $\rho_0 \simeq m_i n_i$ and with no equilibrium bulk flow or current $\mathbf{u}_0 = 0$ and $\mathbf{J}_0 = 0$. The wave vector $\mathbf{k}$ is also in the $y$ direction. Linearizing the Hall-MHD equations (13)-(18) about the given equilibrium with $\mathbf{k}$ parallel to $\mathbf{B}_0$ and solving for the dispersion relation gives

$$ \omega^2 = k^2 v_{A0}^2 \left( 1 + \frac{k^2 d_{i0}^2}{2} + \sqrt{k^2 d_{i0}^2 + \frac{k^4 d_{i0}^4}{4}} \right), \quad (19) $$

where $\omega$ is the wave frequency, $v_{A0} = B_0/(4\pi\rho_0)^{1/2}$ is the Alfvén speed, and $d_{i0} = (m_i c^2/(4\pi n_0 q_i^2))^{1/2}$ is the ion inertial length.

In the limit of $kd_{i0} \to 0$, this dispersion relation reduces to $\omega^2 \to k^2 v_{A0}^2$, which is simply the shear Alfvén wave from ideal-MHD. In the other extreme, consider the limit $kd_{i0} \to \infty$. The 1 in Eq. (19) becomes negligible compared to $k^2 d_{i0}^2/2$ outside the square root, and the $k^2 d_{i0}^2$ is negligible compared to $k^4 d_{i0}^4/4$ inside the square root, so to low order the dispersion relation becomes

$$ \omega^2 \to k^4 v_{A0}^2 d_{i0}^2 \quad (\text{as } kd_{i0} \to \infty). \quad (20) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Sketch of how the Hall effect impacts linear waves in MHD. The wave structure for (a) shear Alfvén waves, (b) parallel propagating whistler waves, and (c) kinetic Alfvén waves. The equilibrium magnetic field $B_{y0}$ is the dashed black arrow. The perturbed magnetic field $B_1$ are the black arrows. The perturbed ion bulk flow $\mathbf{u}_i$ are the red arrows, which occur in (a) because the wavelength is much larger than the ion gyroradius. In (b) and (c), the wavelength is at or below ion gyroscales, so the bulk flow is due to electron motion in the blue arrows. In (c), there is a large out-of-plane equilibrium magnetic field $B_{z0}$, so the magnetic perturbation $B_{z1}$ changes the magnetic pressure which to first order and requires a change to the gas pressure $p$.</figcaption></figure>

This is the dispersion relation for the so-called “whistler wave.”

Physically, the whistler wave is the sub-ion gyroscale counterpart of the shear Alfvén wave. The plasma properties in a shear Alfvén wave are sketched in Fig. 2(a). The equilibrium magnetic field $B_{y0}$ is the dashed black arrow. A magnetic perturbation transverse to the equilibrium $B_{x1}$ is the large black arrow. Since the size of the wave is far larger than the ion gyroradius, the frozen-in ions and electrons feel a restoring force analogous to plucking a guitar string, generating bulk flows in the $x$ direction (shown for ions as the red arrows). The shear Alfvén wave is a linearly polarized wave.

The whistler wave, sketched in panel (b), occurs when the magnetic field varies between ion and electron gyroscales. Then, the ions are not frozen-in and do not respond to the plucked magnetic field line, but the electrons are frozen in and feel a restoring force. The $x$ component of the perturbed electron flow $u_{ex1}$ is sketched as the blue arrows. Since the electrons have a bulk velocity but the ions do not, this means there is a net current density, which induces an oscillating magnetic field out of the plane, labeled as $B_{z1}$. Similarly, the varying $B_{x1}$ requires a $u_{ez1}$ to sustain the current. This turns the wave into a circularly polarized wave. The polarization is right handed, i.e., a receiver sees the magnetic field rotate in a counterclockwise direction.

There is a crucial difference between the shear Alfvén wave and the whistler wave. The shear Alfvén wave has a phase speed $\omega/k = v_{A0}$, which is a constant independent of $\omega$ and $k$, so a shear Alfvén wave travels at the same speed regardless of its wavelength. It is an example of a non-dispersive wave that retains its waveform as it propagates (like a light wave in vacuum or a sound wave in a neutral fluid). In contrast, the phase speed for the whistler wave, from Eq. (20), is $\omega/k = kv_{A0}d_{i0}$. Thus, the phase speed is faster for shorter wavelength waves. This is an example of a dispersive wave, since a wave packet does not retain its shape. The dispersive nature is what gives the wave its name; it makes a characteristic “whistle” from high frequencies descending to low frequencies when it is detected at a location away from where it was generated. For example, when a lightning strike in one hemisphere excites a whistler wave on Earth’s magnetic field, the high frequencies arrive before the low frequencies in the other hemisphere.

This analysis reveals the critical scale at which the Hall effect becomes important for transverse bending of the magnetic field. The three terms in Eq. (19) proportional to $d_{i0}$ are absent if the Hall term is left out of the governing equations (i.e., in ideal-MHD). By
comparing the first term to the other terms, we see they become important when $kd_{i0}$ is on the order of 1. This means that the Hall term becomes important at length scales below $d_{i0}$.

The ion inertial scale is related to a form of the ion gyroradius. To see this, note that an ion moving at the Alfvén speed $v_{A0}$ that gyrates around a magnetic field of strength $B_0$ has a gyroradius $r_{Li}$ given by

$$r_{Li} = \frac{v_{A0}}{\Omega_{ci0}} = d_{i0}, \qquad (21)$$

where $\Omega_{ci0} = q_i B_0 / m_i c$ is the ion gyrofrequency. Thus, the scale at which the Hall effect becomes important for transverse perturbations is when gradient scales are comparable to the ion inertial scale. Using characteristic scales for Earth's dayside magnetosheath, $d_{i0} \simeq 70$ km; for the solar corona, $d_{i0} \simeq 2$ m. Thus, we have a feel for length scales over which we need to use Hall-MHD instead of ideal-MHD in two important space applications.

### 3.4.2 The Kinetic Alfvén Wave

Now, we consider the kinetic Alfvén wave. These waves are almost completely longitudinal, but not perfectly longitudinal. To allow for a strong analogy to the whistler wave, consider a uniform magnetic field that has a very large z component $B_{z0} \gg 0$ and very small y component $B_{y0}$, as sketched in Fig. 2(c). As in panels (a) and (b), the wave propagates in the y direction. If $B_{z0}$ were zero, it would be the whistler wave. Physically, if the wavelength is between ion and electron gyroscales, this perturbation reacts similar to a whistler in that it sets up an electron flow that bends the magnetic field out of the plane. The perturbed magnetic field therefore has a component in the same or opposite direction as $B_{z0}$. Where $B_{z0}$ and the perturbed magnetic field are parallel, the magnetic pressure $(B_{z0} + B_{z1})^2 / 8\pi$ is greater than the initial pressure $B_{z0}^2 / 8\pi$ to first order in the perturbed field. (Note, for the whistler wave, the pressure difference $B_{z1}^2 / 8\pi$ is second order in the perturbed magnetic field, so it is negligible.) Similarly, where $B_{z1}$ opposes $B_{z0}$, the magnetic pressure decreases to first order. Now, if the plasma is overall low $\beta$ (since we took $B_{z0}$ to be large), this magnetic pressure imbalance cannot be maintained, so the gas pressure has to change in order to balance pressure. Electrons move along the magnetic field and ions move across the magnetic field to move from the high magnetic pressure region to the low magnetic pressure region, setting up a high and low gas pressure region as denoted in panel (c). This describes the physics of the kinetic Alfvén wave.

To find the dispersion relation of the kinetic Alfvén wave, we back up to the full dispersion relation of any wave in Hall-MHD. Linearizing Eqs. (13) - (18) around a uniform magnetic field, density, and pressure with an arbitrary linear perturbation and solving gives the following dispersion relation (Rogers et al. 2001):

$$\omega^6 - (c_{ms}^2 + v_{Ak}^2 + k^2 v_{Ak}^2 d_{i0}^2)k^2\omega^4 + [c_{ms}^2 + c_s^2(1 + k^2 d_{i0}^2)]k^4 v_{Ak}^2\omega^2 - k^6 v_{Ak}^4 c_s^2 = 0, \qquad (22)$$

where $v_A^2 = B_0^2 / (4\pi\rho_0)$ is the total Alfvén speed, $c_s^2 = \gamma\rho_0 / \rho_0$ is the sound speed, $c_{ms}^2 = c_s^2 + v_A^2$ is the total fast magnetosonic speed, and $v_{Ak}^2 = B_{y0}^2 / (4\pi\rho_0)$ is the Alfvén speed based only on $B_{y0}$. In the $k^2 d_{i0} \to 0$ long wavelength limit, this equation reduces to the dispersion relation for ideal-MHD waves. In the limit of large $kd_{i0}$, one solution is a high frequency solution which is approximately given by balancing the first two terms in Eq. (22), which gives
$\omega^2 \simeq k^4 v_{Ak}^4 d_{i0}^2$, the whistler wave dispersion relation in Eq. (20). There is also a medium frequency solution which arises from balancing the middle two terms in Eq. (22). This ratio in general is

$$ \omega^2 \approx \frac{c_{ms}^2 + c_s^2(1 + k^2 d_{i0}^2)}{c_{ms}^2 + v_{Ak}^2 + k^2 v_{Ak}^2 d_{i0}^2} k^2 v_{Ak}^2. \quad (23) $$

In the limit in which $c_{ms}^2 \ll c_s^2 k^2 d_{i0}^2$ and $c_{ms}^2 \gg v_{Ak}^2 + k^2 v_{Ak}^2 d_{i0}^2$, the resulting dispersion relation is $\omega^2 \simeq (c_s^2/c_{ms}^2)k^4 d_{i0}^2 v_{Ak}^2$, which in the $v_A^2 \gg c_s^2$ (low $\beta$) limit gives $\omega^2 \simeq (c_s^2/v_A^2)k^4 d_{i0}^2 v_{Ak}^2$. Since $k^2 v_{Ak}^2 = k_\parallel^2 v_A^2$, this becomes

$$ \omega^2 \approx k_\parallel^2 k^2 v_A^2 \rho_s^2, \quad (24) $$

where $\rho_s^2 = c_s^2/\Omega_{ci}^2$ is the ion Larmor radius based on the sound speed. This is the dispersion relation for the kinetic Alfvén wave. As with the whistler wave, the kinetic Alfvén wave is dispersive with $\omega/k \propto k_\parallel$, so it gets faster for smaller wavelengths. The length scale at which the Hall term becomes important for the kinetic Alfvén wave is $\rho_s$ (as opposed to $d_{i0}$ for the whistler wave).

## 3.5 A Numerical Algorithm for Hall-MHD

Including the Hall electric field has a significant impact on numerical simulations relative to MHD simulations. One way to think of why this is the case is that the waves in ideal-MHD are non-dispersive, so waves at any scale from the large scale down to the computational grid scale travel at the same speed. In Hall-MHD, as discussed in the previous section, both whistler and kinetic Alfvén waves are dispersive, so waves at the grid scale (which needs to be sub-ion gyroscale to capture the Hall electric field) are considerably faster than waves at the large scale. This makes the Hall-MHD equations “stiff” – one must use a much smaller time step in Hall-MHD than in ideal-MHD, leading to a significant increase in the run time and expense of the simulation.

There are numerous algorithms that can be used to numerically evolve the equations in time. We provide one in detail, and mention references to other algorithms that have been used. We highlight the F3D code (Shay et al. 2004), which has been used for many years to study magnetic reconnection. First, the evolution equations are written in conservative form as

$$ \frac{\partial n}{\partial t} + \nabla \cdot \mathbf{J}_i = 0, \quad (25) $$

$$ \frac{\partial \mathbf{J}_i}{\partial t} + \nabla \cdot \left( \frac{\mathbf{J}_i \mathbf{J}_i}{n} + \frac{p \mathbf{I}}{m_i} + \frac{\mathbf{B} \mathbf{B}}{4\pi m_i} - \frac{B^2 \mathbf{I}}{8\pi m_i} \right) = 0, \quad (26) $$

$$ \frac{\partial p}{\partial t} + \nabla \cdot (\mathbf{u} p) + (\gamma - 1) p (\nabla \cdot \mathbf{u}) = 0, \quad (27) $$

$$ \frac{\partial \mathbf{B}}{\partial t} + c \nabla \times \mathbf{E} = 0, \quad (28) $$

with auxiliary equations $\mathbf{E} = \mathbf{J} \times \mathbf{B}/nec - \mathbf{J}_i \times \mathbf{B}/nc$ and $\mathbf{J} = (c/4\pi) \nabla \times \mathbf{B}$. Here, $n \simeq \rho/m_i$ is the number density of ions (approximately equal to the number density of electrons due to quasi-neutrality) and $\mathbf{J}_i = n\mathbf{u}$ is the ion flux density.
In F3D, these equations are stepped forward using the trapezoidal leapfrog technique, a predictor-corrector method that is well-equipped to handle conservative equations (Zalesak 1979, 1981). Each of the above equations can be written as a conservative equation of the form (Guzdar et al. 1993)

$$ \frac{\partial \psi}{\partial t} + \nabla \cdot \mathcal{F} - D \nabla^2 \psi + F(\xi) = 0, \qquad (29) $$

where $\psi$ is the plasma variable in question, $\mathcal{F}$ is a suitably defined flux, $D$ is a second order diffusion coefficient which can be added to the equations to represent resistivity, viscosity, or a numerical dissipation to improve code stability, and $F(\xi)$ is a suitably defined sink/source term in terms of any other plasma variables $\xi$. (Equation (28) is not exactly in this form, but an analogous expression holds.)

To write the numerical algorithm, we use the standard notation where a superscript $n$ on a plasma variable refers to the time step in question, so the initial values are set at $n = 0$, the first time step is $n = 1$, and so on. To evolve $\psi^n$ to $\psi^{n+1}$ in a time step $\Delta t$, the trapezoidal leapfrog algorithm is

$$ \psi^{n+1/2} = \frac{\psi^{n-1} + \psi^n}{2} + \Delta t \left[ -\nabla \cdot \mathcal{F}^n + D \nabla^2 \psi^{n-1} - F(\xi^n) \right] \qquad (30) $$

$$ \psi^{n+1} = \psi^n \Delta t \left[ -\nabla \cdot \mathcal{F}^{n+1/2} + D \nabla^2 \psi^n - F(\xi^{n+1/2}) \right]. \qquad (31) $$

The first equation uses the data at the $n$'th time step and the data at the previous $n-1$'st time step to “predict” $\psi$ half a time step in the future. Then, the flux and source terms are evaluated at this intermediate time step to evolve $\psi$ the next half time step to the desired step $n+1$. To go from $n=0$ to $n=1$, data is needed at $n=-1$, which is simply taken to be the same as the data at $n=0$. This algorithm is second order in the time step $\Delta t$, meaning that the error from the algorithm is approximately a coefficient times the square of the time step. This is an example of an “explicit” time stepping algorithm because the data at the future time step is found completely using known data at the current or previous time steps.

The above shows how the equation is stepped forward temporally, but one also needs to calculate spatial derivatives. The approach F3D uses to calculate spatial derivatives is with a finite difference technique, which means the spatial derivatives are simply approximated by the derivative over the size of a grid scale instead of over an infinitesimal distance. For example, one approximation for the partial derivative in the $x$ direction on a grid with grid scale $\Delta x$ is $\partial \psi_j / \partial x \simeq (\psi_{j+1} - \psi_{j-1}) / 2\Delta x$, where the $j$ subscript refers to the index of the spatial cell for which the derivative is desired. This is a second-order scheme because the error relative to the exact derivative scales like $\Delta x^2$. The F3D code uses the approximation

$$ \frac{\partial \psi_j}{\partial x} \approx \frac{2}{3\Delta x} (\psi_{j+1} - \psi_{j-1}) - \frac{1}{12\Delta x} (\psi_{j+2} - \psi_{j-2}), \qquad (32) $$

which is fourth order (the error scales like $\Delta x^4$), and the cost for this higher order derivative is that it requires data from two adjacent cells on each side rather than one. Analogous expressions hold for spatial derivatives in the $y$ and $z$ directions. Second order derivatives are given by the associated fourth order accurate approximation $\partial^2 \psi_j / \partial x^2 \simeq \left[ -(2/3)(\psi_{j+1} - \psi_{j-1}) - (1/12)(\psi_{j+2} + \psi_{j-2}) \right] / (\Delta x)^2$.
The F3D code is written in Fortran 90 and is parallelized using Message Passing Interface (MPI) for use on supercomputers. The computational domain is rectangular with a fixed, regular grid. It can be run in two or three dimensions. When in two dimensions, the vectors can have an out of plane component even though all quantities are invariant in the out-of-plane direction; this is often referred to in the literature as “2.5 dimensional”. The F3D code does not explicitly enforce that $\nabla \cdot \mathbf{B} = 0$, but it has been demonstrated that the value is small when the initial magnetic field is divergence free. Numerous aspects of the F3D code are on user-controlled switches that can turn terms or effects on or off, and the initial plasma variable profiles and values are controlled by the user. Other features of F3D that are used to go beyond the Hall-MHD model will be treated in Sect. 3.6.

To run the F3D code, the user chooses the simulation domain size and desired grid scale to resolve the relevant physics, which is typically at least 5 times smaller than the relevant ion gyroradius. As F3D is an explicit finite difference code, the time step $\Delta t$ can be no larger than allowed by the so-called Courant-Friedrichs-Lewy (CFL) condition (Press et al. 1992), which requires $\Delta t \le \Delta x/v_{fastest}$, where $v_{fastest}$ is the fastest speed that can occur in the system. For Hall-MHD, the fastest speed is typically the whistler or kinetic Alfvén wave speed at or near the grid scale, but can be the fast magnetosonic speed for some ambient plasma conditions. 2D reconnection simulations with F3D are typically performed with time step about 40% of the CFL condition.

Any employed diffusion coefficients then need to be chosen. Often, the resistivity is not used for Hall-MHD, but a fourth-order diffusion numerical dissipation of the form $-D_4(\partial^4\psi/\partial x^4 + \partial^4\psi/\partial y^4 + \partial^4\psi/\partial z^4)$ on the right hand side of Eq. (29) is included in F3D to damp structures at the grid scale while minimally affecting larger scale structures. In order to preserve the five point stencil used for the other finite differences, second order accuracy in the fourth order derivatives is employed, with $\partial^4\psi_j/\partial x^4 \approx [\psi_{j+2} - 4\psi_{j+1} + 6\psi_j - 4\psi_{j-1} + \psi_{j-2}]/(\Delta x)^4$. The appropriate diffusion coefficient scales with $D_4 \sim v_{fastest}[\pi/(\Delta x)^3]$. An appropriate value for this coefficient is when it is large enough to control numerical issues at the grid scale while not impacting the large scale physics. A good approach to optimize this value is to run multiple simulations with only varying $D_4$, and finding a range of values for which the numerics are good and the large scale features are only weakly dependent on $D_4$; it is typically within an order of magnitude of the scaling prediction.

We have focused on F3D as an example of a Hall-MHD code because of its algorithmic simplicity and because it has long been used to study reconnection. There are drawbacks to the code. As a finite difference code, it does not capture shocks, and therefore if the number density gets fairly small in any given simulation the code typically crashes. The code performs well up to about 1000-2000 processors on high powered supercomputers; MHD codes without the Hall effect can be made to scale much better to 10s of thousands of processors. It is also restricted to a rectangular geometry with a regular grid.

There are a number of other Hall-MHD codes that have been used to study magnetic reconnection, some of which we gather here. Some codes used to study Hall-MHD reconnection in a rectangular domain have included VOODOO (Huba 2003), HMHD (Lottermoser and Scholar 1997), another code called HMHD (Huang et al. 2011), and the UI Hall-MHD code (Ma and Bhattacharjee 2001). Codes that have been used to study Hall-MHD in the context of planetary magnetospheres (including Earth’s) are a multi-fluid code (Winglee 2004), Block Adaptive Tree Solar-wind Roe Upwind Scheme (BATS-R-US) (Tóth et al. 2008), and Gkyell (Wang et al. 2018). Codes used in the tokamak geometry include NIMROD (Glasser et al. 1999) and M3D-C1 (Jardin et al. 2008).
## 3.6 Further Extensions of Hall-MHD

Here, we briefly discuss fluid model extensions beyond ideal-MHD that contain the Hall electric field and other terms. We start with terms that were dropped from Eq. (10).

### 3.6.1 Hall-MHD with Electron Inertia

One extension of Hall-MHD is to retain the electron inertia term $m_e du_e/dt$ from Eq. (10). We call this model “Hall-MHD with electron inertia;” it is often called the “two-fluid model” in the literature, but we refrain from this nomenclature since a completely different set of equations is also typically given the same name.

Using the electron inertia term as is would lead to a new variable $\mathbf{u}_e$ with a time derivative in the model. Instead, the standard approach is to recognize that the prefactor $m_e$ is small for electron-ion plasmas, so the only way this term important is if the electrons are moving very fast. The ions would be too slow to keep up in such a case, so on time scales where this term is important, we can treat the ions as approximately stationary. Then, the electron bulk flow velocity $\mathbf{u}_e$ is related to the current density via $\mathbf{u}_e = -\mathbf{J}/n_e e$. Since $\mathbf{J}$ and $n_e$ are already included in the Hall-MHD description, the set of equations remains closed. Such a simplification is convenient but not absolutely necessary. For example, Sect. 4.4 describes a hybrid code in which the ion flows and density effects in the electron inertia term are included.

Analytically, the inertial electric field from dividing Eq. (10) by $-e$ is given by $-(m_e/e)du_e/dt$, and replacing $\mathbf{u}_e$ by $-\mathbf{J}/n_e e$ gives $(m_e/e^2)d(\mathbf{J}/n_e)/dt$. Then, Ohm's law in Hall-MHD with electron inertia becomes

$$ \mathbf{E} + \frac{\mathbf{u} \times \mathbf{B}}{c} = \frac{\mathbf{J} \times \mathbf{B}}{n_e e c} + \frac{m_e}{e^2} \frac{d}{dt} \left( \frac{\mathbf{J}}{n_e} \right). \quad (33) $$

It is important to note that the same approximation $\mathbf{u}_e \approx -\mathbf{J}/n_e e$ is used in the convective derivative term, so that $d/dt \approx \partial/\partial t - (1/n_e e)\mathbf{J} \cdot \nabla$.

One needs to treat the $\partial/\partial t$ term on the right hand side of Eq. (33). To do so, we eliminate $\mathbf{E}$ in Faraday's law using Eq. (33); some algebra reveals that the equation becomes

$$ \frac{\partial \mathbf{B}'}{\partial t} = -c \nabla \times \mathbf{E}', \quad (34) $$

where $\mathbf{B}' = (1 - d_e^2 \nabla^2)\mathbf{B}$ is an auxiliary magnetic field and $\mathbf{E}' = \mathbf{J} \times \mathbf{B}'/n_e c - \mathbf{J}_i \times \mathbf{B}/n_e c$ is an auxiliary electric field, and the factor of $n$ in the inertia term is treated as a constant on the time scales of interest. Coupling this equation with Eqs. (25)-(27) gives a closed set of equations.

To solve these equations, we note that Eq. (34) looks just like the usual Faraday's law except for the primes, so the same numerical technique can be used to solve for $\mathbf{B}'$. Once $\mathbf{B}'$ is found, one uses that variable as the known source term in the equation $(1 - d_e^2 \nabla^2)\mathbf{B} = \mathbf{B}'$ to solve for $\mathbf{B}$. This is an elliptic differential equation with many algorithms that can be used to solve it. The simplest may be a relaxation technique (Press et al. 1992), but it is relatively slow, especially when used for codes that have been parallelized for use on a supercomputer. A faster version of relaxation is called multigrid (Trottenberg et al. 2000). F3D employs the Fast Fourier Transform approach.

We now briefly discuss the physics introduced by the electron inertia term and the advantages for including it in simulations. By comparing the electron inertia term to the Hall
term, we determine the condition under which it is important to retain the electron inertia term. We know the Hall electric field is important at scales below the ion inertial scale $d_i$, but vanishes at the X-line where $\mathbf{B} = 0$. The electron inertia term can be important at small scales, so we seek the scale at which it becomes comparable to the Hall electric field. Setting them equal gives $J_y B_x / nec \sim (m_e / e^2) [(\mathbf{J} / ne) \cdot \nabla] (J_z / n)$. In the scaling sense, we use $B_x \sim B_{up}$, $J_z \sim c B_{up} / (4\pi\delta)$ where $\delta$ is the scale at which the two terms are equal, so $J_y B_{up} / nec \sim (m_e / e^2) (J_y / ne\delta) (c B_{up} / (4\pi\delta n))$, which simplifies to $\delta^2 \sim d_e^2$, where $d_e^2 = m_e c^2 / (4\pi ne^2)$ is the electron inertial scale. Thus, at length scales below $d_e$, the electron inertia term can be important.

Consequently, one reason researchers include the electron inertia term into their Hall-MHD model is to aim to capture electron scale physics more accurately than without it. We point out, however, that the MHD model itself was derived with the assumption that $m_e / m_i$ is small, and therefore including the electron inertia term as we have done does not actually provide a self-consistent treatment of sub-$d_e$ scale physics. The two fluid model or the kinetic approach is needed to more accurately capture electron scale physics.

Then why include electron inertia? The answer is that it helps with the numerics. To see this, we consider waves in the Hall-MHD with electron inertia system, generalizing the treatment in Sects. 3.4.1 and 3.4.2. The dispersion relation for perfectly parallel propagating waves in Hall-MHD with electron inertia, generalizing Eq. (19), becomes

$$ \omega^2 = \frac{k^2 v_{A0}^2}{D_e} \left( 1 + \frac{k^2 d_{i0}^2}{2D_e} + \sqrt{\frac{k^2 d_{i0}^2}{D_e} + \frac{k^4 d_{i0}^4}{4D_e^2}} \right), \quad (35) $$

where $D_e = 1 + k^2 d_e^2$. When $kd_e$ is negligible, this dispersion relation reduces to the Hall-MHD result in Eq. (19). When $kd_e \gg 1$, the waves become electron cyclotron waves with $\omega^2 = \Omega_{ce}^2$, where $\Omega_{ce} = eB_0 / m_e c$. The reason this is useful numerically is that the whistler wave is dispersive, and in Hall-MHD the phase speed goes to infinity as the wavelength goes to zero. In Hall-MHD with electron inertia, the whistler rolls over to the electron cyclotron wave which does not propagate, so there is a maximum speed of the waves. This means that the time step required to run the simulation does not go to zero, and the simulations with very high resolution are less stiff and therefore cheaper to carry out.

### 3.6.2 Electron Magnetohydrodynamics (EMHD)

The electron-MHD (EMHD) model (Kingsep et al. 1990) is used when the large (MHD) scales are not of interest, and only the scales between electron and ion scales are of interest. In such a limit, the MHD terms concerning ion velocity are dropped. Because of this, the density and pressure can no longer change (on the time scales of interest), so the only remaining governing equation is Faraday's law, with the electric field given solely by the Hall electric field:

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}, \quad (36) $$

$$ \mathbf{E} = \frac{\mathbf{J} \times \mathbf{B}}{nec}. \quad (37) $$

Using Ampère's law $\mathbf{J} = (c/4\pi)\nabla \times \mathbf{B}$, it is common to combine these equations into a single equation for $\mathbf{B}$,

$$ \frac{\partial \mathbf{B}}{\partial t} = -\frac{c}{4\pi ne} \nabla \times [(\nabla \times \mathbf{B}) \times \mathbf{B}]. \quad (38) $$
This closed vector equation comprises the EMHD model. One can include electron inertia in a manner analogous to the full Hall-MHD model: $\partial \mathbf{B}' / \partial t = -(1/n_e e) \nabla \times (\mathbf{J} \times \mathbf{B}')$, where $\mathbf{J} = (c/4\pi) \nabla \times \mathbf{B}$ and $\mathbf{B}' = (1 - d_e^2 \nabla^2) \mathbf{B}$.

### 3.6.3 The Electron Pressure Term in Hall-MHD

We now consider the electron pressure term in Eq. (10). There, electron pressure is written as a tensor $\mathbf{p}_e$, which is the most general form following directly from the Vlasov/Boltzmann equation in kinetic theory. Retaining it in Ohm's law would give a term on the right hand side of Eq. (33) of the form $-(1/n_e e) \nabla \cdot \mathbf{p}_e$. We now consider a few commonly used approximations to simplify this.

**Isotropic Electron Pressure** In MHD, the total pressure is assumed to be isotropic, so that $\mathbf{p} = p\mathbf{I}$. In this case, the term entering Ohm's law is $-(1/n_e e) \nabla p_e$. When substituted into Faraday's law [Eq. (16)], this term gives a contribution of $(c/e) \nabla \times (\nabla p_e / n) = (c/e) \nabla (1/n_e) \times \nabla p_e = -(c/\mathbf{en}_e^2) \nabla n_e \times \nabla p_e$. This contribution to the electric field is called the “Biermann battery” (Biermann 1950) and it is often of great importance in reconnection in high energy density plasmas (Fox et al. 2012). In most space applications that employ a fluid model, the electrons are assumed to be adiabatic with $p_e/n_e^\gamma$ equal to a constant or isothermal with $p_e/n_e$ equal to a constant. In either limit, the Biermann battery term vanishes identically.

One may presume from this result that the electron pressure gradient term does not have any contribution to Hall-MHD, but this is not true. Including the scalar electron pressure gradient in the Hall-MHD Ohm's law gives

$$ \mathbf{E} + \frac{\mathbf{u}_e \times \mathbf{B}}{c} = -\frac{1}{n_e e} \nabla p_e, \qquad (39) $$

where we write $\mathbf{u}_e \approx \mathbf{u} - \mathbf{J}/n_e e$ for the bulk electron velocity. Taking the cross product of this equation with $\mathbf{B}$ and solving for the bulk flow velocity $\mathbf{u}_{e,\perp}$ perpendicular to the magnetic field $\mathbf{B}$ gives

$$ \mathbf{u}_{e,\perp} = c \frac{\mathbf{E} \times \mathbf{B}}{B^2} + \frac{c}{n_e e B^2} \nabla p_e \times \mathbf{B}. \qquad (40) $$

This equation implies the electron perpendicular bulk flow velocity is a sum of the $E \times B$ drift $c\mathbf{E} \times \mathbf{B}/B^2$ and the electron diamagnetic drift speed $\mathbf{u}_{e,\parallel} = c\nabla p_e \times \mathbf{B}/n_e e B^2$. This important result shows that the electron diamagnetic drift is captured in Hall-MHD provided the electron pressure is non-zero. An important implication of this is that magnetic flux convects at the sum of the $E \times B$ and electron diamagnetic drift speed when the electron pressure is retained in Ohm's law.

**Gyrotropic Electron Pressure** The next level of approximation for the electron pressure tensor is motivated by the fact that electrons in a magnetic field often have a different temperature in the directions parallel and perpendicular to the magnetic field, i.e., the electrons are gyrotropic. The electron pressure parallel to the magnetic field is $p_{e,\parallel}$ and perpendicular to the magnetic field is $p_{e,\perp}$. The general way to write a gyrotropic electron pressure tensor $\mathbf{p}_{e,g}$ for a magnetic field in an arbitrary direction $\hat{\mathbf{b}} = \mathbf{B}/B$ is (Parker 1957a)

$$ \mathbf{p}_{e,g} = p_{e,\perp}\mathbf{I} + (p_{e,\parallel} - p_{e,\perp})\hat{\mathbf{b}}\hat{\mathbf{b}}. \qquad (41) $$
Using this form in the generalized Ohm's law for the electron pressure, and an analogous term in the pressure gradient force in the momentum equation allows for the modeling of a plasma with a gyrotropic pressure.

Having a closed set of equations requires closures on the parallel and perpendicular pressures. The most widely known closure is the Chew-Goldberger-Low (CGL) closure, which assumes that there is no heat flux and the plasma is magnetized (Chew et al. 1956). A direct calculation using the Vlasov equation and these assumptions gives the CGL “double adiabatic laws”

$$ \frac{d}{dt} \left( \frac{p_{\perp}}{nB} \right) = 0, \quad \frac{d}{dt} \left( \frac{p_{\|} B^2}{n^3} \right) = 0, \qquad (42) $$

A second model of great importance to reconnection is the Egedal closure (Egedal et al. 2013). This arises in the upstream and downstream regions of magnetic reconnection regions. The key physics is that the magnetic fields in this region are mirror fields that can trap electrons. The presence of an electric field heats them parallel to the electric field, leading to elongated gyrotropic distributions in the parallel direction. For large magnetic field strength, the equations reduce to the isothermal equation of state. For small magnetic field strength, the equations reduce to the CGL equations. The two limits were interpolated to find a closure that could be implemented into a fluid code; see Egedal et al. (2013) for details.

# 3.7 Examples of Reconnection Simulation Results with the Hall Electric Field

## 3.7.1 Ion-Coupled Reconnection

The Hall-MHD model holds a significant place of importance in the history of reconnection simulations. From the earliest days of reconnection research, it was known that the Sweet-Parker model (Sweet 1958; Parker 1957b) was too slow to explain solar flares (Parker 1963), and it is also too slow to explain magnetotail reconnection (Parker 1973) and the sawtooth crash in tokamaks (Edwards et al. 1986; Yamada et al. 1994). It was discovered that using a localized resistivity leads to reconnection fast enough to explain the observed reconnection rates (Ugai and Tsuda 1977; Sato and Hayashi 1979), but the resistivity model was not derived from first principles and it remains inconclusive whether it can be. Then, reconnection with the Hall electric field was found to also produce rates comparable to observed values without relying on ad hoc terms (Aydemir 1992; Wang and Bhattacharjee 1993; Kleva et al. 1995; Ma and Bhattacharjee 1996). The “GEM Challenge” study ((Birn et al. 2001) and references therein), one of the most cited papers ever about reconnection, compared comparable simulations with different simulation models, and all models containing the Hall electric field led to fast reconnection with a reconnection rate approximately 0.1, as shown in Fig. 3. We note that this figure contains results from Hybrid and PIC simulations, which are described in Chaps. 4 and 5 respectively. It was only recently that an explanation of why the Hall electric field contributes to make the reconnection rate 0.1 was presented (Liu et al. 2022). Thus, Hall-MHD represents the minimal first-principles physics model that reproduces the reconnection rate achieved in kinetic models.

Another characteristic feature of collisionless reconnection is the quadrupolar structure of the out-of-plane magnetic field. For anti-parallel reconnection, it was posited that the Hall electric field would cause there to be a quadrupolar out-of-plane magnetic field due to the in-plane currents within the ion diffusion region (where the ions decouple from the magnetic field but the electrons remain frozen-in) (Sonnerup 1979). Although this is a completely nonlinear effect, it is analogous to the out-of-plane magnetic field generation in the
<figure><img src="image_5.png" /><figcaption>**Fig. 3** “GEM Challenge” result showing Hall-MHD simulations faithfully obtain the rate of change of reconnected flux obtained in kinetic models. The legend describes the simulation approach for each curve; resistive-MHD is far slower than the other models. Adapted from Birn et al. (2001)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Pressure and out-of-plane magnetic field in Hall MHD reconnection simulations. The legend describes the simulation approach for each curve; resistive-MHD is far slower than the other models. Adapted from Rogers et al. (2003)</figcaption></figure>

whistler wave shown in Fig. 2(b) (Drake and Shay 2007). It can also be thought of as the out-of-plane current dragging the reconnecting field out of the reconnection plane (Mandt et al. 1994). Hall-MHD simulations of reconnection produce a quadrupolar out-of-plane magnetic field (Huba and Rudakov 2002), as shown in the lower left plot in Fig. 4 from simulations in Rogers et al. (2003). When there is a strong out-of-plane (guide) magnetic field, the quadrupolar structure persists, but a quadrupolar structure in the gas pressure also arises with opposite polarity (Kleva et al. 1995). The physical reason is analogous to the gas pressure perturbation formation in kinetic Alfvén waves as shown in Fig. 2(c). These quadrupolar structures arise in Hall-MHD simulations of reconnection, as shown in the right two panels of Fig. 4 from a simulation with guide field three times as strong as the reconnecting magnetic field. Despite the ability to produce quadrupolar structure in these quantities, it is now known that the detailed structure of the quadrupolar structures is not precisely the same as observed in kinetic simulations (Shay et al. 2007; Karimabadi et al. 2007) or magnetospheric observations (Phan et al. 2007). It was shown that the inclusion of an electron pressure anisotropy leads to better agreement with kinetic modeling (Ohia et al. 2012). An electron pressure anisotropy with the CGL relations and without the Hall term can also produce reconnection rates comparable to Hall reconnection (Cassak et al. 2015).

Another aspect of Hall-MHD reconnection that is not captured in collisional (Sweet-Parker) reconnection is how reconnection that is localized in the out-of-plane direction
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Evolution of (a) out-of-plane electric field $E_y$ at the X-point ($x = 0, z = 0$) and (b) root-mean-square value of the normal component of magnetic field $B_z$ evaluated over the length of the reconnecting current sheet (between two outflow regions). Vertical dashed lines mark the times at which $E_y$ attains its peak value $E_y^{peak}$ for different values of the current sheet half-thickness $\epsilon$. (c) Scaling of $E_y^{peak}$ with $\epsilon$ (in log-scale) from simulations (blue circles) and a fit $E_y^{peak} = 0.05/\epsilon^{1.15}$ (red line). Adapted from Jain and Sharma (2015b)</figcaption></figure>

spreads in that direction. Spreading has been studied in a number of Hall-MHD studies (Huba and Rudakov 2002; Shay et al. 2003; Huba and Rudakov 2003; Karimabadi et al. 2004; Nakamura et al. 2012; Shepherd and Cassak 2012; Arencibia et al. 2021, 2023) and EMHD studies (Jain et al. 2013) and they agree with kinetic simulations of reconnection spreading (Lapenta et al. 2006). In particular, it was shown that anti-parallel reconnection does not spread in resistive-MHD, but it does spread in Hall-MHD and kinetic models (Nakamura et al. 2012; Arencibia et al. 2021).

## 3.7.2 Electron-Only Reconnection

Collisionless magnetic reconnection in electron scale current sheets has been investigated using fixed ion models (both fluid and kinetic) for several decades before the observational discovery of electron-only reconnection by MMS (Bulanov et al. 1992; Mandt et al. 1994; Drake et al. 1994, 1997; Shay and Drake 1998; Attico et al. 2000; Chacón et al. 2007; Zocco et al. 2009; Jain and Sharma 2009; Jain et al. 2012; Jain and Büchner 2015; Jain and Sharma 2015a,b). It was, however, not termed as electron only reconnection as it was not conceived at that time that reconnection without ion participation is possible. It was rather termed as “early phase of reconnection” (Jain and Sharma 2009; Jain et al. 2012; Jain and Sharma 2015a,b). Later, electron only reconnection in Earth’s magnetotail was also interpreted as an early phase of the standard reconnection (Hubbert et al. 2022; Farrugia et al. 2021) and there is ongoing discussion if the observed reconnection events without ion coupling are electron-only reconnection events or just an early phase of standard reconnection (Lu et al. 2020, 2022; Wang et al. 2020b; Yi et al. 2022).

Jain and Sharma (2009) were the first to propose that the early phase of reconnection after its onset in electron scale current sheets will be dominated by electron dynamics without coupling to ion dynamics and carried out EMHD simulations to study the physics of the early phase. This study was followed by further theoretical and simulation studies, both in 2-D and 3-D, using EMHD model and comparison with space observations by Cluster spacecraft (Jain et al. 2012; Jain and Sharma 2015a,b; Jain and Büchner 2014a,b, 2015; Jain et al. 2017b,a). These and other EMHD studies are also relevant for electron only reconnection.

Figure 5 shows results from the 2-D EMHD (x-z plane) simulations of magnetic reconnection in electron scale current sheets of different half-thicknesses ($\epsilon$) (Jain and Sharma 2015b). In these simulations, the equilibrium current density is $J = -n_0e u_{ey0} \hat{y} =$
$(B_0 c / 4\pi \epsilon) \text{sech}^2(z/\epsilon) \hat{y}$ corresponding to the equilibrium anti-parallel magnetic field $\mathbf{B} = B_0 \tanh(z/\epsilon) \hat{x}$. The density $n_0$ is uniform and thus the current is due to the electron flow $u_{ey0}$. The results are shown in the normalized units: length by electron inertial length $d_e = c/\omega_{pe} = c/(4\pi n_0 e^2 / m_e)^{1/2}$, magnetic field by $B_0$ and time by $\omega_{ce}^{-1} = (eB_0 / m_e c)^{-1}$. Reconnection rate, measured by out-of-plane electric field $E_y$ at the X-point and shown in Fig. 5a, reaches its peak value when the growth of the rms values of $B_z$ begins to slow down, consistent with the Faraday law which gives $\partial E_y / \partial x = -1/c \partial B_z / \partial t$. The peak reconnection rate $E_y^{\text{peak}}$ scales with $\epsilon$ as $E_y^{\text{peak}} = 0.05 / \epsilon^{1.15}$ and drops from $E_y^{\text{peak}} = 0.05 v_{\text{Ae}} B_0 / c$ to $E_y^{\text{peak}} = 0.01 v_{\text{Ae}} B_0 / c$ as the $\epsilon$ increases from $\epsilon = d_e$ to $\epsilon = 4 d_e$. This range of reconnection rate in units based on ion Alfvén speed is $E_y^{\text{peak}} = 0.43 - 2.15 v_{\text{Ai}} B_0 / c$ (using ion $m_i / m_e = 1836$) which is much larger than the value of the reconnection rate (0.1 $v_{\text{Ai}} B_0 / c$) for the standard ion-coupled reconnection. Reconnection rates for the electron only reconnection have been reported in the similar range by MMS observations (Burch et al. 2020) and PIC simulations (Sharma Pyakurel et al. 2019). A recent analytical model (Liu et al. 2025b) that incorporates the dispersive properties of Alfvén waves (Sect. 3.4) was proposed to explain these higher electron-only reconnection rates.

Note that the EMHD simulation results in Fig. 5 are independent of the strength of the guide magnetic field because a uniform guide field does not appear in 2-D EMHD equations (Jain and Büchner 2015). However, in 3-D, guide field can introduce current aligned instabilities in addition to the tearing instability as has been predicted by 3-D EMHD eigen value analysis (Jain and Büchner 2015) and simulations (Jain et al. 2017a). In 3-D, current aligned electron Kelvin-Helmholtz instabilities can grow in electron scale current sheets even in the absence of guide magnetic field (Jain and Büchner 2014a,b; Greess et al. 2021). These EMHD studies are relevant for the MMS observations of the electron shear flow generated electron Kelvin-Helmholtz vortices within the diffusion region of collisionless magnetic reconnection (Zhong et al. 2022, 2018; Hwang et al. 2019).

## 3.8 The Future of Hall-MHD

Interestingly, many computational plasma physicists in the reconnection community are moving away from the Hall-MHD model and its fluid extensions to study the small-scale properties of reconnection. The Hall-MHD model was crucial for understanding the minimal physics that gives rise to a 0.1 reconnection rate. As questions have moved to other aspects of reconnection including heating and particle acceleration, many researchers opt for kinetic models to more realistically capture the small scale physics than can be done with Hall-MHD. Treatments of particle acceleration and heating at large scales (Arnold et al. 2021) do not require the Hall electric field.

As example of an avenue of modern reconnection research and modeling where the Hall-MHD model remains highly beneficial is in the MHD-EPIC approach to global magnetospheric modeling (Daldorff et al. 2014), as is discussed more fully in Sect. 6. In this approach, the fluid model is used in regions of the magnetosphere where no small important scale physics takes place, which allows for faster run times. In regions where small scale physics does take place, the code couples to a particle-in-cell code that captures this physics. The numerical results between the two models are passed back and forth to each other across their boundaries. The Hall-MHD model is well suited to be used in the transition region between the PIC and MHD models to facilitate a more accurate transition between the two models. The Hall-MHD model has been used to study global magnetospheric systems for Earth and other planets and moons to great effect (Paty and Winglee 2004; Dorelli et al. 2015; Dong et al. 2019; Li et al. 2023b), and it is anticipated that further advances will
continue to be made with the Hall-MHD approach for systems too large to employ global kinetic codes.

# 4 Hybrid Simulations

The Earth's magnetosphere is a complex plasma system characterized by a multitude of multiscale processes governing the interaction of the solar wind with the Earth's dipole magnetic field. Modeling small-scale turbulent processes in the foreshock and magnetosheath requires inclusion of ion kinetic effects and Hall physics (Karimabadi et al. 2014; Omelchenko et al. 2021a). Furthermore, kinetic treatment of hot and cold ion populations is greatly needed for improved modeling of ionospheric outflows and their impact on the magnetopause and magnetotail. To study magnetic reconnection, one also needs to incorporate finite electron-mass effects (e.g. Biskamp 2000; Birn and Priest 2007; Gonzalez and Parker 2016), or mimic these effects with *ad hoc* (resistivity) models.

The necessity to account for smaller and faster scales in kinetic simulations in a manner that would guarantee their numerical accuracy and computational efficiency creates challenges in global modeling of the Earth's magnetosphere. Since describing plasma kinetics with pure “first-principles” models is still not feasible, various approximations have been developed as candidates for future “beyond MHD” operational modeling, with multiple levels of physical fidelity included. However, to what degree kinetic processes may influence the “fluid-like” behavior of the magnetosphere on global scales still remains an open question. As we argue below, many of these challenges can be addressed by self-consistent hybrid modeling, where Maxwell's equations are solved in the quasi-neutral Darwin limit, ion species are treated kinetically, and the plasma electrons are approximated as an inertialess fluid. These hybrid models can be broken into two categories, according to the computational techniques used to represent kinetic ions: Particle-in-Cell (PIC) models (Winske et al. 2003; Lipatov 2002) and Vlasov models (von Alfthan et al. 2014). In what follows we discuss only the hybrid-PIC approach because it has already been applied successfully to perform three-dimensional (3D) simulations of global plasma systems that range from the Earth's magnetosphere (e.g. Lin and Wang 2005; Omelchenko et al. 2021a), planets (Herčík et al. 2013; Jarvinen et al. 2020), and small space bodies (Fatemi et al. 2017; Kallio et al. 2019) to compact laboratory plasmas (Omelchenko and Sudan 1997; Lin et al. 2008; Thoma et al. 2013; Omelchenko 2015; Omelchenko and Karimabadi 2022). The hybrid-Vlasov approach (von Alfthan et al. 2014) is relatively new and considerably more computationally expensive, with production runs being still restricted to quasi-3D setups (Pfau-Kempf et al. 2020).

An important issue to grasp reconnection correctly is the consideration of the finite electron inertia, as it has been shown by EMHD simulations. This could reveal not only the properties of electron-only reconnection, proposed by (Jain and Sharma 2009) and recently discovered by MMS in the magnetosheath (Phan et al. 2018), but also the transition from electron- to ion reconnection. These effects can be treated by hybrid-approach with kinetic ions and an inertial electron fluid (Sect. 4.4).

## 4.1 Model Equations: Massless Electrons

The standard hybrid model (Winske et al. 2003) assumes plasma quasi-neutrality, neglects the displacement current in Maxwell's equations, and treats ions as full-orbit macro-particles (in the PIC approach) moving in self-consistent electric and magnetic fields. The plasma
electrons are approximated as an inertialess fluid with scalar pressure described by either an adiabatic law or evolution equation. Together with a self-consistent PIC method for the ion components, this leads to a set of hybrid equations that include Ampere's law in the magnetostatic limit, Faraday's law, and an algebraic expression for electric field (generalized Ohm's law) with the Hall, electron pressure gradient, and resistive terms (e.g. Omelchenko and Karimabadi 2012):

$$ \frac{d\mathbf{x}_i}{dt} = \mathbf{v}_i, \qquad (43) $$

$$ m_i \frac{d\mathbf{v}_i}{dt} = q_i (\mathbf{E} + \frac{\mathbf{v}_i \times \mathbf{B}}{c}), \qquad (44) $$

$$ \nabla \times \mathbf{B} = \frac{4\pi}{c} \mathbf{J}, \quad \mathbf{J} = \mathbf{J}_e + \mathbf{J}_i, \qquad (45) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}, \qquad (46) $$

$$ \mathbf{E} = \frac{\mathbf{J}_e \times \mathbf{B}_t}{en_e c} - \frac{\nabla p_e}{en_e} + \eta \mathbf{J}, \quad \mathbf{B}_t = \mathbf{B} + \mathbf{B}_{ext}, \qquad (47) $$

$$ en_e = \rho_i, \qquad (48) $$

$$ p_e = n_e T_e \sim n_e^\gamma. \qquad (49) $$

Eqs. ((43) and (44)) are the equations of motion for each ion particle. $q_i$ is the charge on each ion, and in these equations the electric and magnetic fields are interpolated from the grid onto the each particle's position. In Eqs. (47)-(49) $n_e$, $\mathbf{J}_e$ are the electron number and current density, respectively; $p_e$ is the electron pressure, here assumed to governed by Eq. (49)) with an adiabatic constant of $\gamma$; $T_e$ is the electron temperature; $\rho_i$, $\mathbf{J}_i$ are the total ion charge density and current density (found by summing up individual particles around each grid point); $\mathbf{E}$ is the electric field; $\mathbf{B}$, $\mathbf{B}_{ext}$ are the “self-generated” ($\mathbf{B}|_{t=0} = 0$) and “external” (steady state) magnetic fields, respectively.

The applied plasma resistivity, $\eta$ is either constant or chosen to be a function of plasma parameters (e.g. Lin et al. 2007; Omelchenko et al. 2021a). The resistive term in the generalized Ohm's law (Eq. (47)) may (i) describe finite conductivity of plasma or space bodies (e.g., the Moon (Fatemi et al. 2017; Omelchenko et al. 2021b)), (ii) imitate finite electron inertia effects in magnetic reconnection events, and (iii) enable fast magnetic field diffusion at low-density (“vacuum”) cells, $n_e \le n_{min}$, where $n_{min}$ is a small density cutoff value (Omelchenko et al. 2021b). Failure to properly treat low-density regions in hybrid simulations may lead to non-physical results (Omelchenko et al. 2021b; Poppe 2019).

## 4.2 Key Physics Beyond MHD in Hybrid Models

The “mesoscale” hybrid model occupies the middle ground between the “large-scale” fluid and “micro-scale” first-principles modeling paradigms. For global magnetospheric simulations, the hybrid model enables a number of “beyond MHD” capabilities, as explained below.

**Modeling turbulent processes in the foreshock and magnetosheath.** Unless large resistive damping (or smoothing) is applied, the hybrid model accurately captures the Hall
physics for mesh cell sizes, $\Delta \le d_i$, where $d_i = c/\omega_{pi}$ and $\omega_{pi}$ are the local ion inertial length and plasma frequency, respectively. The Hall effects phase out on coarser meshes, $\Delta \gg d_i$, where the Alfvén term becomes greater than the Hall term in Eq. (47) and the whistler mode frequency, $\propto \Delta^{-2}$ becomes lower than the Alfvén mode frequency, $\propto \Delta^{-1}$. In fact, in this case the Hall term can completely be removed from the electric field in Faraday’s law (Eq. (46)) and kept only in the equations of ion motion (Karimabadi et al. 2004). The ability of a hybrid code with the Hall term to run stably on coarser meshes ($\Delta \gtrsim d_i$) may depend on the numerical implementation of Faraday’s law (Omelchenko and Karimabadi 2012).

Global hybrid codes have been used to address the ultra-low frequency (ULF) physics of the curved bow shock on the ion inertial/Larmor radius scales as the physics of the bow shock is predominantly determined by kinetic physics associated with charged particles from the solar wind. Of particular interest are the foreshock waves and diffuse ion distributions (Wang et al. 2009) and transient perturbations originating from the wave-particle processes in the quasi-parallel shock or due to the shock interaction with incoming solar wind discontinuities, including hot flow anomalies (Lin 2002; Lin et al. 2022b), foreshock bubbles (Omidi et al. 2010; Wang et al. 2020a), foreshock cavities (Lin and Wang 2005; Blanco-Cano et al. 2011), and high-speed jets (Omelchenko et al. 2021a; Palmroth et al. 2018b). The 3D hybrid simulations with ANGIE3D (AuburN Global hybrId CodE in 3-D) link the foreshock perturbations to the surface perturbations and kinetic-scale shear Alfvén waves (KAWs) at the magnetopause through mode conversion from the incoming compressional waves (Lin and Wang 2005; Shi et al. 2013), as well as the subsequent excitation of toroidal-mode field line resonances in the magnetosphere (Shi et al. 2021). It has also been shown that 3D models are essential for addressing the nonlinear physics of mode coupling and ion diffusion at the magnetopause (Lin et al. 2012).

The whistler mode plays a significant role in regulating turbulence in the magnetosheath and mediating magnetic reconnection in the Earth’s magnetosphere (e.g. Dorelli and Birn 2003; Drake et al. 2008). Hybrid simulations generally have to resolve the quadratic dispersion of this mode, $\omega \propto k^2$. This requirement may create computational bottlenecks in simulations of strongly inhomogeneous magnetospheric and laboratory plasmas, where whistler timescales typically span several orders of magnitude (Omelchenko and Karimabadi 2012, 2022). If not accurately integrated in time (or resistively damped), the spurious short-wavelength oscillations may grow explosively unstable from noise and terminate simulation (Lin et al. 2008). It should also be noted that although particle noise in hybrid-PIC simulations typically degrades their physical resolution, the Lagrangian (particle) approach enables transport of ion species with less numerical diffusion compared to the Eulerian approach to solving the Vlasov equation on velocity meshes (von Alfthan et al. 2014).

**Collisionless reconnection at the magnetopause and in the tail plasma sheet.** The physics of magnetic reconnection in the magnetosphere can be investigated by carrying out global hybrid simulations with an ad-hoc current-dependent resistivity. For the dayside magnetopause, the modeling topics include the structures of ion diffusion region and outflow regions (Tan et al. 2011), global evolution of flux transfer events (FTEs) and magnetic flux ropes (Omidi and Sibeck 2007; Guo et al. 2020, 2021d), propagation of kinetic Alfvén waves and Poynting flux from reconnection (Wang et al. 2019), ion cusp precipitation and energy spectrum (Omidi and Sibeck 2007; Tan et al. 2012), the triggering of reconnection by solar wind discontinuities (Omidi et al. 2009; Pang et al. 2010; Guo et al. 2021c), and magnetosheath turbulence (Chen et al. 2021; Ng et al. 2021). Likewise, high-latitude reconnection tailward of the cusp under northward IMF has also been simulated (Lin and Wang 2006; Guo et al. 2021e). The 3D global physics of storm-time magnetotail reconnection, fast flow
and entropy bubbles, the Hall-effects control of dawn-dusk asymmetry (Lin et al. 2014; Lu et al. 2016; Lin et al. 2017), and the associated global Alfvénic coupling between the magnetotail and the ionosphere under southward IMF have been simulated using the ANGIE3D code (Cheng et al. 2020). An Attempt has also been made to investigate the subsequent connection of fast flows to the ring current and radiation belt by combining ANGIE3D with the Comprehensive Inner Magnetosphere-Ionosphere (CIMI) model (Lin et al. 2021).

**Inclusion of multi-species plasma ion populations of solar wind origin and improved representation of ionospheric outflow.** In general, global hybrid models may include multiple ion species for representing solar wind and ionospheric outflow plasmas. Ionospheric outflow ions should be treated kinetically and self-consistently in order to properly account for their impact on the Earth's magnetosphere. Multi-fluid MHD models do not account for ion resonance acceleration and cyclotron effects, especially for heavy ions (Toledo-Redondo et al. 2021). Self-consistent 3D hybrid simulations of the impact of oxygen outflow on the magnetotail configuration and stability have recently been performed with the HYPERS code (Mouikis et al. 2021; Omelchenko et al. 2023).

**Modeling local reconnection and electron scale physics.** In MHD simulations, magnetic reconnection is often a result of mesh-dependent diffusion that is difficult to control numerically. The hybrid-PIC model is inherently more robust in this regard because ions are modeled as Lagrangian particles. As a result, reconnection onset and dynamics are controlled by the Hall physics and parameter-dependent resistivity. The hybrid-PIC model is also known to accurately reproduce reconnection rate when the ion inertial and cyclotron scales are properly resolved (Stanier et al. 2015). Further modifications of the hybrid model, which for instance may incorporate finite electron mass effects (e.g. Omelchenko et al. 2021c), could increase physical fidelity of reconnection modeling in the future.

**Modeling non-MHD waves in a global context.** The standard hybrid model supports ion cyclotron, whistler, and kinetic-Alfvén wave modes, which play an important role in regulating plasma turbulence in the Earth's magnetosphere and impact its global behavior. Modern observations report streams of non-Maxwellian ions that excite plasma turbulence through numerous kinetic instabilities that cannot be modeled within MHD. Global hybrid modeling naturally incorporates the ion kinetic effects into global models of the Earth's magnetosphere, which helps advance our understanding of the effects of turbulent plasma dynamics on global physical processes.

**Other applications.** 3D hybrid codes in space physics have been used to simulate shock-driven ion acceleration (Caprioli 2014; Guo and Giacalone 2013), solar wind turbulence (Franci et al. 2018; Roytershteyn et al. 2015), the Moon's wake (Fatemi et al. 2017; Kallio et al. 2019; Omelchenko et al. 2021b), planetary magnetospheres (Jarvinen et al. 2020), and small space bodies (Alho et al. 2019; Delamere 2009). Comparing results from full-scale 3D hybrid simulations of small space bodies with satellite observations provides yet another important route for validation and further extension of the hybrid approach to plasma modeling. Importantly, the recent advances in ionosphere-magnetosphere coupling (Lin et al. 2022b), code optimization (Dong et al. 2021) and multiscale computing (Omelchenko et al. 2021a) have greatly improved the prospects for hybrid simulations to become a key factor to consider in the overall theory of global solar wind-magnetosphere-ionosphere interactions.

## 4.3 Numerics

Equations (47)-(49), together with the self-consistent equations of motion for the ion macro-particles, may present computational challenges when used for modeling complex 3D plasma systems, such as the Earth's magnetosphere. Below we discuss some recent computational advances aimed at overcoming these problems.
### 4.3.1 Spatial Scales

Hybrid-PIC simulations typically intend to resolve the spatial scales of the order of the ion inertial length, $d_i$ and ion cyclotron radius, $r_{ci}$. The physical validity regime of the hybrid model ranges from large MHD scales down to $kr_{ci} \sim 1$ and $\omega t \sim 1$. The actual physical resolution of a hybrid simulation is largely determined by (i) how well these characteristic lengths are resolved on a mesh, (ii) how many particles are used. For the typical solar wind proton inertial length, $d_i \sim 100$ km, the Earth's radius, $R_E \sim 60d_i$. To encompass the whole magnetosphere, the computational mesh in a global simulation should cover the magnetopause with a typical standoff distance, $R_{MP} \sim 10R_E \sim 600d_i$ and the magnetotail stretching from the Earth to far distances, $R \sim 100R_E \sim 6000d_i$. The need to accurately account for the “far-field” inflow and outflow boundary conditions in the presence of a magnetic dipole may additionally require multiplying these dimensions by a factor of 2-3. Approximating such large 3D domains with uniform meshes with cell sizes of the order of $\sim 1d_i$ is prohibitively expensive for the hybrid model because of the need to advance ions and fields at all cells on kinetic scales.

To overcome these restrictions, several options are available. First, one may increase the cell size beyond $1d_i$ at the expense of lower accuracy in resolving the Hall physics (Omelchenko and Karimabadi 2012). Second, one may artificially increase the ratio between the solar wind ion inertial length $d_i$ ($d_i \sim r_{ci}$ for the outer magnetosphere regions with ion $\beta \sim 1$) and the magnetopause distance $R_{MP}$, in order to better accommodate the available computation resources while still choosing a sufficiently large value of $R_{MP}/d_i$ for assuring the separation between the global and local-kinetic scales (Omidi et al. 2004). Both approaches efficiently “downscale” the Earth’s magnetosphere. For instance, ANGIE3D (Lin et al. 2014; Lin and Wang 2005) does it by artificially reducing the solar wind plasma density, i.e. by inflating the characteristic inertial ion length and proportionally increasing the Alfvén speed. Alternatively, H3D (Karimabadi et al. 2014), hybrid-VPIC (Dong et al. 2021), and HYPERS (Omelchenko and Karimabadi 2012) employ the physical ion inertial length but scale the realistic magnetopause standoff distance down by a factor of 4-6 by using a weaker magnetic dipole.

Regardless of a chosen magnetosphere scaling method, present-day 3D hybrid codes typically use $R_{MP}/d_i \ge 100$. One of the largest 3D hybrid simulations to date was performed with HYPERS for $R_{MP}/d_i \simeq 160$ on a uniform mesh with approximately 1000 $\times$ 2000 $\times$ 2000 cells (Omelchenko et al. 2021b). To speed up global 3D simulations, hybrid codes may employ nonuniform meshes, among which ‘stretched’ (logically mapped) Cartesian meshes are the simplest. Nonuniform meshes typically maintain high resolution in a central domain of interest, while expanding cells towards domain boundaries (Omelchenko et al. 2021a; Lin and Wang 2005) to guarantee that the dipole field vanishes at the inflow/outflow (GSM X) and lateral (GSM Y and Z) boundaries so that robust local boundary conditions can be implemented. Sometimes, for simplicity, the lateral domain boundaries may be assumed to be periodic (e.g. Turc et al. 2015; Müller et al. 2011). This simplification, however, makes a global simulation valid for shorter simulation periods, until reflected particles or electromagnetic perturbations reach the periodic boundaries. To improve the counting statistics for macro-particles, splitting techniques may be used to enhance energetic particle distributions in the dayside magnetosphere (Omelchenko et al. 2021a) and magnetotail (Lin et al. 2007; Wang et al. 2009).

To further reduce the number of computational cells in global simulations, one may employ curvilinear (e.g., spherical) meshes (e.g. Dyadechkin et al. 2013; Guo et al. 2021b).
For example, to capture the short-wavelength physics along the shock normal and the magnetopause, early hybrid simulations, focusing on the dayside regions, employed cylindrical (2D) (Swift 1996; Lin et al. 1996; Lin 2002) and spherical (3D) (Lin and Wang 2005) coordinate systems. The spherical coordinate lines, however, have a singularity on the polar axis, which was handled by rotating the polar coordinates to the equator and omitting a conic region around them, while keeping the physical polar regions inside the domain (Lin and Wang 2005). A 2D hybrid simulation of the magnetotail also used curvilinear coordinates to accommodate the tail geometry (Swift and Lin 2001; Lin 2002). Similar to the need of assuring proper numerical resolution for resolving the kinetic scales along the curved or oblique boundary surfaces on the Cartesian meshes, special care is also necessary for the curvilinear meshes, especially when their coordinate lines are not orthogonal (Swift and Lin 2001). Compared to Cartesian meshes, curvilinear meshes may introduce additional discretization errors due to their (i) typically lower orders of numerical approximation, (ii) anisotropic particle-mesh weighting. These errors lead to various numerical artefacts and non-conservation of particle momentum ("self-forces"). As a result, it is necessary to benchmark results from simulations obtained with curvilinear meshes with similar simulations performed with Cartesian meshes (Dyadechkin et al. 2013).

Some global hybrid codes employ adaptive mesh refinement (AMR) (Leclercq et al. 2016; Müller et al. 2011). Hybrid AMR simulations, however, may suffer from spurious particle "self-forces" and wave reflections that occur at the mesh refinement interfaces. To mitigate these artefacts, AMR algorithms are typically complemented with smoothing procedures, which, however, should be performed with caution in order to avoid affecting underlying physics in the regions of interest.

### 4.3.2 Temporal Scales

In addition to the “slow” Alfvénic (MHD) timescales, hybrid codes need to follow the “fast” ion kinetic, ion cyclotron and whistler timescales. This requirement typically makes global hybrid simulations numerically “stiff” in the near-Earth space, where timesteps, required for numerical accuracy and stability, may become prohibitively small (Omelchenko and Karimabadi 2012). To partially mitigate these effects, the kinetic ions may be replaced in this region by a dense fluid (Swift 1996; Lin et al. 2021). Hybrid-PIC simulations inherently generate spurious oscillations with large wave numbers, $k \sim 1/\Delta$ and high frequencies, $\omega \sim k^2 \sim 1/\Delta^2$. If not properly integrated or resistively damped, these noisy oscillations may explosively grow and abort simulation (Lin et al. 2008). Deleterious instabilities may be avoided by applying “noise filtering” (smoothing) or/and various “flux-limiting” techniques for electric and magnetic fields. These modifications, however, need to be implemented with caution, as the may produce artificial solutions not supported by the hybrid model.

For accuracy, typical full-orbit particle solvers (“pushers”) require that time steps, $\Delta t_p$, should be small enough that $\Omega \Delta t_p \ll 1$, where $\Omega$ is the local ion gyro-frequency. Using the same time step for all particles may create another numerical bottleneck in global hybrid simulations of the Earth’s magnetosphere. For instance, particle time steps of the order of $\Omega_0 \Delta t_p \sim 0.05$ (where $\Omega_0$ is the ion gyro-frequency computed with respect to the IMF strength, $B_{IMF}$) fairly well describe ion gyro-motion in the solar wind (e.g. Turc et al. 2015; Guo et al. 2021b). At the same time, gyro-orbits and drifts of ions with $\Omega \gtrsim 10 \Omega_0$ (e.g. found in the cusp or some parts of the magnetosheath), will not be reproduced with accuracy. As a remedy, in some simulations, sub gyro-orbit time steps may be employed in these (large magnetic field) regions of the magnetosphere (Lin et al. 2014).

To summarize, predicting optimum global time steps for the particles and fields in global hybrid simulations is difficult in practice. This challenge has been addressed by replacing
time stepping with an asynchronous approach to time integration, which combines discrete-event simulation (DES) with elements of artificial intelligence: Event-driven Multi-Agent Planning System (EMAPS) (Omelchenko and Karimabadi 2006b, 2022). EMAPS enables time advance of individual particles and local fields on meshes of arbitrary topology by integrating them on their self-adaptive timescales, similar to Conway's Game of Life, where simulation elements evolve asynchronously based on a set of local interaction rules, rather than being updated synchronously at global time steps (Gardner 1970). Thus, EMAPS effectively performs the role of an intelligent "simulation time operating system". This approach was first applied to model 1D collisionless plasma shocks (Omelchenko and Karimabadi 2006a) and fluids (Omelchenko and Karimabadi 2006b, 2007). Implemented in HYPERS (Omelchenko and Karimabadi 2012), EMAPS has enabled efficient and accurate global 3D hybrid simulations of the Earth's magnetosphere (Omelchenko et al. 2021a,b, 2023).

Global 3D hybrid simulations of the Earth's magnetosphere are typically performed for simulation periods, $\Omega_0 t \sim 100-500$, where $\Omega_0$ is the IMF based proton cyclotron frequency. Assuming $B_{IMF} = 5 nT$, these simulations formally span relatively short (compared to MHD) magnetospheric times, $t < 20$ min. For convenience, in order to present physical results in "magnetospheric hours", some modelers multiply this simulation time by a model scaling factor (Lin et al. 2022b). Although this scaling is useful for comparing "macro-scale" simulation phenomena with observations, it is not appropriate for describing ion kinetic effects, e.g. those that drive the "magnetokinetic" formation of high-speed jets (Omelchenko et al. 2021a).

### 4.3.3 Plasmasphere and Ionosphere

Currently, hybrid codes cannot model global magnetospheric convection lasting many hours or days, e.g. a steady-state process of magnetotail loading and unloading. Therefore, global hybrid models typically assume a simple perfectly conducting or resistive ionosphere, where dipole magnetic field lines are "tied up" to the inner boundary (zero electric field) or allowed to diffuse due to its finite resistivity, respectively. model (e.g. Lin et al. 2021),

To avoid computing fast kinetic timescales, a cold, incompressible, dense ion fluid may be assumed to co-exist together with low-density particle ions in the inner magnetosphere within the distance of plasmasphere, where the plasma density is high (Swift 1996; Lin and Wang 2005; Lin et al. 2014). In ANGIE3D, this region is bounded by the near-Earth (inner) boundary, which is located at a radial distance at $r \approx 3.5 R_E$ in the inner magnetosphere. The field-aligned currents, calculated near this inner boundary and mapped along the geomagnetic dipole field lines down to the ionospheric altitude (1000 km), are used as input to the ionospheric potential equation solved on a sphere (Lin et al. 2014, 2021):

$$ \nabla \cdot (-\boldsymbol{\Sigma} \cdot \nabla \Phi) = J_{\parallel} \sin I, \qquad (50) $$

where $\boldsymbol{\Sigma}$ is the conductance tensor, $\Phi$ is the electric potential, $J_{\parallel}$ is the mapped field-aligned current density, and $I$ is the inclination of the dipole field at the ionosphere. The static analytical model of Hall and Pederson conductance that accounts for EUV and diffuse auroral contributions can be used for the conductance tensor. Similarly to the global MHD models, the ionospheric electric field is mapped along the dipole lines back to the inner magnetospheric boundary, to serve as a boundary condition for the cold ion fluid (Lin et al. 2021).

### 4.4 Hybrid Simulations with Inertial Electrons: Hybrid-PIC Code CHIEF

The MMS mission investigates physical processes like magnetic reconnection, shock waves and turbulence, which span from ion to electron scales. Simulation studies of these processes
should ideally cover full kinetic physics from ion to electron scales for which the necessary present and near-future computational resources are prohibitively expensive. Therefore simulation models, which cover different scale ranges and physical phenomena, are used.

Hybrid-kinetic plasma simulation model, introduced in the previous Sect. 4.1, treats ions as kinetic species and electrons as a massless fluid. This restricts their applicability to physical processes in which not only electron kinetic effects are not important but also to the scales exceeding by far the electron scales. Hybrid-kinetic codes with inertia-less electrons, discussed in Sect. 4.1, can, therefore, be used to simulate global phenomena and in some cases for specifically limited physics studies of magnetic reconnection, plasma turbulence and shock waves.

The validity of hybrid-kinetic model can, however, be extended down to electron length scales, viz., to electron inertial length by considering electrons as an inertial fluid (Jain et al. 2023). Since the electron kinetic physics is still ignored such plasma model might computationally be more feasible than the fully kinetic model and describe larger scale phenomena and plasma process like magnetic reconnection, plasma turbulence and shock formation in collisionless plasmas, in which electron scale structures develop.

### 4.4.1 Model Equations: Inertial Electrons

Hybrid-kinetic model treats ions as kinetic species and electrons as an inertial fluid. In hybrid-kinetic simulation codes, ion dynamics can be described by solving either the ion's Vlasov equation using Eulerian methods or the equations of motion for ion macro-particles using semi-Lagrangian Particle-in-Cell (PIC) method. Solving Vlasov equation is computationally more expensive. Here we discuss the hybrid-PIC codes which treat ions as Lagrangian macro-particles modelled via the PIC method. Following are the governing equations of hybrid-PIC model.

$$ \frac{d\mathbf{x}_i}{dt} = \mathbf{v}_i, \qquad (51) $$

$$ m_i \frac{d\mathbf{v}_i}{dt} = e(\mathbf{E} + \frac{\mathbf{v}_i \times \mathbf{B}}{c}), \qquad (52) $$

$$ \nabla \times \mathbf{E} = -\frac{1}{c} \frac{\partial \mathbf{B}}{\partial t}, \qquad (53) $$

$$ \nabla \times \mathbf{B} = \frac{4\pi}{c} \mathbf{J}, \qquad (54) $$

$$ \mathbf{J} = e(n_i \mathbf{u}_i - n_e \mathbf{u}_e), \qquad (55) $$

$$ n_i = n_e \qquad (56) $$

$$ \mathbf{E} = -\frac{\mathbf{u}_e \times \mathbf{B}}{c} - \frac{1}{en} \nabla p_e - \frac{m_e}{e} \left( \frac{\partial \mathbf{u}_e}{\partial t} + (\mathbf{u}_e \cdot \nabla)\mathbf{u}_e \right) + \eta \mathbf{J}, \qquad (57) $$

$$ p_e = C n_e^\gamma \qquad (58) $$

The electric and magnetic fields (E and B respectively) in Maxwell's equations, Eqs. (53) and (54), are coupled to the plasma dynamics via the total current density $\mathbf{J} = e(n_i \mathbf{u}_i - n_e \mathbf{u}_e)$ resulting from the bulk motion of ions (number density $n_i$, bulk velocity $\mathbf{u}_i$) and electrons (number density $n_e$, bulk velocity $\mathbf{u}_e$). Ion's number density $n_i$ and the bulk velocity $\mathbf{u}_i$ is obtained from their positions $\mathbf{x}_i$ and velocities $\mathbf{v}_i$ governed by Eqs. (51) and (52). Electron
dynamics is governed by quasi-neutrality condition, Eq. (56), momentum equation of the inertial electron fluid, Eq. (57), and equation of state relating electron scalar pressure $p_e$ with electron number density $n_e$, Eq. (58). Here, $e$ is the fundamental charge, $m_i$ ion mass, $m_e$ electron mass, $\eta$ collisional resistivity, $\gamma$ the adiabatic constant and $C$ is a constant (to be determined from initial conditions). Equations (51)-(58) are the fundamental equations of the hybrid-kinetic model with inertial electron fluid (Jain et al. 2023). These equations differ from the equations of hybrid-kinetic model with inertia-less electron fluid only by the electron inertial terms proportional to $m_e$ on the RHS of Eq. (57). Addition of electron inertial terms in Eq. (57) makes the numerical solution of these equations much more involved in comparison to the case of inertia-less electron fluid. The algebraic calculation of electric field from Eq. (57) is not as straightforward as in the case of the inertia-less electron fluid. One needs to now calculate time derivative of $\mathbf{u}_e$ or find some other way to obtain electric field. The calculation of magnetic field also now requires numerical solution of additional elliptic partial differential equations arising because of the finite electron inertia.

In majority of the hybrid-kinetic codes with inertial electrons, evolution of magnetic field is followed by solving an evolution equation for the generalized vorticity $\mathbf{W} = \nabla \times \mathbf{u}_e - e\mathbf{B}/m_e c$ obtained by taking curl of Eq. (57) and using Eq. (53). This equation is,

$$ \frac{\partial \mathbf{W}}{\partial t} = \nabla \times [\mathbf{u}_e \times \mathbf{W}] - \nabla \times \left( \frac{\nabla p_e}{m_e n} \right) - \nabla \times \left( \frac{e \eta}{m_e} \mathbf{J} \right). \quad (59) $$

The magnetic field is then calculated by solving an elliptic partial differential equation (PDE) which is obtained by substituting for $\mathbf{u}_e$ from Eq. (54) and (55), $\mathbf{u}_e = \mathbf{u}_i - c\nabla \times \mathbf{B}/(4\pi en)$, in the expression for $\mathbf{W} = \nabla \times \mathbf{u}_e - e\mathbf{B}/m_e c$.

$$ \frac{c}{4\pi e} \nabla \times \left( \frac{\nabla \times \mathbf{B}}{n} \right) + \frac{e\mathbf{B}}{m_e c} = \nabla \times \mathbf{u}_i - \mathbf{W}. \quad (60) $$

Some of the hybrid-kinetic codes make approximations of electron inertial terms in Eq. (57) and (60) to simplify their numerical solutions (Lipatov 2002; Shay et al. 1998; Kuznetsova et al. 1998). Spatial density variations are neglected in Eq. (60) (Shay et al. 1998; Kuznetsova et al. 1998). The electric field was then calculated from the generalized Ohm's law by neglecting the electron inertial term with time derivatives of the electron fluid velocity (Kuznetsova et al. 1998). These approximations are valid when length scale of variations is much larger than the electron inertial length. For a detailed discussion of these approximations, see Muñoz et al. (2018). These hybrid-kinetic codes which partially included electron inertial effects have mainly been used to study collisionless magnetic reconnection (Shay et al. 1999; Kuznetsova 2000; Kuznetsova et al. 2001). In particular Shay et al. (1998) used an evolution equation for a scalar electron pressure while Kuznetsova et al. (1998) included the full electron pressure tensor to take into account the non-gyrotropic effects.

More recently a hybrid-kinetic code CHIEF (Code Hybrid With Inertial Electron Fluid) was developed (Muñoz et al. 2018). This code solves Eqs. (57) and (60) without making any of the electron inertia related approximations used by other codes. The details of the numerical algorithm to solve Eqs. (51)-(58) are discussed by Muñoz et al. (2018). CHIEF was used to simulate kinetic plasma turbulence and it was found that the electron inertia related approximations are not valid in electron scale current sheets formed in the turbulence (Jain et al. 2022; Muñoz et al. 2023).

Some hybrid-kinetic codes with electron inertia calculate electric field from an elliptic PDE instead of Eq. (57) (Amano et al. 2014; Valentini et al. 2007). The elliptic PDE for the
electric field is obtained by taking curl of Faraday's law, Eq. (53), and utilizing Eqs. (54) and (57). Electron inertia effects were considered in the elliptic equation for the electric field while, still, the electron inertia term was ignored that contains the divergence of the electric field. Two dimensional simulations of kinetic plasma turbulence have shown that this approximation is not valid from ion to electron scales (Jain et al. 2022).

### 4.4.2 Hybrid-PIC Code CHIEF for Magnetic Reconnection Studies

Hybrid-PIC code CHIEF can be used to study magnetic reconnection with a large guide magnetic field in which case electron non-gyrotropy is weak and reconnection electric field is expected to be balanced by the electron inertial terms in generalized Ohm's law. Indeed, PIC simulations of magnetic reconnection have shown that electron inertial terms are significant to balance the reconnection electric field when guide field is large (Hesse and Winske 1998; Hesse et al. 2016; Liu et al. 2013; Pritchett 2005).

A particular reconnection study for which CHIEF code can be employed is electron-only reconnection in which ions do not couple to electrons (Sharma Pyakurel et al. 2019), examples of which are discussed in Sect. 3.7.2. Electron-only reconnection has recently been discovered in space observations (Phan et al. 2018). In many of electron only reconnection events observed by MMS, guide magnetic field is significantly larger than the asymptotic value of the reconnecting component magnetic field (Phan et al. 2018; Zhou et al. 2021; Man et al. 2020; Stawarz et al. 2022). In the statistical survey of electron only reconnection events in Earth's magnetosheath, guide field in majority of the events is 1 to 10 times larger than the reconnecting component of the magnetic field (Stawarz et al. 2022). MMS observations of electron scale reconnection in Earth's magnetosheath show that the electron non-gyrotropy, which balances the reconnection electric field in weak or zero guide field case, reduces with the increasing strength of the guide magnetic field (Wilder et al. 2018). No evidence of agyrotropy was found in another MMS observations of large guide field reconnection (Eriksson et al. 2016). For the electron only events observed with large guide magnetic field and/or absence of non-gyrotropy, electrons can be modeled as an inertial electron fluid. Hybrid-PIC code CHIEF treats electrons as inertial fluid without any approximations and can be used to study electron-only reconnection.

### 4.4.3 Outlook

Hybrid-kinetic simulations with electron inertia provide a computationally less expensive (in comparison to fully kinetic simulations) tool to study magnetic reconnection with guide field in which bulk electron inertia is the dominant mechanism breaking the frozen-in condition of magnetic field. There have been some hybrid-kinetic simulations studies with electron inertia of guide field magnetic reconnection (Kuznetsova et al. 1998; Kuznetsova 2000; Shay et al. 1998; Califano et al. 2020; Muñoz et al. 2023). More studies are, however, required to address still many open questions about the guide field magnetic reconnection. The EMHD limit (stationary ions) of the hybrid-kinetic model with electron inertia is particularly useful to study the nature of reconnection in electron scale current sheets and is relevant for the recently discovered electron-only reconnection (Phan et al. 2018). At the same time, simulations of kinetic plasma turbulence from ion to electron scales using hybrid-kinetic model with electron inertia will shed light on the conditions under which electron scale current sheets form and reconnect with or without ion coupling in the turbulence.
# 5 Fully Kinetic Particle-in-Cell Simulations

## 5.1 Introduction

The particle-in-cell (PIC) method is a simulation method in which the plasma is treated as a collection of particles (electrons and ions), where each species is typically composed of up to $10^{12}$ particles in 3D cases. In PIC simulations, the motions of individual particles and the evolution of electric and magnetic fields are solved self-consistently. The electric and magnetic fields are defined on discrete grid points. In this subsection, we will explain an explicit PIC simulation, where all the quantities are updated based on the quantities obtained in the previous time step. Let us assume that the total particle number in a simulation is $N_p$ (in other words, $N_p/2$ for ions, and $N_p/2$ for electrons). The equation of motion for the j-th particle's position $\mathbf{x}_j(t)$ and momentum $\mathbf{p}_j(t)$, where j represents an integer between 1 to $N_p$, is discretized in time, while Maxwell's equations for electromagnetic fields $E(\mathbf{x}, t)$ and $B(\mathbf{x}, t)$ are discretized in both space and time, using a grid spacing $\Delta x$ (assuming that the grids are uniform in all the coordinates, i.e. $\Delta x = \Delta y = \Delta z$) and a time step $\Delta t$, respectively. They are given as

$$ \frac{\mathbf{x}_j(t_n) - \mathbf{x}_j(t_{n-1})}{\Delta t} = \frac{\mathbf{p}_j(t_{n-1/2})}{m_j \gamma_j(t_{n-1/2})}, \qquad (61) $$

$$ \frac{\mathbf{p}_j(t_{n+1/2}) - \mathbf{p}_j(t_{n-1/2})}{\Delta t} = q_j \left[ \mathbf{E}(\mathbf{x}_j(t_n), t_n) + \frac{\mathbf{p}_j(t_n)}{c m_j \gamma_j(t_n)} \times \mathbf{B}(\mathbf{x}_j(t_n), t_n) \right], \qquad (62) $$

$$ \frac{\mathbf{E}(\mathbf{x}, t_{n+1}) - \mathbf{E}(\mathbf{x}, t_n)}{\Delta t} = -4\pi \mathbf{J}(\mathbf{x}, t_{n+1/2}) + c \nabla_f \times \mathbf{B}(\mathbf{x}, t_{n+1/2}), \qquad (63) $$

$$ \frac{\mathbf{B}(\mathbf{x}, t_{n+3/2}) - \mathbf{B}(\mathbf{x}, t_{n+1/2})}{\Delta t} = -c \nabla_f \times \mathbf{E}(\mathbf{x}, t_{n+1}), \qquad (64) $$

where $m_j$ is a mass, $\gamma_j$ is the Lorentz factor, $q_j$ is a charge, $c$ is the speed of light, and the operator $\nabla_f \times$ represents the finite difference version of the curl operation. The time is discretized to be $t_a = a\Delta t$, where a represents an integer $n$ or a half integer $n + 1/2$. Field quantities such as $\mathbf{B}(\mathbf{x}, t_n)$ are defined at a grid position $\mathbf{x}$, while quantities such as $\mathbf{B}(\mathbf{x}_j(t_n), t_n)$ are defined at the position of a particle $\mathbf{x}_j(t_n)$. Note that $\mathbf{B}(\mathbf{x}_j(t_n), t_n)$ in the right-hand side of Eq. (62) represents the mean of $\mathbf{B}(\mathbf{x}_j(t_n), t_{n+1/2})$ and $\mathbf{B}(\mathbf{x}_j(t_n), t_{n-1/2})$. As seen in Eqs. (61) and (62), the position $\mathbf{x}_j$ is computed at integer time, $t = t_n$, while the momentum $\mathbf{p}_j$ is computed at half-integer time, $t = t_{n+1/2}$. This time staggering gives the second-order accuracy, i.e. the error is $O(\Delta t^2)$. In the same way, for the spatial discretization for $E(\mathbf{x}, t)$ and $B(\mathbf{x}, t)$, the Yee lattice (Yee 1966) is used, in which each component of electric and magnetic fields is defined as in Fig. 6. Also, each component of the current density $\mathbf{J}(\mathbf{x}, t)$ is defined at the same position as $E(\mathbf{x}, t)$. These fields, which are spatially and temporarily staggered, are advanced using Eqs. (63) and (64), keeping the second-order accuracy in space and time. Equation (64) guarantees that the Gauss's law for the magnetic field, $\nabla_f \times \mathbf{B} = 0$, where $\nabla_f$ represents the finite difference version of the divergence, is satisfied when it is satisfied at $t = 0$.

To solve Eq. (62), the Boris method (Boris 1970) is commonly used. This method has three steps: (1) The momentum is updated from $\mathbf{p}_j(t_{n-1/2})$ to $\mathbf{p}_j(t_n)^*$, using only the electric field $E(\mathbf{x}_j(t_n), t_n)$ for a half time step $\Delta t/2$. (2) The momentum vector $\mathbf{p}_j(t_n)^*$ is rotated to be $\mathbf{p}_j(t_n)^{**}$ using only the magnetic field $\mathbf{B}(\mathbf{x}_j(t_n), t_n)$ for a full time step $\Delta t$. (3) The
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Yee lattice and electric and magnetic fields in a cell in a 3D case, where the length of each side is $\Delta x$. Electric fields $E$ are defined at the midpoint of each side of the cube, while magnetic fields $B$ are defined at the center of each face of the cube. In a 2D case, all the quantities are defined in the x-y plane, projecting each position onto the cell in the x-y plane</figcaption></figure>

rotated momentum is further updated from $p_j(t_n)^{**}$ to $p_j(t_{n+1/2})$, using the electric field $E(\mathbf{x}_j(t_n), t_n)$ for another half time step $\Delta t/2$.

In the PIC method, each particle is not a point particle, but it has a finite size to reduce noise. The shape of a particle depends on simulation codes, but the most commonly used shape function is a triangular function, $S_x(x-x_j) = (1-|x-x_j|/\Delta x)$ when $|x-x_j|/\Delta x < 1$ and zero otherwise, where only the $x$ component is considered. In 2D and 3D simulations, the $y$ and $z$ components of the shape functions are multiplied, as $S(\mathbf{x}-\mathbf{x}_j) = S_x(x-x_j)S_y(y-y_j)S_z(z-z_j)$ for 3D. Using these shape functions, the charge density $\rho(\mathbf{x}, t)$ is computed as $\rho(\mathbf{x}, t) = \sum_j q_j S(\mathbf{x}-\mathbf{x}_j)$. This way of charge assignment is reversed to compute the electric field exerted from each grid point to a particle's position. To avoid the self-force (the force due to the electric field generated by the particle itself), we must first average the electric fields defined on half-integer grids to obtain the mean electric field on each integer grid, before assigning the electric fields to the particle. The magnetic fields are assigned from grids to the particle position in the same way.

The current density can also be calculated using $J(\mathbf{x}, t) = \sum_j q_j v_j S(\mathbf{x} - \mathbf{x}_j)$, where $v_j$ is the velocity, but the calculated $J(\mathbf{x}, t)$ does not satisfy the continuum equation, $[\rho(\mathbf{x}, t_{n+1}) - \rho(\mathbf{x}, t_n)]/\Delta t + \nabla_f \cdot J(\mathbf{x}, t_{n+1/2}) = 0$; therefore, the electric field calculated using Eq. (63) with this $J(\mathbf{x}, t)$ does not satisfy the Gauss's law, $\nabla_f \cdot E(\mathbf{x}, t_{n+1}) = 4\pi\rho(\mathbf{x}, t_{n+1})$. This means that we must either correct the electric field $E(\mathbf{x}, t_{n+1})$ to satisfy the Gauss's law, or use another method to compute $J(\mathbf{x}, t)$. A technique for the former is explained in Birdsall and Langdon (1991). For the latter, for example, Villasenor and Bunemann (1992) developed a rigorous charge conservation method for 2D and 3D PIC simulations, which guarantees that both the Gauss's law and the continuum equation are satisfied at the same time, when they are satisfied at $t=0$.

The time step $\Delta t$ must satisfy the Courant–Friedrichs–Lewy condition, $\Delta t < \Delta x/(cN_d^{1/2})$, where $N_d$ represents the dimensionality ($N_d = 1, 2, \text{ or } 3$). Also, the grid spacing $\Delta x$ should be close to the Debye length $\lambda_D$, otherwise strong numerical heating occurs. Even when those conditions are satisfied, if particles are relativistic, a numerical Cherenkov
instability can occur and the noise field becomes extremely large. When this occurs, a noise reduction method such as by Godfrey (1980) can be used.

To study magnetic reconnection, it is important to separate the spatiotemporal scale of protons and that of electrons, which is controlled by the ion to electron mass ratio, $m_i/m_e$. In general, to reduce the computing time, full PIC simulations use an artificial mass ratio, such as $m_i/m_e = 25$ and 100, which is smaller than the real mass ratio $m_i/m_e = 1836$. There are two effects if we use a smaller mass ratio. One is that the thickness of the electron-scale current layer near the reconnection X-line, which is of the order of the electron skin depth $d_e = c/\omega_{pe}$, becomes thicker than that for a realistic case. The other is that the time scale of the electron physics becomes longer than the realistic case. For example, if we choose the mass ratio 25, the spatial scale separation between the ion scale ($d_i$) and the electron scale ($d_e$) is 5 times, and the temporal scale separation between the ion scale ($\Omega_i^{-1}$) and the electron scale ($\Omega_e^{-1}$) is 25 times. Even though these are much smaller than those in the realistic case, the choice of the mass ratio 25 is acceptable in order to see the scale separation physics that occurs in a study of magnetic reconnection.

Various boundary conditions can be implemented including periodic, conducting wall, and open boundaries (Daughton and Scudder 2006; Ohtani and Horiuchi 2009). For open boundaries, particles that reach the boundaries are removed, and new particles are injected into the simulation box at each time. How to inject particles depends on simulation codes. For example, in Daughton and Scudder (2006), the flux of injecting particles is calculated assuming that the spatial derivative of the distribution function at the boundary is zero in the normal direction. Also, electromagnetic fluctuations can pass through the boundaries.

Some research groups use PIC codes with the adaptive mech refinement (AMR) technique (Fujimoto 2011; Innocenti et al. 2013). In AMR PIC simulations, the simulation region is subdivided based on the required spatial resolution: the regions where the small-scale physics becomes important are solved using fine grids, while the other regions outside the fine regions are solved using coarse grids. If more smaller-scale resolution is required, further finer levels of grids are produced. This technique can reduce the required particle number for fully kinetic simulations.

## 5.2 Magnetotail Reconnection

In this section, we introduce simulation results of magnetic reconnection obtained by the standard fully kinetic PIC simulations with uniform grids. In the Earth's magnetotail, the strength of magnetic field across the current sheet is symmetric, and the guide field ($B_y$ field in the GSM coordinates) is small in many reconnection events. Many authors have been studying symmetric magnetic reconnection with zero guide-field, using the Earth's magnetotail parameters (Hoshino 1987; Pritchett et al. 1991; Horiuchi and Sato 1994; Dreher et al. 1996; Zhu and Winglee 1996; Hesse et al. 1996). In these simulations, the initial plasma is set up based on a Harris equilibrium (Harris 1962a): the magnetic field $B_x = B_0 \tanh(z/w)$, and the density $n = n_0 \text{sech}^2(z/w) + n_b$, where $B_0$ is the asymptotic magnetic field, $w$ is the sheet thickness, $n_0$ is the peak density of the current sheet, and $n_b$ is the background density. Also the conditions for the Harris equilibrium, $B_0^2/(8\pi) = n_0(T_i + T_e)$ and $|V_{di} - V_{de}| = [2c/(weB_0)](T_i + T_e)$, and $V_{di}/V_{de} = -T_i/T_e$, are satisfied, where $T_i$ and $T_e$ are the ion and electron temperatures, respectively, $e$ is the elementary charge, and $V_{di}$ and $V_{de}$ are the y-directional drift velocity in the current sheet component of ions and electrons, respectively.

The following describes an example of a 2D PIC simulation of magnetotail reconnection, Hesse et al. (2018a). The system size is $L_x \times L_z = 102.4d_i \times 51.2d_i$, where
<figure><img src="image_9.png" /><figcaption>**Fig. 7** 2D PIC simulation result for magnetotail reconnection. (a) Electron fluid velocity $V_{\text{ex}}$, (b) electron temperature $T_e$, (c) zoom-in view of $T_e$, and (d) the heating term in Eq. (65). Adapted from Hesse et al. (2018a, 2019)</figcaption></figure>

$d_i$ is the ion skin depth, $c/\omega_{pi}$ with $\omega_{pi}$ being the ion plasma frequency based on $n_0$ ($\omega_{pi} = (4\pi n_0 e^2/m_i)^{1/2}$), and 3200 × 3200 grids are used. The mass ratio is $m_i/m_e = 100$, the sheet thickness is $w = 0.5d_i$, the temperature ratio is $T_i/T_e = 5$, the density ratio $n_b/n_0 = 0.2$, and the ratio of the plasma frequency (based on $n_0$) to the electron cyclotron frequency (based on $B_0$) is $\omega_{pe}/\Omega_e = (4\pi n_0 e^2/m_e)^{1/2}/[eB_0/(m_ec)] = 2.0$, which gives the ratio of the light speed to the Alfvén speed (based on $B_0$ and $n_0$, $v_{A0} = B_0/(4\pi m_i n_0)^{1/2}$) to be $c/v_{A0} = 20.0$. The x boundaries are periodic, and the z boundaries are conducting walls. To initiate magnetic reconnection, a perturbation is added to the magnetic field as $\delta B_x = (a_0\pi/L_z)\cos(2\pi x/L_x)\sin(\pi z/L_z)$ and $\delta B_z = -(a_02\pi/L_x)\sin(2\pi x/L_x)\cos(\pi z/L_z)$, which gives a reconnection X-line at the origin $x=0$ and $z=0$. Here $a_0$ is the amplitude of the perturbation. The total number of particles used in the simulation is $7 \times 10^{10}$.

After the simulation starts, the reconnection electric field $E_y$ is generated in the diffusion region near the X-line, where both the ion and electron motions are decoupled from the magnetic field line motion, which allows a pair of magnetic field lines across the current sheet, one is in the positive z region ($B_x > 0$) and the other is in the negative z direction ($B_x < 0$), are reconnected and energy conversion occurs from the magnetic energy to the kinetic and thermal energies of ions and electrons. The reconnection rate is measured as $E_y/(B_0v_{A0}/c)$, and in this simulation it is around 0.2 (Hesse et al. 2018b). Both ion and electron outflows are produced from the X-line toward the positive and negative x directions. Figure 7(a) shows the electron fluid velocity $V_{\text{ex}}$, where the bipolar positive and negative $V_{\text{ex}}$ peaks appear along the $z=0$ in $40 < x/d_i < 60$, and each peak value reaches near the electron Alfvén speed, $v_{\text{Ae}} = B_0/(4\pi m_e n_0)^{1/2} = (m_i/m_e)^{1/2}v_{\text{A0}}$. Outside the region of $40 < x/d_i < 60$, there are strong inflows, which also reach near $v_{\text{Ae}}$, toward the X-line along the separatrices. Because of these strong counter streaming electron flows (outflows and inflows), an electrostatic instability occurs that produces waves propagating along the separatrices toward the X-line, and electrons are heated due to wave-particle interactions.
Figure 7(b) and (c) show the electron temperature in the entire box and a zoom-in view that includes a separatrix. Electrons are heated inside the separatrices (panel (b)), and the zoom-in view (panel (c)) shows that there are two solitary structures due to the nonlinear evolution of the wave (at $x = 64d_i$ and $66d_i$ along the separatrix) where electron temperature significantly enhances.

The locations of the instability, along the separatrices, correspond to the boundary of the high electron temperature, which suggests the importance of the electrostatic instability to heat electrons. To understand the effect of the instability on the heating, Hesse et al. (2018a) analyzed the pressure equation:

$$ \frac{\partial p}{\partial t} = -\nabla \cdot (V p) - \frac{2}{3} \sum_{l} P_{ll} \frac{\partial}{\partial x_l} V_l - \frac{1}{3} \sum_{l,i} \frac{\partial}{\partial x_i} Q_{lii} - \frac{2}{3} \sum_{l,i(l \neq i)} P_{li} \frac{\partial}{\partial x_i} V_l, \quad (65) $$

where all the quantities are for electrons (the subscript $e$ is omitted): $p$ is the scalar pressure, $V$ is the fluid velocity, and $P_{ij}$ and $Q_{ijk}$ are the pressure tensor and the heat tensor, respectively. The first two terms represent the compression effect, the third term is due to the heat flux, and the fourth term represents the quasi-viscous effect due to the off-diagonal components of the pressure tensor. Figure 7(d) shows the contribution of each term in Eq. (65), integrated over the region of $A > A_0$, where $A$ is the $y$ component of the magnetic flux function (i.e., $B_x = \partial A/\partial z$ and $B_z = -\partial A/\partial x$). Note that $A = 0$ at the outermost $z$ boundaries, and $A$ is decreasing as we approach $z = 0$. Figure 7(d) indicates that the quasi-viscous term (the fourth term) is the dominant term to provide the pressure increase, leading heating in the reconnection region.

## 5.3 Magnetopause Reconnection

PIC simulations of magnetopause reconnection can include certain challenges due to asymmetries in the densities, temperatures and magnetic field strengths of the abutting plasmas (Sonnerup et al. 1986; Cassak and Shay 2007). Specifically, magnetospheric plasma is usually relatively sparse, hot, and threaded by a strong magnetic field, while magnetosheath plasma (which arises from shocked solar wind) is denser, cooler, and includes a somewhat weaker field.

From a simulation perspective, the density asymmetry – which can exceed an order of magnitude – can be particularly problematic. PIC simulations are inherently noisy. The random fluctuations tend to follow Poissonian statistics with an amplitude scaling as $1/\sqrt{N_{pc}}$, with $N_{pc}$ the number of (macro) particles per computational cell. If variations in the number of macroparticles directly translate to variations in density, a 16 : 1 ratio between the magnetosheath and magnetospheric plasma densities will produce noise levels an unacceptable four times larger in the latter than the former. One obvious approach – throwing more particles at the problem – can quickly become computationally burdensome. An alternative is the use of particle weighting, in which each particle is assigned a weight that determines its significance in the calculation of particle moments (e.g., charge and current density). Doing so allows for an initially uniform distribution of particles with a roughly constant noise level. More sophisticated algorithms allow for the splitting and joining of particles as the simulation progresses to account for the development of density variations and to address computational load imbalances.

A 2D simulation of the magnetopause with p3d, a PIC code employing weighted particles (Zeiler et al. 2002), was presented in Swisdak et al. (2018). In its normalization, a reference magnetic field strength $B_0$ and density $n_0$ define the velocity unit
$v_{A0} = B_0/(4\pi m_i n_0)^{1/2}$. Times are normalized to the inverse ion cyclotron frequency $\Omega_{i0}^{-1} = m_i c/(eB_0)$, lengths to the ion inertial length $d_{i0} = c/\omega_{pi0}$ (where $\omega_{pi0} = (4\pi n_0 e^2/m_i)^{1/2}$ is the ion plasma frequency), electric fields to $v_{A0}B_0/c$, and temperatures to $m_i v_{A0}^2$.

The initial conditions closely mimic those observed during the diffusion region encounter described in Burch et al. (2016). In the system considered here, $B_0$ and $n_0$ correspond to their asymptotic magnetosheath values: $B_0 = 23$ nT and $n_0 = 11.3$ cm$^{-3}$. The simulation uses an LMN coordinate system in which the reconnecting field parallels the L axis (roughly north-south), the M axis runs roughly east-west, with dawnward positive, and the N axis points radially away from Earth and completes the right-handed triad. The computational domain has dimensions ($L_L, L_N$) = (40.96, 20.48) with periodic boundary conditions used in all directions. While particles can move in the M direction, variations in physical quantities are not permitted: $\partial/\partial M = 0$.

The reconnecting component of the field $B_L$ and the ion and electron temperatures, $T_i$ and $T_e$, vary as functions of N with hyperbolic tangent profiles of width 1. The asymptotic values of $n, B_L, T_i$, and $T_e$ in code units are 1.0, 1.0, 1.37, and 0.12 in the magnetosheath and 0.06, 1.70, 7.73, and 1.28 in the magnetosphere. Pressure balance determines the initial density profile. The guide field $B_M = 0.099$ is much smaller than $B_L$ (i.e., the reconnection is nearly anti-parallel) and initially uniform. While not an exact kinetic equilibrium, the unperturbed configuration is in force balance and would not undergo significant evolution during the timescales of interest. Instead, a small initial perturbation is introduced to trigger reconnection onset.

The ion-to-electron mass ratio is chosen to be 100, which is sufficient to separate the electron and ion scales (the electron inertial length $d_{e0} = 0.1 d_{i0}$). The normalized speed of light is $c = 15$ so that $\omega_{pe}/\Omega_e = 1.5$ in the asymptotic magnetosheath and $\approx 0.2$ in the asymptotic magnetosphere; the observed ratios are larger, $\approx 46$ and 7, and as a consequence the simulation's Debye length is larger than in the real system. However, since the development of reconnection does not appreciably depend on physical effects at the Debye scale the expected impact is minimal. The spatial grid has resolution $\Delta = 0.01$ in normalized units while the Debye length in the simulation's magnetosheath, $\approx 0.03$, is the smallest physical scale. To ameliorate numerical noise, particularly in the low-density magnetosphere, each grid cell initially contains 3000 weighted macroparticles.

Figure 8 shows results. The asymmetry in the field strength is apparent in the distribution of the field lines, with the separatrices extending much farther (in the N direction) into the plasma of the magnetosheath (top) than the magnetosphere (bottom). The Hall electric and magnetic fields (panels c and e) differ significantly from the case of symmetric reconnection, with the former concentrated almost exclusively on the magnetospheric side while the latter is almost completely dipolar rather than quadrupolar. Due to the use of weighted particles, the numerical noise is similar on both sides.

## 5.4 The 3D Nature of Magnetic Reconnection

The third dimension pointing out of the 2D reconnection plane introduces numerous plasma instabilities (e.g. Daughton et al. 2011; Graham et al. 2025). In addition to these secondary instabilities, inherent 3D nature of reconnection X-line is also omitted in 2D pictures, which includes the effect of limited X-line extent (Shay et al. 2003; Liu et al. 2019; Huang et al. 2020; Pyakurel et al. 2021; Huang et al. 2024), its tendency of spreading (Huba 2003; Shay et al. 2003; Lapenta et al. 2006; Nakamura et al. 2012; Shepherd and Cassak 2012; Li et al. 2020; Arencibia et al. 2023; Li et al. 2023c; Lin et al. 2025), and its orientation preference (Sonnerup 1974; Swisdak and Drake 2007; Hesse et al. 2013; Aunai et al. 2016; Liu et al. 2013, 2018).
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Simulation results overplotted with magnetic field lines. (a) The **J**·**E** term from Poynting's theorem; (b) In-plane electron flow field; (c) $E_N$, the normal component of the electric field; (d) $E_\|$, the component of the electric field parallel to the magnetic field; (e) $B_M - B_{M,0}$, the change in the out-of-plane component of the magnetic field from its (spatially constant) initial value; (f) $S_L$, the horizontal component of the Poynting flux. Adapted from Swisdak et al. (2018)</figcaption></figure>

We briefly highlight the property of reconnection X-lines with a short extent here. To sustain a current sheet, electrons and ions drift in opposite directions. This fact introduces the asymmetry along the X-line (current) direction. To reveal this effect, Liu et al. (2019) and Huang et al. (2020) studied magnetic reconnection with the X-line spatially confined in the current direction. They included thick current layers to prevent the reconnection from spreading out of the two ends of a thin current sheet that has a thickness on an ion inertial ($d_i$) scale. The $x$ component of the magnetic field is given as $B_x = B_0 \tanh[z/L(y)]$, where the half-thickness $L(y) = L_{min} + (L_{max} - L_{min})[1 - f(y)]$ and $f(y) = [\tanh((y+w_0)/S) - \tanh((y-w_0)/S)]/[2\tanh(w_0/S)]$, $L_{min} = 0.5d_i$, $L_{max} = 4d_i$, and $S = 5d_i$. The parameter $w_0$, which controls the $y$-extent of the thin current region ($L_{y-thin}$), is varied from $w_0 = 2d_i$ to $20d_i$, corresponding to the $y$-extent of the thin current region from $L_{y-thin} \sim 4d_i$ ($w_0 = 2d_i$) to $L_{y-thin} \sim 30d_i$ ($w_0 = 20d_i$). The density is $n = n_0 \text{sech}^2[z/L(y)] + n_b$, and $n_b = 0.3n_0$. The size of the system is $L_x \times L_y \times L_z = 32d_i \times 64d_i \times 16d_i$. The periodic boundary condition is used for $x$ and $y$, and the conducting walls are placed in the $z$ boundaries. Over $2.6 \times 10^{10}$ particles for each species are used.

The resulting reconnection is shown in Fig. 9, which is for $L_{y-thin} \sim 30d_i$. Liu et al. (2019) found that the reconnection rate and the outflow speed drop significantly when the extent of the thin current sheet, $L_{y-thin}$, is less than $O(10d_i)$. When the thin current sheet extent is long enough, it consists of two distinct regions: a suppressed reconnecting region (on the ion-drifting side, marked in Fig. 9(b)) exists adjacent to the active region where reconnection proceeds normally as in a 2D case with a typical fast rate value $\approx 0.1$. The extent of this suppression region is $O(10d_i)$, and it suppresses reconnection when $L_{y-thin}$ is comparable or shorter. The time scale of current sheet thinning toward fast reconnection can be translated into the spatial scale of this suppression region, because the electron drifts
<figure><img src="image_11.png" /><figcaption>**Fig. 9** 3D PIC simulation results. (a) 3D view of reconnection with a limited X-line extent, where the thin current sheet region extends $30d_i$ in y (Huang et al. 2020). The mass ratio in the simulation is 25. (b) The current density on the $x=0$ plane (left) and magnetic field $B_z$ on the $z=0$ plane (middle and right) (Liu et al. 2019). The mass ratio is 75. The gray shaded area represents the “suppressed reconnecting region”. Adapted from Huang et al. (2020) and Liu et al. (2019)</figcaption></figure>

inside the ion diffusion region transport the reconnected magnetic flux (that is critical in driving outflows and furthers the current sheet thinning) away from this region. This is a consequence of the Hall effect in 3D.

Huang et al. (2020) incorporated the length scale of this suppression region $O(10d_i)$ to quantitatively model the reduction of the reconnection rate and the maximum outflow speed observed in the short X-line limit. The average reconnection rate drops because of the limited active region (where the current sheet thins down to the electron inertial scale) within the X-line. The outflow speed reduction correlates with the decrease of the $J \times B$ force, which can be modeled by the phase shift between the J and B profiles, also as a consequence of the flux transport out of the reconnection plane.

While the existence of this suppression region may explain the shortest possible azimuthal extent of dipolarizing flux bundles at Earth (Liu et al. 2015), it may also explain the dawn-dusk asymmetry observed at the magnetotail of Mercury (Sun et al. 2016, 2022), which has a global dawn-dusk extent much shorter than that of Earth.

## 5.5 Particle Acceleration

There have been quite remarkable advances in using PIC simulations to understand particle acceleration processes in magnetic reconnection, discussed in Oka et al. (2023), Drake et al. (2025), and Guo et al. (2023) of this collection. The simulation provided energetic particle flux, spectra and even detailed distributions that can be compared with in situ ob-
<figure><img src="image_12.png" /><figcaption>Fig. 10 An example of particle trajectory analysis. Adapted from Oka et al. (2010)</figcaption></figure>

servations. We introduce several key diagnostics recently used for gaining insight in particle energization.

First, it has been a common practice to output particle trajectories to study the acceleration process (e.g., Hoshino et al. 2001; Drake et al. 2006; Fu et al. 2006; Oka et al. 2010; Guo et al. 2015). These have led to the identification of different acceleration mechanisms, as discussed in Oka et al. (2023) and Guo et al. (2023) of this collection. Figure 10 shows a representative particle trajectory adapted from Oka et al. (2010). This particle is first accelerated by an X-line (a), then further energized due to electric field during anti-reconnection between two merging island (c). The acceleration persists after the particle is ejected out of the X-line region. In addition, one can output the electric and magnetic fields and other quantities associated with particles, to complement the understanding of acceleration mechanisms.

The limitation of just showing several particle trajectories, even with the best effort, is that the “representative” examples are usually cherry-picking results and it is difficult to evaluate the relative importance of each mechanism. There has been recent effort to evaluate the acceleration processes over a large number of trajectories, (Guo et al. 2019, 2021a; Kilian et al. 2020; French et al. 2022; Li et al. 2023a).

Another method, developed and widely used over the last decade, is to study the collective energy gain, such as guiding center and pressure-restrained terms (discussed in Oka et al. 2023, this collection) using ensemble averaged moments (Dahlin et al. 2014, 2015; Li et al. 2015, 2017, 2018, 2019b; Du et al. 2018). For example, the acceleration due to the curvature (Fermi) and gradient (betatron) drifts can be evaluated under guiding center approximation. Figure 11 shows an example under the guiding-center approximation, and shows the curvature drift term is the main acceleration term. Moreover, it is possible to collect the energy dependent acceleration rates by considering particles with different energy,
<figure><img src="image_13.png" /><figcaption>**Fig. 11** (a) An example of guiding-center drift analysis. Particle energization due to different drift currents for electrons (top) and ions (bottom). $j_c$ is due to particle curvature drift. $j_g$ is due to particle grad-$B$ drift. $j_m$ is due to magnetization. $j'' = j_c + j_g + j_m$. $\dot{K}_e$ and $\dot{K}_i$ are the energy change rates for electrons and ions, respectively. They are all normalized by $m_e c^2 \omega_{pe}$. (b) Similar to (a) but shows energy dependent values (Li et al. 2019b) at a given time. Adapted from Li et al. (2017)</figcaption></figure>

so the energization can be studied in a energy-dependent fashion (Dahlin et al. 2017; Guo et al. 2014; Li et al. 2018, 2019b). Figure 11b shows such an example.

## 5.6 Simulations of Magnetic Reconnection in Shock Waves

Magnetic reconnection can occur in current sheets generated in plasma turbulence, and PIC simulations have also been applied to turbulent environments (Wu et al. 2013; Matthaeus et al. 2016; Haggerty et al. 2017; Shay et al. 2018; Vega et al. 2020; Adhikari et al. 2021; Rueda et al. 2021), including turbulence in Kelvin Helmholtz vortices in the magnetopause flank region (Nakamura and Daughton 2014; Nakamura et al. 2017, 2022), and the transition region in shock waves (Matsumoto et al. 2015; Bohdan et al. 2017, 2020; Bessho et al. 2019, 2020, 2022, 2023; Ng et al. 2022).

Here, let us review 2D and 3D PIC simulation studies of magnetic reconnection in the shock turbulence. Bessho et al. (2019) used a 2D domain to study a quasi-parallel shock under the parameters in the Earth's bow shock. The size of the simulation domain is $L_x \times L_y = 375d_i \times 51.2d_i$, where the ion skin depth $d_i$ has 40 grids. The plasma is uniform at $t = 0$, both ions and electrons are Maxwellian with their temperatures $T_i$ and $T_e$, respectively, and the magnetic field is given as $B = [B_0 \cos\theta, B_0 \sin\theta, 0]$, where $\theta$ is the shock angle with respect to the $x$ axis. Periodic boundaries are used in the $y$ direction, and conducting walls are placed in the $x$ direction. To all the plasma particles, a negative drift speed, $-v_d$, in the $x$ direction is given, and a uniform positive $z$ component of electric field, as $E_z = v_d B_0 \sin\theta/c$, is set in the domain. At the right boundary, $x = L_x$, new particles for both ions and electrons are injected, using the same temperatures as the initially loaded particles, with the negative drift speed $-v_d$. At the left boundary, $x = 0$, all the particles are specularly reflected, and the incident particles and the reflected particles generate counter-streaming beams, which cause a beam instability. As a result, a non-linear wave grows near the left boundary, and a wave steepening occurs. Eventually, a shock wave forms, propagating toward the positive $x$ direction.

In the 2D simulation, the following parameters are used: the electron and ion beta $\beta_e = \beta_i = 1$, the ratio of the plasma frequency to the electron cyclotron frequency $\omega_{pe}/\Omega_e = 4$, the shock angle $\theta = 25^\circ$, and the mass ratio $m_i/m_e = 200$. With these parameters, the
<figure><img src="image_14.png" /><figcaption>**Fig. 12** PIC simulations of reconnection in shocks. (a) 2D simulation domain and the current density $J_z$. (b) Electron fluid velocity $V_{ex}$. (c) Ion fluid velocity $V_{ix}$. (d) 3D simulation domain and the current density $J_z$. Adapted from Bessho et al. (2019) and Ng et al. (2022)</figcaption></figure>

electron thermal speed becomes $v_{Te} = 14.4v_A$. The drift speed is set to be $v_d = 9v_A$. In the simulation (the downstream rest frame), the shock speed is $2.4v_A$, which corresponds to the Alfvén Mach number of the shock wave $M_A = 11.4$. In other words, the shock speed in the laboratory frame is $11.4v_A$, which is less than $v_{Te}$, consistent with the Earth's bow shock.

In the simulation, the shock transition region shows a non-resonant ion-ion beam instability due to the interactions between the ions reflected by the shock and the incident ions, and many current sheets are generated, some of which show signatures of magnetic reconnection. In Fig. 12(a), the color shows the current density $J_z$ in the 2D simulation domain, where the black curves are magnetic field lines projected onto the x-y plane, and magenta X marks represent the positions of reconnection X-lines. One of the reconnecting current sheet is zoomed up in Fig. 12(b) and (c), where the electron fluid velocity $V_{ex}$ and the ion fluid velocity $V_{ix}$ are shown. There is one magnetic island above the current sheet, and there are bipolar electron jets generated from the X-line. In contrast, the ion velocity plot does not show ion jet structures, and the ions are passing through the reconnection region with a negative $V_{ix}$. Therefore, this region is a site of electron-only reconnection, where only electrons are participating in reconnection, while ions cannot respond to the strong gradient of magnetic fields in the thin current sheet, whose thickness is less than the ion skin depth $d_i$. Electron-only reconnection has been observed in the Earth's magnetosheath (Phan et al. 2018; Gingell et al. 2021; Stawarz et al. 2022) and the transition region of the Earth's bow shock (Wang et al. 2019; Gingell et al. 2019, 2020). Note that the shock transition region has a negative $B_z$ magnetic field, $B_z \sim -4B_0$ (not shown), and the reconnecting magnetic field is the same order. Therefore, in the 2D simulation, reconnection in the shock transition region is guide-field reconnection.

Ng et al. (2022) performed a 3D PIC simulation to study reconnection in the shock transition region. The simulation parameters are: $\beta_e = \beta_i = 1.41$, $\omega_{pe}/\Omega_e = 4$, $m_i/m_e = 100$, $\theta = 30^\circ$, $v_d = 10v_A$, and the system size $L_x \times L_y \times L_z = 200d_i \times 50d_i \times 20d_i$. The z direction is set to be a periodic boundary. Figure 12(d) shows the current density $J_z$. In the 3D simulation, the current direction can be not only in the z direction, but also in the y direction; therefore, the reconnection plane does not have to be in the x-y plane as in the
2D simulation, and some current sheets show reconnection with a weak guide field, even though the shock transition region has a large negative $B_z$.

# 6 Embedded PIC: MHD-AEPIC

## 6.1 Overview

Due to the large separation between the kinetic scales and the size of Earth's magnetosphere, it is highly computationally expensive to apply a purely kinetic code for simulating global magnetospheric dynamics. Various hybrid methods have been proposed to incorporate kinetic physics into global simulations while keeping the computational costs feasible. Traditional hybrid codes model the electron species as a fluid and simulate the ions with either macro-particles or a grid-based Vlasov solver. These hybrid models reduce the separation between the kinetic scales and the global scale by removing the electron kinetic scales from the model so that it becomes feasible to apply them to Earth's magnetosphere.

The Magnetohydrodynamic with Adaptively Embedded Particle-in-Cell (MHD-AEPIC) model represents another type of hybrid approach to incorporate kinetic effects into global models (Daldorff et al. 2014; Shou et al. 2021; Chen et al. 2023). This type of hybrid model couples a kinetic code with a global fluid model, and only applies the kinetic code to simulate part of the simulation domain, where kinetic physics is crucial while using the fluid model to simulate the rest of the domain. Compared to a purely kinetic model, this type of hybrid model reduces the computational cost by reducing the domain size for the kinetic code, and it is best suited for applications where the important kinetic physics is localized. The MHD-AEPIC model, and its precursor, the Magnetohydrodynamic with Embedded Particle-in-Cell (MHD-EPIC) model, are the first two-way coupled models that work for global applications. Since then, similar coupled models have been developed by different independent teams. For example, Makwana et al. (2017) also developed a model that couples a PIC code with an MHD code, and Rieke et al. (2015) tried to couple a Vlasov solver with a two-fluid code.

## 6.2 Methodology

### 6.2.1 Development History

The original MHD-EPIC model was developed by Daldorff et al. (2014), in which the semi-implicit particle-in-cell code iPIC3D (Markidis et al. 2010) is coupled with the global fluid model BATS-R-US (Powell et al. 1999) through the Space Weather Modeling Framework (SWMF) (Tóth et al. 2005). The model has been successfully applied to study magnetic reconnections in the magnetospheres of Ganymede (Tóth et al. 2016; Zhou et al. 2019, 2020), Earth (Chen et al. 2017, 2020; Wang et al. 2022a,b), Mercury (Chen et al. 2019) and Mars (Ma et al. 2018). A PIC region has to be a box in the MHD-EPIC model. To cover the kinetic regions of interest, the MHD-EPIC model supports applying multiple independent kinetic regions (Tóth et al. 2016) in the same simulation domain, and it also allows rotating a box so that the corresponding PIC region does not have to be aligned with the global grid (Chen et al. 2020). These two features expand the capabilities of the MHD-EPIC model. However, not all the kinetic regions of interest can be covered by one or a few boxes. If the kinetic region moves at the global spatial scale during a simulation, the PIC box has to be very large to cover the whole region of interest, which is computationally expensive.
<figure><img src="image_15.png" /><figcaption>**Fig. 13** A schematic shows the improvement of the MHD-AEPIC (right) model from the MHD-EPIC (left) model. Adapted from Chen et al. (2023)</figcaption></figure>

To overcome these difficulties, the MHD-AEPIC model has been developed, which allows a dynamic PIC region of any shape (see Fig. 13).

To support dynamic PIC regions, two different PIC codes, the Adaptive Mesh Particle Simulator (AMPS) (Shou et al. 2021) and the FFlexible Exascale Kinetic Simulator (FLEKS) (Chen et al. 2023), have been developed as the PIC component of the MHD-AEPIC model. Since FLEKS is more widely used for MHD-AEPIC simulations, we focus on the FLEKS code in this section.

## 6.2.2 Coupling Algorithm

MHD-AEPIC supports coupling with single-fluid MHD, multi-species MHD, multi-ion MHD (Glocer et al. 2009b), and five- and six-moment multi-fluid models (Huang et al. 2019). The most widely used fluid component is the single-fluid Hall MHD with a separate electron pressure equation, and we briefly describe the coupling algorithm for this case here.

In an MHD-AEPIC simulation, the PIC code covers part of the whole simulation domain. The MHD model provides the initial conditions for the PIC code at the beginning of the coupled simulation. Once the initialization is done, both the PIC and MHD models update independently for one or a few time steps until the next coupling time point is reached. During the coupling, the MHD model provides the boundary conditions for the PIC code, and the PIC code provides the updated magnetic field and plasma quantities to overwrite the overlapped MHD region. Since the MHD and PIC codes solve different sets of equations, conversion between the MHD and PIC variables is needed. When calculating PIC variables from MHD variables, we need densities, velocities, and pressures for both electron and ion species, and they are calculated as follows:

* Charge neutrality is assumed, so both the electron and ion densities can be easily obtained from total MHD density.

* From the MHD magnetic field, the current density can be calculated. Since the sum of electron and ion momentum is the total MHD momentum, and the velocity difference between electrons and ions produces the current, the electron and ion velocities can be obtained.

* Since we usually solve both ion and electron pressure equations on the MHD side, the electron and ion pressures can be obtained directly to initialize thermal PIC macro-particles.
Once the electron velocity is obtained, it is used to calculate the electric field **E** for PIC from the generalized Ohm's law:

$$ \mathbf{E} = -\frac{\mathbf{U}_e \times \mathbf{B}}{c}, \qquad (66) $$

where **B** is the MHD magnetic field and $\mathbf{U}_e$ is the electron bulk velocity including the Hall term. We note that no matter Hall physics is included or not into the MHD model, the equation above is applied to calculate the initial and boundary electric field for PIC.

Calculating MHD variables from PIC variables is more straightforward: we simply sum up the mass, momentum and energy of the electron and ion macro-particles to obtain the plasma variables required. We refer the readers to Daldorff et al. (2014) for more details. Currently, the PIC codes used for MHD-EPIC/MHD-AEPIC coupling have to use a Cartesian mesh, but the MHD model BATSU-R-US can use non-uniform Cartesian or non-Cartesian grids. The interpolation between the PIC and MHD grids is done by a second-order linear interpolation.

### 6.2.3 Particle-in-Cell Algorithm

The embedded PIC model is a particular version of the PIC models discussed in detail in Sect. 5. The original MHD-EPIC implementation used the iPIC3D model while MHD-AEPIC uses FLEKS. Both MHD-EPIC and FLEKS are semi-implicit (Brackbill and Lapenta 2008; Lapenta 2017; Chen and Tóth 2019), meaning that the electric field is solved for by an implicit scheme. We choose the semi-implicit PIC algorithm because it has a relaxed stability constraint so that the Debye length does not have to be resolved and the stability constraint for the time step is based on the thermal speed instead of the speed of light. Based on our numerical experiments, we found the stability of the PIC code is extremely important for a successful MHD-EPIC/MHD-AEPIC simulation. To improve the stability, we designed the Gauss's Law satisfying Energy-Conserving Semi-Implicit Method (GL-ECSIM) (Chen and Tóth 2019), which is based on the Energy-Conserving Semi-Implicit Method (ECSIM) by Lapenta (2017). GL-ECSIM shares the same energy conservation property as ECSIM, i.e., the total energy of the system can be exactly conserved with proper parameters. In practice, we found the code is more stable with parameters that slowly dissipate the total energy numerically. In addition, satisfying Gauss's law (charge conservation) is also crucial for the stability and accuracy of the PIC code. GL-ECSIM applies a novel method to satisfy Gauss's law by adjusting particle positions at the end of each cycle. The details of GL-ECSIM can be found in Chen and Tóth (2019).

In a long MHD-AEPIC simulation, the macro-particle number per cell may vary significantly due to the transport of particles. The uneven distribution of particle numbers can cause load imbalance and reduce computational efficiency. To alleviate this problem, we designed particle splitting and merging algorithms for FLEKS. A particle splitting (merging) algorithm is applied to split (merge) particles when the number of particles per cell is below (above) a threshold (Chen et al. 2023).

### 6.2.4 Kinetic Region Adaptation

The most important improvement of MHD-AEPIC over MHD-EPIC is the adaptive PIC region. Although the PIC grid is still Cartesian, its cells can be switched on or off so that the active cells can fit any shape of kinetic regions. We note that the PIC cells can be activated or deactivated dynamically during a simulation. The active PIC region can be defined either
based on geometric or physical criteria. For physics-based adaptation, BATS-R-US calculates the physical criteria and sends the corresponding grid information to FLEKS to turn on or turn off cells.

### 6.2.5 Kinetic Scaling

In some applications, the difference between kinetic and global spatial and temporal scales makes it difficult, if not impossible, to resolve the kinetic scales in an MHD-EPIC, or even MHD-AEPIC simulation. Fortunately, the large separation of scales can be exploited, and the kinetic scales can be increased by changing the mass per charge ratio without affecting the global dynamics (Tóth et al. 2017). This technique is not needed or even applicable for Ganymede and Mercury simulations, where the kinetic and global scales are not very different. On the other hand, kinetic scaling is applicable and extremely useful for modeling Earth's magnetosphere. We typically increase the kinetic scales by a factor of 4 to 16. See Tóth et al. (2017) for more detail.

## 6.3 Applications

The MHD-EPIC/MHD-AEPIC model has been applied to investigate the physical processes and consequences of both magnetopause and magnetotail reconnection. In these simulations, the PIC code is usually used to cover either the magnetopause or the magnetotail current sheet, where reconnection happens. Since the initial and boundary conditions of the PIC code are obtained from the MHD model, the physical parameters inside the PIC region, such as the plasma quantities and the shape of the current sheet, are more realistic than those in a standalone PIC simulation. On the other hand, the information from the PIC code is also fed back to the MHD model so that we can evaluate the global consequences of the kinetic magnetic reconnection.

Here we briefly describe a few applications of the MHD-EPIC/MHD-AEPIC model to study Earth's magnetosphere. Chen et al. (2017) studied both the kinetic features of magnetopause reconnection and the evolution of flux transfer events (FTEs) shown in Fig. 14. Near the reconnection site, the simulation successfully produced key kinetic features of asymmetric magnetic reconnection, such as the crescent electron phase space distribution and the lower hybrid drift instability. Due to the multiple X-line reconnections inside the PIC code, FTEs are generated quasi-periodically at low latitudes, then propagate toward the cusps. We briefly describe the evolution of the FTEs here:

* During the growth of an FTE, its cross-section increase, and its length extends along the dawn-dusk direction (from $t = 100$ s to $t = 150$ s in Fig. 14).

* Since its ambient plasma flow speed varies, an FTE may become tilted ($t = 240$ s in Fig. 14).

* There may be multiple FTEs on the magnetopause, and a few FTEs can merge into one (from $t = 320$ to $t = 660$ s in Fig. 14).

* FTEs can be dissipated at high latitudes due to the reconnection between the FTE magnetic field and the cusp field lines (Fig. 5 of Chen et al. 2017).

Chen et al. (2020) simulated the GEM dayside kinetic processes challenge event, and compared simulation results with both MMS observations and ground-based SuperDARN observations. The MHD-EPIC simulation shows there are usually multiple X-lines at the magnetopause, and the expanding speed of the X-line endpoints is comparable with SuperDARN observations.
<figure><img src="image_16.png" /><figcaption>**Fig. 14** Evolution of FTEs. Viewed from the Sun, a series of snapshots are shown with magnetic field lines colored by ion velocity $u_{iz}$[km/s]. Adapted from Chen et al. (2017)</figcaption></figure>

The MHD-EPIC model has also been applied to study magnetotail reconnection. Wang et al. (2022b) found the MHD-EPIC simulation can produce global-scale magnetospheric sawtooth-like oscillations periodically even under steady solar wind conditions, while the ideal- and Hall-MHD simulations do not produce such variations. It suggests that kinetic reconnection physics may play an important role in driving sawtooth oscillations. Recently, the development of the MHD-AEPIC model has enabled us to simulate storm events with a dynamic PIC region that covers the highly dynamic magnetotail reconnection sites.

Currently, research employing MHD-AEPIC focuses on modeling extreme geomagnetic storm events. Extreme events occur infrequently, which makes it difficult to validate MHD models employing simple numerical diffusion to approximate reconnection physics. Using
a higher-fidelity model, such as MHD-AEPIC, can improve the reliability of simulations of extreme events.

# 7 Kglobal: Particle Acceleration Self-Consistently Embedded in a Fluid Model

Solar flares convert magnetic energy into particle energy via magnetic reconnection. Observations of power-law tails in particle distribution functions imply that a large fraction of the released energy goes to energetic (i.e., non-thermal) electrons and ions (Warmuth and Mann 2016). However, the particle spectra found in particle-in-cell (PIC) simulations of reconnection in the relevant regime typically do not form power-laws, except in the limit of extremely low upstream plasma $\beta$ (Dahlin et al. 2015, 2017; Zhang et al. 2021). Why? With structures extending $\sim 10^4$ km and a Debye length of $\sim 1$ cm (for $n \sim 10^{10}$ cm$^{-3}$ and $T_e \sim 100$ eV), the corona spans ten orders of magnitude in physical scale. Explicit PIC models must resolve kinetic scales and hence can only simulate a tiny fraction of the macroscopic domain. The dependence of the Larmor radius on energy means nonthermal particles can quickly acquire orbits that approach the size of the simulation domain, halting further energy gains.

In contrast, MHD simulations study macroscopic domains with a fluid description that averages over small spatial and temporal scales. Following test particles in the MHD fields produces information about how particles gain energy but, without feedback coupling the particles and the fields, runaway energy gain can occur so that the system as a whole does not conserve total energy. It is possible to embed PIC models into MHD descriptions at selected locations, but such models presume that particle energy gain occurs in the vicinity of magnetic nulls, which is not consistent with the development and interaction of macroscale magnetic islands or the development of turbulence in large-scale current layers.

The *kglobal* model incorporates the physics necessary to explore particle energization from both the PIC and MHD descriptions (Drake et al. 2019; Arnold et al. 2019). The fundamental question is whether kinetic-scale boundary layers play an essential role in particle energy gain – or if they can be ordered out of the equations to facilitate simulations of macroscale systems. Kinetic boundary layers control the regions where $E_\|$, the component of the electric field parallel to the magnetic field, is non-zero. However, Fermi reflection rather than $E_\|$ is the dominant driver of energetic particles (Dahlin et al. 2016; Li et al. 2019a). Particle energy gain from Fermi reflection takes place over macro-scale regions and occurs even where $E_\| = 0$. As a consequence, kinetic-scale boundary layers are not required to describe the non-thermal energization in macroscale systems.

In order to keep the physics most important for describing particle energization while still being able to model macroscale systems, *kglobal* combines aspects of these descriptions. Guiding center particles move through a computational domain that includes a grid for fluid quantities. These particles feed back on the fluid through their gyrotropic pressure tensor. They can be small in number density but can contribute a pressure comparable to the pressure of the reconnecting magnetic field. The entire system conserves total energy.

The basic version of *kglobal* includes three species: fluid ions, fluid electrons, and particle electrons (the latter of which form the nonthermal population). We note that in a very strict sense the model may be viewed as multifluid plus guiding center methods. However,
the motion of the fluid electrons is severely constrained: their density is determined by quasi-neutrality, their parallel velocity from the requirement that parallel currents vanish, and their perpendicular velocity equal to the $E \times B$ flow (as it is for every species). The electromagnetic fields follow the usual MHD equations

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \nabla \times \mathbf{E}_{\perp} \qquad \mathbf{E}_{\perp} = -\frac{1}{c} \mathbf{v}_i \times \mathbf{B} \qquad (67) $$

and the ion fluid satisfies the usual MHD continuity equation

$$ \frac{\partial n_i}{\partial t} + \nabla \cdot n_i \mathbf{v}_i = 0 \qquad (68) $$

and energy equation

$$ \frac{d}{dt} \left( \frac{P_i}{n_i^\gamma} \right) = 0 \qquad (69) $$

The ion momentum equation takes the form

$$ \rho_i \frac{d\mathbf{v}_i}{dt} = \frac{1}{c} \mathbf{J} \times \mathbf{B} - \nabla P_i - \nabla_{\perp} P_{ef} - m_e n_{ef} v_{ef}^2 \boldsymbol{\kappa} + en_i E_{\parallel} \mathbf{b} - (\nabla \cdot \mathbf{T}_{ep})_{\perp} \qquad (70) $$

in which the left-hand side and first terms on the right-hand side are the same as in MHD ($P_i$ and $P_{ef}$ are the ion and fluid electron pressure, respectively). However the final terms on the right-hand side include the curvature $\kappa = \mathbf{b} \cdot \nabla \mathbf{b}$, large-scale parallel electric field $E_{\parallel}$, and particle electron stress tensor $T_{ep}$ and quantify the self-consistent back-reaction of the particles on the system.

The perpendicular motion of the particle electrons is given by the conservation of the first adiabatic invariant

$$ \mu_{ep} = \frac{p_{ep\perp}^2}{2B} = \text{const.} \qquad (71) $$

while the parallel motion satisfies

$$ \frac{d}{dt} p_{e\parallel} = p_{e\parallel} \mathbf{v}_E \cdot \boldsymbol{\kappa} - \frac{\mu_e}{\gamma_e} \mathbf{b} \cdot \nabla B - eE_{\parallel} \qquad (72) $$

This equation includes a contribution from the parallel electric field given by

$$ E_{\parallel} = -\frac{1}{n_i e} \left( \mathbf{B} \cdot \nabla \left( \frac{m_e n_{ef} v_{ef}^2}{B} \right) + \mathbf{b} \cdot \nabla P_c + \mathbf{b} \cdot \nabla \cdot T_{ep} \right) \qquad (73) $$

Finally, the fluid electron density enforces quasi-neutrality

$$ n_{ef} = n_i - n_{ep} \qquad (74) $$

the parallel flow eliminates parallel currents

$$ n_{ef} v_{ef\parallel} = n_i v_i \parallel - n_{ep} v_{ep\parallel} \qquad (75) $$
<figure><img src="image_17.png" /><figcaption>**Fig. 15** Energetic electron spectra from kglobal. A log-log plot of the electron differential density F(W) versus energy (W) at multiple times from a reconnection simulation with a guide field $B_g/B_0 = 0.25$. A power-law develops after $t/\tau_A \sim 3-5$. Inset: The late-time F(W) for several guide fields, illustrating the dependence on the ratio of the guide-to-ambient magnetic field. Adapted from Arnold et al. (2021)</figcaption></figure>

and the pressure equation takes the usual form

$$ \frac{d}{dt} \left( \frac{P_{ef}}{n_{ef}^{\gamma}} \right) = 0 \qquad (76) $$

A full derivation of these equations is given in Drake et al. (2019) and Arnold et al. (2019).

Simulations with these equations pass several tests. They describe the linear propagation of stable, circularly polarized Alfvén waves and the linear growth of firehose modes. The latter plays an important role in controlling the feedback of energetic particles during magnetic reconnection since magnetic tension is suppressed on the approach to firehose marginal stability. In addition, they accurately capture the dynamics of electron acoustic waves and describe the suppression of transport of hot electrons parallel to the ambient magnetic field. The inclusion of the large scale $E_\parallel$ is important in describing the development of return currents that form as hot electrons escape from regions of electron acceleration in macroscale energy release events such as flares (Egedal et al. 2012).

Reconnection simulations with kglobal have produced power-law spectra of energetic electrons that extend nearly three decades in energy, while simultaneously generating the super-hot thermal electrons characteristic of flare observations (Arnold et al. 2021). Figure 15 shows the electron energy spectrum for a typical simulation. Electrons in the initial Maxwellian distribution (black curve) transform into a nonthermal spectrum in a few Alfvén crossing times ($\tau_A$). Consistent with observations, the total energy content of the nonthermal electrons can exceed that of the hot thermal electrons even though the number density does not. The strength of the ambient out-of-plane guide field strongly impacts the energy content and power-law index of the nonthermal electrons (see inset of Fig. 15): the guide field increases the radius of curvature of a reconnected field line, thereby weakening Fermi reflection (Drake et al. 2006). In contrast, the size of the global system has relatively little influence.

The governing equations of kglobal can be extended to include the contributions of non-thermal (particle) ions. Unlike for electrons, whose small mass can be used to simplify the equations, the ion inertia can not be neglected and must be included. Recent work has incorporated these equations into the computational model and early results reveal the simultaneous development and evolution of extended electron and proton power law distributions (Yin et al. 2024a,b).
<figure><img src="image_18.png" /><figcaption>**Fig. 16** Sketch of the typical sampling of the plasma distribution function $f(x, v)$ adopted in PIC (left) and Eulerian (methods). Adapted from Finelli (2022)</figcaption></figure>

# 8 Vlasov

## 8.1 Overview

Eulerian Vlasov-Maxwell numerical simulations are a useful tool for investigating fundamental kinetic-scale plasma processes, such as turbulence and magnetic reconnection, as well as the interaction between the solar wind and planetary magnetospheres.

Thanks to the clean description of the plasma dynamics in the entire phase space at the expense of a larger computational cost, Eulerian algorithms complement well Particle-In-Cell (PIC) codes. The almost noise-free description of velocity space is generally guaranteed by the discretization of the plasma distribution function on a six-dimensional phase-space grid characterized by collocation points in both physical and velocity space. On the other hand, PIC methods suffer from the intrinsic stochastic shot noise which becomes especially relevant at small scales and in cases where the number of particles per cell is not large. However, in Eulerian methods, setting a six-dimensional grid in the entire phase space dramatically increases the computational cost. The bottleneck is generally constituted by the memory necessary to store the plasma distribution function, as shown in the following simple example.

The main difference when sampling the plasma distribution function in PIC and Eulerian approaches is depicted in Fig. 16. In PIC methods, the grid is defined only on the physical space **x**, and the distribution function is sampled through macroparticles (blue circles), each one representative of a large number of effective plasma particles. In the Eulerian approach, the grid is defined on the entire phase space (**x**, **v**) and the distribution function is known on this ensemble of grid points (red circles). In both PIC and Eulerian methods, the electromagnetic fields and the moments of the distribution function (e.g., density, bulk speed, etc.) are defined on the physical-space grid **x**.

Let us imagine a generic phenomenon occurring in a plasma composed of protons and electrons that requires a physical-space grid discretized with 512³ points. Since the electromagnetic fields have the same memory requirements in both PIC and Eulerian methods, we neglect them here for simplicity. PIC simulations indicate that the plasma dynamics is overall well described with ~1000 particles per cell. The memory (assuming double-precision
variables) required to store particles positions and velocities is ~ 6 TB. Similarly, Eulerian Vlasov simulations with a homogeneous Cartesian grid in velocity space, generally adopt at least ~ 51³ velocity-space points to well describe the fine details of the velocity plasma distribution function and preserve mass and entropy conservation, thus requiring about ~ 260 TB of memory. Hence, PIC simulations can usually be performed in larger physical-space computational boxes compared to Eulerian ones. Moreover, 3D simulations are more easily achievable with PIC methods while they remain often prohibitive for Eulerian codes, or require more advanced methods such as adaptive mesh refinement or sparse velocity space techniques to become tractable.

In the following, we will introduce two Vlasov-Maxwell algorithms that have been intensively used by the scientific community — namely the Vlasiator code (von Alfthan et al. 2014; Palmroth et al. 2018a) and the Hybrid Vlasov-Maxwell HVM code (Valentini et al. 2007) — and, then, discuss the results of global and local simulations of magnetic reconnection based on these two codes.

## 8.2 Models and Algorithms

Picturing a virtual journey from large to small scales, the first relevant Vlasov-Maxwell model widely adopted for performing Eulerian simulations is the hybrid Vlasov-Maxwell model. The hybrid model considers protons as a kinetic species, while electrons are a background fluid. It is a low-frequency approximation of the full Vlasov-Maxwell system of equations that assumes quasi-neutrality and neglects the displacement current (Mangeney et al. 2002). Faraday's law is used to evolve the magnetic field, while the electric field is provided by the generalized Ohm's law which includes the Hall term, the electron pressure gradient term, and possibly the terms related to electron inertia (see, e.g., Valentini et al. 2007 for further details about the generalized Ohm's law).

The following paragraphs present Vlasiator (von Alfthan et al. 2014; Palmroth et al. 2018a) and the Hybrid Vlasov-Maxwell (HVM) (Valentini et al. 2007) codes as particular examples of algorithms based on the Vlasov approach. The reader is invited to refer to recent review papers (Califano and Cerri 2022; Palmroth et al. 2018a) for more details on how to solve the Vlasov equation numerically.

The hybrid Vlasov-Maxwell model has been adopted for global simulations of the interaction between the solar wind and the Earth's magnetosphere through the Vlasiator algorithm (von Alfthan et al. 2014; Palmroth et al. 2018a). The Vlasiator code adopts a splitting algorithm to decompose the six-dimensional Vlasov equation into a set of two three-dimensional advection equations (Strang 1968) in physical and velocity space, respectively. The solution of each advection equation is obtained by a semi-Lagrangian method (Zerroukat and Allen 2012), which relieves from the strict limitations to the time step length posed by the CFL condition<fnref n="1" /> in velocity space acceleration in regions of strong magnetic field. Position space is discretized on a cell-adaptive Cartesian grid (Honkonen et al. 2013; Ganse et al. 2023) and at each position in space the velocity-space grid is stored on a uniform Cartesian grid. Vlasiator developed a sparse velocity-space method in which only regions of the velocity distribution function above a set phase-space density are stored and propagated (von Alfthan et al. 2014), yielding a gain of two orders of magnitude in terms of memory and computations. This technique made two-dimensional (2D position space periodic in the third dimension, 3D velocity space) magnetospheric simulations possible,

<footnote n="1">¹The CFL condition dictates that an element of the distribution function $f$ may not cross more than one velocity space grid spacing in one time step under the acceleration due to the Lorentz force.</footnote>
as well as quasi-three dimensional simulations with a very limited extent in the third dimension (Pfau-Kempf et al. 2020). The implementation of adaptive mesh refinement in position space, allowing to focus resolution on regions of interest while saving computations in less-resolved regions, is what made full three-dimensional, global magnetospheric simulations achievable with Vlasiator on modern, bleeding-edge supercomputers (Grandin et al. 2023; Ganse et al. 2023; Palmroth et al. 2023). The electric and magnetic fields are propagated using an upwind constrained transport method (Londrillo and Del Zanna 2004) with divergence-free magnetic field reconstruction (Balsara 2009) on a uniform Cartesian grid matching the finest refinement level of the Vlasov spatial grid, requiring a dedicated coupling scheme between the grids (Papadakis et al. 2022). Note that, despite the sparse velocity-space grid, Vlasiator also requires a huge amount of RAM memory similar, as order of magnitude, with the simple estimate discussed above.

The hybrid Vlasov-Maxwell model has also been adopted for local simulations of plasma turbulence at sub-proton scales using the HVM code (Valentini et al. 2007) (see Califano and Cerri 2022 for a recent review). The HVM code reduces the six-dimensional Vlasov equation to a set of six one-dimensional advection equations. Each equation is then solved through the van Leer method (van Leer 1977). Fields are computed through the Current-Advance Method (CAM). The grid is homogeneous in both physical and velocity space. Periodic boundary conditions are implemented in physical space. In velocity space, the proton distribution function is set to zero after a large number of thermal speeds, thus ensuring that mass conservation is preserved. Numerical resolution in physical space is usually about 4-5 gridpoints per ion skin depth, while velocity-space grid resolution is of about 0.2–0.25 ion thermal speed. The latter condition is an optimal compromise to (i) guarantee that velocity-space filamentation and distortion of the ion distribution function are well observed in velocity space, thus allowing to investigate energy conversion in the entire phase space, (ii) preserve the entropy conservation to an excellent value (≈ 0.1%), and (iii) maintain simulations numerically feasible.

Moving towards smaller scales, fully-kinetic Vlasov-Maxwell Eulerian algorithms have been recently implemented to describe electron-scale dynamics. Given the larger computational cost of Eulerian simulations with respect to PIC methods, the former are generally more recent than the latter and possibly implement different assumptions to simplify the Maxwell equations. The full Maxwell system has been retained in several codes (Umeda et al. 2009, 2010; Delzanno 2015; Ghizzo et al. 2017; Juno et al. 2018; Pezzi et al. 2019a; Allmann-Rahn et al. 2022). However, different approximations of the Maxwell equations have been proposed to alleviate the CFL constraint which sets a very small time step when the wave phase speed approaches the speed of light. In this regard, Wiegelmann and Büchner (2001) neglect the displacement current while allowing for charge separation, while Tronci and Camporeale (2015) ignore both the displacement current and the charge separation. Yet a different approach neglects only the transverse part of the displacement current responsible for ordinary mode propagating at the speed of lights (Schmitz and Grauer 2006a; Pezzi et al. 2019a; Shiroto 2023). Finally, an original approach has been developed based on the Vlasiator model, whereby a small section of interest from an ion-hybrid Vlasiator run is used to initialize an electron-hybrid setup in which the ions are kept static while the electron distribution function evolves (Battarbee et al. 2021). This so-called eVlasiator approach has successfully reproduced properties of electron distributions observed in the vicinity of reconnection diffusion regions (Alho et al. 2022).
## 8.3 Examples of Applications Focused on the Study of Magnetic Reconnection

### 8.3.1 Local Simulations

In this section, we will present key results relevant to magnetic reconnection that have been obtained with the Hybrid Vlasov Maxwell code (HVM). Further treatments, based on fully-kinetic Eulerian Vlasov-Maxwell simulations, will not be covered in detail in this article, but the reader is referred to the works of Schmitz and Grauer (2006b), Inglebert et al. (2011), Zenitani and Umeda (2014), Sarrat et al. (2017), Pezzi et al. (2019a), as well as Table 2 in the review by Palmroth et al. (2018a) which lists works using Vlasov-based methods in space and astrophysics.

The HVM code, which retains alpha particles (Perrone et al. 2012; Valentini et al. 2016) and inter-particle collisions (Pezzi et al. 2019b), has been used for years to investigate plasma processes occurring at ion kinetic scales. It has been massively employed to study the properties of plasma turbulence (Valentini et al. 2010; Servidio et al. 2012, 2014, 2015; Cerri et al. 2017), showing that turbulent fluctuations generate manifestly non-Maxwellian proton distribution functions (Greco et al. 2012). This emergent velocity-space complexity has been envisioned as a cascade process occurring in velocity space (e.g., Tatsuno et al. 2009; Schekochihin et al. 2016; Servidio et al. 2017): HVM results have allowed to characterize it in a full Vlasov system rather than in the gyrokinetic approximation (Cerri et al. 2018; Pezzi et al. 2018, 2021b). Characterizing non-equilibrium plasma distribution functions is significant to understanding energy transfer and dissipation processes occurring at ion kinetic scales in nearly-collisionless plasmas such as the solar wind (Matthaeus et al. 2020; Cassak et al. 2023), as also reported in different studies based on the HVM code (Sorriso-Valvo et al. 2018; Pezzi et al. 2019c, 2021a; Fadanelli et al. 2021). In the perspective of the Holloway and Dorning (1991) work showing that non-Maxwellian plasmas can support the propagation of undamped plasma waves, the HVM code has been adopted to study the onset of a novel type of electrostatic fluctuations triggered by trapped ions (Valentini et al. 2011a,b, 2014).

One of the first studies employing the HVM code for investigating magnetic reconnection reported the onset of a fast reconnection process obtained as a result of magnetic islands developed by the electromagnetic current filamentation (Califano et al. 2001). In the following years, despite the large number of studies adopting the HVM code, the vast majority of the research work mostly focused on the investigation of fully developed plasma turbulence. More recently, Finelli et al. (2021) investigated the magnetic reconnection in a similar manner as PIC-based studies discussed in Sect. 5, that is, modeling an isolated Harris-like current sheet (Harris 1962b) in equilibrium or pressure balance. Such a current sheet, usually doubled in the physical-space domain to accommodate for periodic boundary conditions, quickly starts reconnecting thanks to an initial perturbation of proton density and/or current (thus magnetic field).

In particular, Finelli et al. (2021) compared results from three different models (i) the HVM model with (isotropic) isothermal electrons including finite electron-inertia, (ii) a modified HVM model, called hybrid-Vlasov-Landau-fluid (HVLF); (iii) a fully-kinetic PIC code (iPIC3D Markidis et al. 2010). The HVLF model is equipped to include anisotropies of the gyrotropic electron pressure with a Landau-fluid (LF) closure for the transport of the gyrotropic electron thermal energy along magnetic field lines (Sulem and Passot 2015). Using these three models, Finelli et al. (2021) performed 2D-3V magnetic reconnection simulations with moderate guide field ($B_g = 0.25 B_0$, where $B_0$ is the asymptotic magnetic field) and with reduced mass ratio $m_p/m_e = 100$. The initial setup consists of a double Harris