A Computational Model for Ion and Electron Energization during Macroscale Magnetic Reconnection

A set of equations are developed that extend 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 non-thermal energization of ions and electrons during magnetic reconnection in macroscale systems. Numerical tests of the propagation of Alfv´en waves and the growth of firehose modes in a plasma with anisotropic electron and ion pressure are presented to benchmark the new model.


INTRODUCTION
Magnetic reconnection is a fundamental process that occurs in a diverse array of astrophysical and laboratory plasmas (Tsuneta 1996;Yamada et al. 1997;Lin et al. 2003;Gosling et al. 2005;Phan et al. 2006;Burch et al. 2016;Gary et al. 2018).It allows for the rapid reconfiguration of the magnetic field, facilitating the release of magnetic energy into kinetic, thermal energy and nonthermal particle energization (Parker 1957;Petschek 1964).An important consequence of magnetic reconnection is the efficient acceleration of charged particles into power-law (i.e., non-thermal) distributions that extend many decades in energy (Lin et al. 2003;Krucker et al. 2010;Oka et al. 2013;Øieroset et al. 2002;Ergun et al. 2020).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 (Yamada et al. 1994) and space weather (Angelopoulos et al. 2020).
Many computational simulations, most notably those performed with particle-in-cell (PIC) and magnetohydrodynamic (MHD) (Biskamp 1986;Karpen et al. 2012;Dahlin et al. 2022) models, have limitations that compromise their utility for studying non-thermal particle energization in such systems.PIC simulations correctly model kinetic physics by resolving kinetic scales (Shay et al. 1998;Hesse et al. 1999;Shay et al. 2007;Daughton et al. 2011;Li et al. 2019;Zhang et al. 2021) 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 (Li et al. 2019;Zhang et al. 2021).
Conversely, MHD models employ a fluid approximation and can simulate large-scale 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 non-thermal particles.The electric and magnetic fields from MHD models have been used to explore non-thermal 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 (Onofri et al. 2006).
The original kglobal model combined aspects of the PIC and MHD descriptions to address electron accel-eration in large-scale astrophysical reconnection (Drake et al. 2019;Arnold et al. 2019).Particle electrons are distributed on the MHD grid, as in PIC models, but are evolved with MHD fields in the guiding-center 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 power-law spectra of non-thermal electrons that aligned with existing observational data (Arnold et al. 2021).
However, magnetic reconnection can simultaneously produce non-thermal electrons and non-thermal protons (Lin et al. 2003;Emslie et al. 2012).Hence, developing a self-consistent 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 self-consistent energization of both ions and electrons during magnetic reconnection.In Section 2 we introduce the equations that describe the model.In Section 3 we present test results for Alfvén wave propagation in a system with a finite pressure anisotropy and the firehose instability.In Section 4 we present an overview of the resulting equations and implications for exploring particle energization during reconnection in macroscale systems.

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.

kglobal with particle electrons
The original model includes three distinct plasma species: the ion fluid (number density n i ), the electron fluid (n ef ), and particle electrons (n ep ) (Drake et al. 2019;Arnold et al. 2019).No conversion between species is allowed.The particle electrons are governed by the guiding-center equations given in Northrop (1963).Their perpendicular velocity is determined by conservation of the magnetic moment: where p ep is the momentum and B is the magnitude of the local magnetic field.Their perpendicular motion is given by the usual MHD E × B drift while they stream parallel to the ambient magnetic field at their parallel velocity.The variation in the parallel velocity reflects contributions from Fermi reflection, magnetic mirroring, and the large-scale parallel electric field with v E the E×B drift velocity, b = B/B a unit vector in the direction of the magnetic field, κ = b •∇ b the field curvature, and γ e the relativistic Lorentz factor.The guiding-center electrons move through the MHD grid and, crucially, do not form a thermal distribution.
Their gyrotropic pressure tensor feeds back on the MHD momentum equation.
The ion fluid forms an MHD-like backbone that satisfies the usual continuity equation and a modified momentum equation where the ion stress tensor with inertial contributions is with P i the ion scalar pressure and I the unit tensor.
The corresponding electron fluid tensor is with P ef the electron fluid scalar pressure.Finally, T ep is the particle electron gyrotropic stress tensor, including inertial contributions, where T ep∥ is the component of the stress tensor along the magnetic field B, and P ep⊥ is the usual perpendicular pressure, where in the frame drifting with v E = cE × B/B 2 there are no perpendicular flows so f = f (x, p e∥ , p e⊥ , t).T ep∥ includes the mean parallel drifts of the particle electrons, with p e∥ being the electron particle parallel momentum.Note that in the ion momentum equation the electron fluid and particle inertia are neglected compared to the ion inertia on the left side of Eq. (4), while the electron fluid and particle pressures are retained on the right side of this equation.This result is based on an ordering in which the ion velocity scales as the Alfvén speed and the electron particle and fluid thermal velocities and parallel streaming velocities are of the order of the electron Alfvén speed, C Ae = m i /m e C A , where the Alfvén speed is defined as C A = B 2 /(4πm i n i ).In this ordering the electron momentum scales as m e C Ae = m e /m i m i C A and so is much smaller than that of the ion fluid.In contrast, the electron pressure is m e C 2 Ae = m i C 2 A , which is the same as the ion pressure.Thus, the pressure terms on the right side of Eq. ( 4) must be retained even though the electron inertia in this equation is neglected.Finally, the ions follow an adiabatic equation of state where the adiabatic index γ is taken to be 5/3.The number density of the electron fluid is given by while the parallel motion follows from the requirement of vanishing parallel current The zero current condition follows from the constraint on the current from Ampère's law where V d is the drift speed associated with the current, d i is the ion inertial scale length and L is the characteristic magnetic scale length.In the kglobal model all kinetic scales are ordered out, which yields the zero current condition in Eq. ( 12).The fluid electrons also follow an adiabatic equation of state The magnetic field is advanced via Faraday's law, with Finally, E ∥ is produced by magnetic-field-aligned gradients in the electron pressure and is given by (Arnold et al. 2019) Thus, Equations ( 1)-( 17) constitute the complete set of equations for the kglobal model with particle electrons.

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 guiding-center 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.
Due to the small mass of the electrons, their inertial terms in the fluid ion momentum equation were neglected in the original kglobal equations.That simplification is not possible for the particle ions in the upgraded model, which complicates the derivation and the final form.However, the perpendicular motion of the fluid and particle ions is still governed by the same E × B drift.As a consequence, the generic form of the ion momentum equation, can be written as with v if and n if the fluid ion velocity and density, v ip and n ip the particle ion velocity and density, and S the usual momentum driver, including convective inertial terms.The differences in the inertia between the fluid and particle ions arises solely from the parallel motion so the total momentum equation takes the form where n i = n if +n ip .The inertial terms parallel to B are easily evaluated, which leaves an equation for the fluid ion momentum that includes the inertia of the particle ions.The details of the evaluation of the parallel inertial terms are found in the Appendix.Most of the other equations in Section 2.1 remain unchanged, aside from minor notation adjustments (e.g., n i → n if + n ip ).Specifically, Equations ( 1) and ( 2) governing the motion of the particle electrons are unchanged.Equations ( 3) and ( 10) become equations for the ion fluid density, and ion fluid pressure, The guiding center equations for the particle ions are essentially identical to those of the electrons in Equations ( 1) and ( 2), and The ions are taken to be non-relativistic.Including the new species also requires changes in the expressions for the electron fluid density and parallel motion The adiabatic equation of state for the electron fluid, Equation ( 14), remains the same as do equations (15), and ( 16) for B and E ⊥ .Equation ( 17) for E ∥ remains the same because only the electrons contribute.
In both the original and upgraded model, electron inertia can be discarded compared with that of the ions in the ion momentum equation, as discussed in Sec.2.1.The distinction between the inertia of collective and individual particle electrons, along with their movements in relation to the ambient magnetic field and their specific interactions with the MHD fluid, is comprehensively discussed in Drake et al. (2019).However, as shown in equation ( 18), the inertia of the particle ions cannot be neglected.The derivation of the new fluid ion momentum equation is detailed in Appendix A, with the final result being where I i∥ is the difference in the parallel inertia between the fluid and particle ions, given by where n i = n ip + n if , and the ion fluid stress tensor, T if , is defined in Equation ( 5).The particle ion pressure tensor is given by with the parallel and perpendicular pressures calculated in the frame of the bulk motion of the particle ions, and with p ip the ion momentum.Note that in the inertia term on the left side of Equation ( 26), n i is the total ion density and not just the fluid ion density n if .In the limit n ip → 0, I i∥ → 0 and T ip → 0 and equation ( 26) reduces to equation (4).The ratios n ip /n if and n ep /n ef must be chosen at t = 0, but are free to evolve in space and time during simulations.The electron spectra reported in previous kglobal reconnection simulations were insensitive to n ep /n ef (Arnold et al. 2021), but the sensitivity of reconnection simulations with both particle electrons and ions to these ratios will also need to be explored.

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 (Drake et al. 2006;Drake & Swisdak 2012;Drake et al. 2019;Arnold et al. 2019).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 ∥ ̸ = P ⊥ ), a problem with a well-established solution (Parker 1957).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 the references (Drake et al. 2019;Arnold et al. 2019).The code is written in normalized units.A reference magnetic field strength B 0 and the initial total ion density n i0 = n ip + n if define the Alfvén speed, C A = B 2 0 /(4πm i n i0 ).Lengths and times are normalized to an arbitrary macroscale length, L, and the Alfvén crossing time, τ A = L/C A , respectively.Electric fields and temperatures are normalized to C A B 0 /c and m i C 2 A , respectively.The kglobal model employs second-order diffusivities and fourth-order hyperviscosities (∝ ν∇ 4 ) in each fluid equation to mitigate noise at the grid scale (Drake et al. 2019).
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πL × πL with 512 × 256 cells with periodic boundary conditions on all sides.Each grid cell is initialized with 200 self-consistent (i.e., not test) particles, 100 ions and 100 electrons.The initial density ratios everywhere are n ip /n if = n ep /n ef = 1/3.The ion-to-electron 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.
We introduce perturbations to the magnetic field and velocity in the form of a circularly polarized Alfvén wave, with a wavelength equal to the length of the sim- ulation box (k x = 2π/L).During a propagation interval of 3τ A , we measured the wave's velocity in order to compare it with the theoretical phase speed C p of an Alfvén wave, which is given by where α = 1 − 4π(P ∥ − P ⊥ )/B 2 , aligns with the Chew-Goldberger-Low (CGL) framework in the linear regime where pressure perturbations are negligible (Chew et al. 1956).Here, P ∥ and P ⊥ denote the combined ion and electron pressures along and perpendicular to the magnetic field, respectively.
The results are shown in Fig 1, where C p , normalized to C A is plotted against α.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 pre-set pressure anisotropy above the instability threshold (α 2 = −0.2).
The system was initialized with small amplitude periodic disturbances with seventeen distinct wavelengths, defined by k = m/2πL, where m is the mode number.The simulations incorporated a kinematic viscosity of ν = 5.0 × 10 −5 .The growth rate, which includes the pressure anisotropy drive and viscous damping is given by γ = kC A |α| − νk 4 .The fourth-order 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 (Tronci et al. 2014).All waves in the kglobal model are below the cyclotron frequency.

CONCLUSIONS AND DISCUSSION
The kglobal model has been developed to explore magnetic-reconnection-driven particle acceleration in macroscale systems.The basis for the model is that particle energization is controlled by the dynamics of reconnection-driven flux ropes and turbulence at the macroscale rather than in kinetic-scale 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 power-law 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?
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 80NSSC20K0627, NASA Grants 80NSSC20K1277 and 80NSSC22K0352 and NSF Grant 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.The derivation of the fluid ion momentum equation is complicated by the need to keep the inertia of the particle ions.We consider four distinct species, fluid ions, fluid electrons, particle electrons, and particle ions with number densities n if , n ef , n ep , and n ip , respectively, that satisfy the charge neutrality constraint where n if is advanced in time with a continuity equation.The motion of the particle species includes the E × B perpendicular to B and parallel dynamics at the relevant local velocity along B. The number densities n ip and n ep are computed on the grid so both particle species also effectively satisfy continuity equations.The electron fluid also satisfies a continuity equation because of charge neutrality and because the three other species satisfy continuity equations.Specifically Eq. (A1) for charge neutrality is used to determine n ef .In addition, each species satisfies a momentum equation: for particle electrons, and for particle ions.In these equations: • m i and m e : Represent the masses of ions and electrons, respectively.
• v ip , v ef , v ep and v if : Denote the velocities of fluid ions, fluid electrons, particle electrons, and particle ions, respectively.
• J if , J ef , J ep , and J if : Represent the current densities corresponding to fluid ions, fluid electrons, particle electrons, and particle ions, respectively.
• P if and P ef : Denote the (scalar) pressures for fluid ions and fluid electrons.
• T ep and T ip : Represent the stress tensor for particle electrons and particle ions.
The choice to use the stress tensor in the particle equations, rather than the mathematically equivalent pressure and convective terms in the fluid equations, simplifies the remainder of the derivation.We emphasize that the electron equations are unchanged from those in the original kglobal model.As in the original model all species are assumed to move perpendicular to the magnetic field with the E × B drift velocity.As discussed in Sec.2.1, the parallel velocity of the electrons is ordered to be much larger than the perpendicular velocity so that the parallel inertia of electros is retained while the perpendicular inertia is discarded: With these approximations equation (A3), reduces to where b is a unit vector parallel to the local magnetic field and κ = b • ∇b is the curvature vector.Similarly, the parallel particle electron velocity becomes Again, Eqs.(A7) and (A8) are unchanged from the original kglobal model.Summing equations (A2), (A5), (A7) and (A8), using the quasi-neutrality of equation (A1), which eliminates E, we obtain where n i = n ip + n if is the total ion density.Terms proportional to the electron mass on the left-hand side (LHS) of Equation (A9) have been discarded as discussed in Sec.2.1.The derivation of the large-scale parallel electric field, E ∥ , which is important in calculating the motion of particle electrons (see Eq. ( 2)), follows exactly as in Arnold et al. (2019).The parallel components of the four momentum equations are summed, the time rate of change of the parallel current is discarded, and we obtain an expression for E ∥ , Note that all of the ion contributions have dropped out and that the result is unchanged from that in Equation 17, except that n i = n ip + n if .Now the expression for the second term on the LHS of equation (A9) can be transformed as follows: Using the parallel component of the momentum equation of fluid ions in Equation (A2), the continuity equation, and the chain rule, we obtain the expression: The stress tensor T and pressure tensor P can be linked via bb From the combination of Equations ( A9), ( 17), and (A12), plus some algebraic manipulations we can obtain the following equivalent expressions: the second of which is Equation (26) of Section 2.

B. ENERGY CONSERVATION
Energy conservation is most simply explored using the equations of motion for each of the four species in Equations (A2)-(A5) rather than Equation (26).Starting from these more primitive equations is valid because nothing is discarded in the manipulation of the ion equations and in the electron equations the only discarded terms are the inertia associated with the perpendicular motion and this was previously shown to lead a set of equations that conserve energy (Arnold et al. 2019).Each of the equations has the generic form where the subscript j denotes the jth species.Taking the dot product with v j and using the continuity equation, we obtain For a generic gyrotropic pressure tensor P j , Equation (B17) takes the form where W = 1 2 m j n j v 2 j + 1 2 is the energy density and is the energy flux.For electrons only the bulk flow kinetic energy associated with the parallel motion needs to be retained.The bulk flow kinetic energy associated with the electron perpendicular motion can be neglected.Finally, Equation (B18) for the total energy of the particle species also follows directly from the guiding center equation for the parallel momentum and µ conservation.Summing the energy equations of the four species and integrating over all space yields the energy equation where J = J if + J ip + J ef + J ep is the total current.From Faraday's law The final energy conservation law then takes the form The derivation of energy conservation from the four individual momentum equations (see Eq. ( B16)) is exact.However, in the total fluid ion momentum equation in Eq. ( 26), the perpendicular momenta in the electron fluid and particle stress tensors and comparable time-dependent inertial terms have been discarded.The discarded terms scale like m e /m i compared with those retained.Also, the parallel electric field in Faraday's law (Eq.( 15)) is discarded compared with E ⊥ since their ratio scales like ρ s /L with ρ s the ion sound Larmor radius.The consequence is that energy conservation in kglobal is not exact but violations are small.An alternative derivation of energy conservation that parallels that in the original kglobal formulation (Arnold et al. 2019), can be obtained by defining a center-of-mass (CM) velocity n i v icm = n ip v ip + n if v if .The equation for the CM velocity is obtained by summing the two ion momentum equations (Eqs.(A2) and (A5)), where the CM pressure tensor is given by where the velocity dependent terms produce an effective presssure associated with the differences between the two ion velocities and the CM velocity.This equation for the CM velocity is then summed with the two electron momentum equations (Eqs.(A3) and (A4)) to produce B26)

Figure 1 .
Figure 1.The observed phase speed of the Alfvén wave, Cp, plotted against the anisotropy parameter α.The blue dots are the results of simulations.The intersection of the vertical and horizontal dotted lines marks the isotropic Alfvén wave.

Figure 2 .
Figure 2. Normalized growth rate of the firehose mode, γτA, versus mode number, m = 2πkL, over a span of unstable m values.Red crosses show values from simulations with equal electron and ion pressure anisotropies; values for unequal anisotropies are shown as green squares.The blue line is the theoretical growth rate.
OF THE FLUID ION MOMENTUM EQUATION