A set of equations is developed that extends the macroscale magnetic reconnection simulation model kglobal to include particle ions. The extension from earlier versions of kglobal, which included only particle electrons, requires the inclusion of the inertia of particle ions in the fluid momentum equation. The new equations will facilitate the exploration of the simultaneous nonthermal energization of ions and electrons during magnetic reconnection in macroscale systems. Numerical tests of the propagation of Alfvén waves and the growth of firehose modes in a plasma with anisotropic electron and ion pressure are presented to benchmark the new model.
I. INTRODUCTION
Magnetic reconnection is a fundamental process that occurs in a diverse array of astrophysical and laboratory plasmas.^{1–7} It allows for rapid reconfiguration of the magnetic field, facilitating the release of magnetic energy into kinetic and thermal energy and nonthermal particle energization.^{8,9} An important consequence of magnetic reconnection is the efficient acceleration of charged particles into powerlaw (i.e., nonthermal) distributions that extend many decades in energy.^{3,10–13} Such acceleration events are of significant interest not only because of their scientific relevance in astrophysical phenomena such as solar flares and magnetospheric substorms but also because of the implications for controlled fusion experiments^{14} and space weather.^{15}
Many computational simulations, most notably those performed with particleincell (PIC) and magnetohydrodynamic (MHD)^{16–18} models, have limitations that compromise their utility for studying nonthermal particle energization in such systems. PIC simulations correctly model kinetic physics by resolving kinetic scales,^{19–24} but the spatial scales associated with solar flares can exceed 10^{4} km, while the Debye length—a measure of the scale of kinetic effects—can be of the order of centimeters, a scale separation of 10^{10}. Hence, computational costs drastically limit the ability of PIC simulations to simulate the extensive range of scales occurring in magnetic reconnection events. A consequence of this inadequate separation of scales is the demagnetization of energetic electrons and ions in PIC simulations of reconnection, which in turn inhibits particle energy gain and prevents the formation of the extended power laws that are observed in nature.^{23,24}
Conversely, MHD models employ a fluid approximation and can simulate largescale phenomena. While capable of representing macroscopic events, MHD simulations frequently assume the plasma is thermal (i.e., represented by a single temperature) and hence by definition cannot produce nonthermal particles. The electric and magnetic fields from MHD models have been used to explore nonthermal energization of test particles. However, these test particles do not interact selfconsistently with the fields, so energy is not conserved, which limits the test particle model's fidelity.^{25}
The original kglobal model combined aspects of the PIC and MHD descriptions to address electron acceleration in largescale astrophysical reconnection.^{26,27} Particle electrons are distributed on the MHD grid, as in PIC models, but are evolved with MHD fields in the guidingcenter approximation in which all kinetic scales including the Larmor radii of particles are ordered out of the equations. Notably, kglobal includes the effects of Fermi reflection within evolving magnetic flux ropes and the pressure anisotropy that feeds back on the MHD fields to reduce the tension force that drives magnetic reconnection. Magnetic reconnection simulations with kglobal produced, for the first time, extended powerlaw spectra of nonthermal electrons that aligned with existing observational data.^{28}
However, magnetic reconnection can simultaneously produce nonthermal electrons and nonthermal protons.^{3,29} Hence, developing a selfconsistent simulation model is vital to elucidate the mechanisms of ion and electron energization during magnetic reconnection and the resulting partitioning of energy between the two species. In this paper, we extend the kglobal equations to include particle ions, which requires the inclusion of the inertia of the particle ions in the bulk fluid ion momentum equation. The resulting equations will facilitate the exploration of the selfconsistent energization of both ions and electrons during magnetic reconnection. In Sec. II, we introduce the equations that describe the model. In Sec. III, we present test results for Alfvén wave propagation in a system with a finite pressure anisotropy and the firehose instability. In Sec. IV, we present an overview of the resulting equations and implications for exploring particle energization during reconnection in macroscale systems.
II. MODEL EQUATIONS
For completeness, we first outline the governing equations for the original kglobal model before turning to the additions made in the upgraded version.
A. kglobal with particle electrons
B. kglobal with particle electrons and ions
In the upgraded version of kglobal, particle ions are introduced as a new species and, as is the case for the particle electrons, are treated in the guidingcenter limit. Similar to the original model, the populations of the four plasma species do not change after initialization. The basic ordering of the electrons and the equations describing their dynamics are basically unchanged from the original model. Minor modifications to the charge neutrality and zero current conditions are discussed below, but the details of the electron dynamics that were discussed in earlier papers are not repeated here.
The ratios $nip/nif$ and $nep/nef$ must be chosen at t = 0, but they are free to evolve in space and time during simulations. The electron spectra reported in previous kglobal reconnection simulations were insensitive to $nep/nef$,^{28} but the sensitivity of reconnection simulations with both particle electrons and ions to these ratios will also need to be explored.
To summarize, the kglobal model with particle electrons and ions consists of Eqs. (1), (2), (6)–(9), (14)–(17), and (20)–(30). Because of the complexity of Eqs. (26) and (27), an alternative is to advance the equation for the centerofmass velocity in Eq. (B11) and calculate the ion fluid velocity by subtracting the ion particle flux.
III. ALFVÉN WAVE PROPAGATION AND GROWTH OF FIREHOSE MODES IN AN ANISOTROPIC PLASMA
Particle energy gain in reconnection is dominantly parallel to the ambient magnetic field, with the implication that the parallel pressure exceeds the perpendicular pressure in the heated plasma as reconnection develops. As a consequence, the cores of magnetic islands approach the marginal firehose condition at late time in simulations. The tension in a bent magnetic field line, which is the fundamental driver of reconnection, goes to zero at the marginal firehose condition. Thus, the pressure anisotropy that develops as reconnection proceeds regulates reconnection and particle energy gain.^{26,27,31,32} Therefore, to ensure that the updated kglobal model properly reproduces the influence of pressure anisotropy on the magnetic field dynamics, we benchmark the model with a circularly polarized Alfvén wave mode propagating along a uniform magnetic field in a system with pressure anisotropy ( $P\u2225\u2260P\u22a5$), a problem with a wellestablished solution.^{8} With increased pressure anisotropy, the reduction in magnetic tension of a propagating Alfvén wave reduces the wave phase speed. In a second benchmarking test on the dynamics of the pressure anisotropy of both particle ions and electrons, we explore the growth of the firehose instability.
The upgraded kglobal model builds upon the original kglobal framework, the details of which can be found in Refs. 26 and 27. The code is written in normalized units. A reference magnetic field strength $B0$ and the initial total ion density $ni0=nip+nif$ define the Alfvén speed, $CA=B02/(4\pi mini0)$. Lengths and times are normalized to an arbitrary macroscale length, $L$, and the Alfvén crossing time, $\tau A=L/CA$, respectively. Electric fields and temperatures are normalized to $CAB0/c$ and $miCA2$, respectively. The kglobal model employs secondorder diffusivities and fourthorder hyperviscosities ( $\u221d\nu \u22074$) in each fluid equation to mitigate noise at the grid scale.^{26}
The simulations were performed within a twodimensional (2D) spatial domain, coupled with a threedimensional velocity space, featuring an initially uniform magnetic field of magnitude B_{0} oriented along the x axis. The domain measures $2\pi L\xd7\pi L$ with $512\xd7256$ cells with periodic boundary conditions on all sides. Each grid cell is initialized with 200 selfconsistent (i.e., not test) particles, 100 ions and 100 electrons. The initial density ratios everywhere are $nip/nif=nep/nef=1/3$. The iontoelectron mass ratio is set to 25, although the results of the simulations are insensitive to this value. The scalar temperature of both fluid species is set to 0.1. For the particle ions and electrons, the perpendicular temperature is initially set to 0.1 while the parallel temperature is adjusted to control the magnitude of the anisotropy.
The results are shown in Fig. 1, where $Cp$, normalized to $CA$, is plotted against $\alpha $. There is excellent agreement between the code results and linear wave theory, highlighting the precision of the kglobal model and suggesting that the impact of ion and electron pressure anisotropy on reconnection dynamics will be properly included in the new model.
In the second benchmark, the rate of growth of the firehose instability was explored with a preset pressure anisotropy above the instability threshold ( $\alpha 2=\u22120.2$). The system was initialized with small amplitude periodic disturbances with seventeen distinct wavelengths, defined by $k=m/2\pi L$, where m is the mode number. The simulations incorporated a kinematic viscosity of $\nu =5.0\xd710\u22125$. The growth rate, which includes the pressure anisotropy drive and viscous damping, is given by $\gamma =kCA\alpha \u2212\nu k4$. The fourthorder hyperviscosity sets the upper bound of instability at small spatial scales. Figure 2 presents the expected theoretical growth rate (represented by a solid blue line) and the growth rates found in the simulations for a range of unstable modes. The red crosses indicate systems in which the electrons and ions contribute equally to the total anisotropy. To explore the effect of unequal pressure anisotropies, we also performed several simulations where the ions contribute 3/4 and the electrons contribute 1/4 of the total anisotropy (green squares). The growth rates measured in the unequal anisotropy simulations closely follow those from equal anisotropy simulations and the results from all the simulations are consistent with the predictions of linear theory. Note that there are no spurious unstable modes observed at short scales such as those that arise in some numerical models when the mode frequency exceeds the cyclotron frequency.^{34} All waves in the kglobal model are below the cyclotron frequency.
IV. CONCLUSIONS AND DISCUSSION
The kglobal model has been developed to explore magneticreconnectiondriven particle acceleration in macroscale systems. The basis for the model is that particle energization is controlled by the dynamics of reconnectiondriven flux ropes and turbulence at the macroscale rather than in kineticscale boundary layers. In the present paper, we develop equations that extend the kglobal to include the dynamics of particle ions. This requires the inclusion of the inertia of the ion particles into the fluid momentum equation. The model correctly captures the development of the pressure anisotropy on the dynamics of the magnetic field, which has a significant influence on reconnection dynamics.
Our future work will center on using the upgraded kglobal model to investigate the relative acceleration and energization of electrons and ions during magnetic reconnection in macroscale systems. Fundamental questions will be addressed by this model. Can reconnection produce extended powerlaw distributions of energetic ions and electrons? What is the relative partition of energy into nonthermal ions and electrons? Can reconnection produce ion and electron energy spectra consistent with in situ measurements from NASA's Magnetospheric Multiscale Mission in the Earth's magnetotail and from the Parker Solar Probe in the near solar environment?
ACKNOWLEDGMENTS
We acknowledge support from the FIELDS team of the Parker Solar Probe (NASA Contract No. NNN06AA01C), the NASA Drive Science Center on Solar Flare Energy Release (SolFER) under Grant No. 80NSSC20K0627, NASA Grant Nos. 80NSSC20K1277 and 80NSSC22K0352, and NSF Grant No. PHY2109083. The simulations were carried out at the National Energy Research Scientific Computing Center (NERSC). The data used to perform the analysis and construct the figures for this paper are preserved at the NERSC High Performance Storage System and are available upon request.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Zhiyu Yin: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). J. F. Drake: Funding acquisition (lead); Methodology (equal); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). M. Swisdak: Methodology (equal); Project administration (equal); Software (equal); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: DERIVATION OF THE FLUID ION MOMENTUM EQUATION

$mi$ and $me$ represent the masses of ions and electrons, respectively;

$vip$, $vef$, $vep$ and $vif$ denote the velocities of fluid ions, fluid electrons, particle electrons, and particle ions, respectively;

$Jif$, $Jef$, $Jep$, and $Jif$ represent the current densities corresponding to fluid ions, fluid electrons, particle electrons, and particle ions, respectively;

$Pif$ and $Pef$ denote the (scalar) pressures for fluid ions and fluid electrons, respectively; and

$Tep$ and $Tip$ represent the stress tensor for particle electrons and particle ions, respectively.