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 non-thermal 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.

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 non-thermal particle energization.8,9 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.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 experiments14 and space weather.15 

Many computational simulations, most notably those performed with particle-in-cell (PIC) and magnetohydrodynamic (MHD)16–18 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,19–24 but the spatial scales associated with solar flares can exceed 104 km, while the Debye length—a measure of the scale of kinetic effects—can be of the order of centimeters, a scale separation of 1010. 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 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 self-consistently 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 large-scale astrophysical reconnection.26,27 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.28 

However, magnetic reconnection can simultaneously produce non-thermal electrons and non-thermal protons.3,29 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 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.

For completeness, we first outline the governing equations for the original kglobal model before turning to the additions made in the upgraded version.

The original model includes three distinct plasma species: the ion fluid (number density ni), the electron fluid (nef), and particle electrons (nep).26,27 No conversion between species is allowed. The particle electrons are governed by the guiding-center equations given in Ref. 30. Their perpendicular velocity is determined by conservation of the magnetic moment:
μep=pep,22B=const.,
(1)
where pep 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,
dpep,dt=pep,vE·κμeγeb·BeE,
(2)
with vE 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,
nit=·(nivi),
(3)
and a modified momentum equation,
mi(nivi)t=·Ti(·Tef)+J×B/c(·Tep)+eniEb,
(4)
where the ion stress tensor with inertial contributions is
Ti=PiI+minivivi,
(5)
with Pi the ion scalar pressure and I the unit tensor. The corresponding electron fluid tensor is
Tef=PefI+menefγefvef,2bb,
(6)
with Pef the electron fluid scalar pressure. Finally, Tep is the particle electron gyrotropic stress tensor, including inertial contributions,
Tep=Tepbb+Pep(Ibb),
(7)
where Tep is the component of the stress tensor along the magnetic field B, and Pep is the usual perpendicular pressure,
Pep=dpepe22meγef,
(8)
where in the frame drifting with vE=cE×B/B2 there are no perpendicular flows, so f=f(x,pe,pe,t). Tep includes the mean parallel drifts of the particle electrons,
Tep=dpepe2meγef,
(9)
with pe 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, CAe=mi/meCA, where the Alfvén speed is defined as CA=B2/(4πmini). In this ordering, the electron momentum scales as meCAe=me/mimiCA and so is much smaller than that of the ion fluid. In contrast, the electron pressure is meCAe2=miCA2, 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
ddt(Piniγ)=0,
(10)
where the adiabatic index γ is taken to be 5/3.
The number density of the electron fluid is given by
nef=ninep,
(11)
while the parallel motion follows from the requirement of vanishing parallel current,
nefvef,=nivi,nepvep,.
(12)
The zero current condition follows from the constraint on the current from Ampère's law,
VdCA=JneCABc4πneCALdiL1,
(13)
where Vd is the drift speed associated with the current, di 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,
ddt(Pefnefγ)=0.
(14)
The magnetic field is advanced via Faraday's law,
Bt=c×E,
(15)
with
E=1cvi×B.
(16)
Finally, E is produced by magnetic-field-aligned gradients in the electron pressure and is given by27 
E=1e(nep+nef)(b··Tep+b·Tef).
(17)
Thus, Eqs. (1)–(17) constitute the complete set of equations for the kglobal model with particle electrons.

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
tmi(nifvif+nipvip)=S,
(18)
with vif and nif the fluid ion velocity and density, vip and nip 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
tminivif+tminip(vipvif)b=S,
(19)
where ni=nif+nip. 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  Appendixes A and  B.
Most of the other equations in Sec. II A remain unchanged, aside from minor notation adjustments (e.g., ninif+nip). Specifically, Eqs. (1) and (2) governing the motion of the particle electrons are unchanged. Equations (3) and (10) become equations for the ion fluid density,
nift=·(nifvif),
(20)
and ion fluid pressure,
ddt(Pifnifγ)=0.
(21)
The guiding-center equations for the particle ions are essentially identical to those of the electrons in Eqs. (1) and (2),
μip=pip,22B=const.
(22)
and
dpip,dt=pip,vE·κμib·B+eE.
(23)
The ions are taken to be non-relativistic. Including the new species also requires changes in the expressions for the electron fluid density,
nef=nif+nipnep,
(24)
and parallel motion,
nefvef,=nifvif,+nipvip,nepvep,.
(25)
The adiabatic equation of state for the electron fluid, Eq. (14), remains the same as do Eqs. (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 models, electron inertia can be discarded compared with that of the ions in the ion momentum equation, as discussed in Sec. II A. 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 Ref. 26. However, as shown in Eq. (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
mi(nivif)t=·Tif(·Tef)+J×B/c·Tip(·Tep)+eniEbIi,
(26)
where Ii is the difference in the parallel inertia between the fluid and particle ions, given by
Ii=tminip(vipvif)b=bmivif(nipnif·nifvifb·nipvipbnifvi,·nipnif)+nipnifbb··Tifbb··Tip+minipvipvifB(Ibb)·Bt,
(27)
where ni=nip+nif, and the ion fluid stress tensor, Tif, is defined in Eq. (5). The particle ion pressure tensor is given by
Tip=Pipbb+Pip(Ibb)+minipvipvip,
(28)
with the parallel and perpendicular pressures calculated in the frame of the bulk motion of the particle ions,
Pip=dpipi22mif
(29)
and
Pip=dpipi2mif,
(30)
with pip the ion momentum. Note that in the inertia term on the left side of Eq. (26), ni is the total ion density and not just the fluid ion density nif. In the limit nip0,Ii0 and Tip0 and Eq. (26) reduces to Eq. (4).

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 center-of-mass velocity in Eq. (B11) and calculate the ion fluid velocity by subtracting the ion particle flux.

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 (PP), a problem with a well-established 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πmini0). Lengths and times are normalized to an arbitrary macroscale length, L, and the Alfvén crossing time, τA=L/CA, respectively. Electric fields and temperatures are normalized to CAB0/c and miCA2, respectively. The kglobal model employs second-order diffusivities and fourth-order hyperviscosities (ν4) in each fluid equation to mitigate noise at the grid scale.26 

The simulations were performed within a two-dimensional (2D) spatial domain, coupled with a three-dimensional velocity space, featuring an initially uniform magnetic field of magnitude B0 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 nip/nif=nep/nef=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 simulation box (kx=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 Cp of an Alfvén wave, which is given by
Cp=CA14π(PP)B2αCA,
(31)
where α=14π(PP)/B2, which aligns with the Chew–Goldberger–Low (CGL) framework in the linear regime where pressure perturbations are negligible.33 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 Cp, normalized to CA, 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.

FIG. 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.

FIG. 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.

Close modal

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×105. The growth rate, which includes the pressure anisotropy drive and viscous damping, is given by γ=kCA|α|νk4. 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.34 All waves in the kglobal model are below the cyclotron frequency.

FIG. 2.

Normalized growth rate of the firehose mode, γτA, vs 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.

FIG. 2.

Normalized growth rate of the firehose mode, γτA, vs 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.

Close modal

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 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.

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available from the corresponding author upon reasonable 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 nif, nef, nep, and nip, respectively, that satisfy the charge neutrality constraint
nif+nipnefnep=0,
(A1)
where nif is advanced in time with a continuity equation. The particle motion is followed by the E×B perpendicular to B and parallel dynamics at their local velocity along B. Their number density, nip and nep, is 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 nef. In addition, each species satisfies a momentum equation:
(minifvif)t+·(minifvifvif)=enifE+Jif×B/cPif
(A2)
for fluid ions,
(menefvef)t+·(menefvefvef)=enefE+Jef×B/cPef
(A3)
for fluid electrons,
(menepvep)t=enepE+Jep×B/c·Tep
(A4)
for particle electrons, and
(minipvip)t=enipE+Jip×B/c·Tip
(A5)
for particle ions. In these equations,
  • 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.

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. II A, the parallel velocity of the electrons is ordered to be much larger than the perpendicular velocity so that the parallel inertia of electrons is retained while the perpendicular inertia is discarded:
vef,vep,vep,vef,.
(A6)
With these approximations, Eq. (A3) reduces to
(menefvefb)t=enefE+Jef×BcPefbB·menefvef,2Bmenefvef,2κ,
(A7)
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
(menepvep,b)t=enepE+Jep×B/c·Tep.
(A8)
Again, Eqs. (A7) and (A8) are unchanged from the original kglobal model. Summing Eqs. (A2), (A5), (A7), and (A8), using the quasi-neutrality of Eq. (A1), which eliminates E, we obtain
mi(nivif)t+mi(nip(vip,vif,)b)t+mie·(Jifvif)=J×B/cPifPefbB·menefvef,2Bmenefvef,2κ·(Tep+Tip),
(A9)
where ni=nip+nif is the total ion density. Terms proportional to the electron mass on the left-hand side (LHS) of Eq. (A9) have been discarded as discussed in Sec. II A.
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 Ref. 27. 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,
eniE=b··Tepb·PefB·menefvef,2B.
(A10)
Note that all of the ion contributions have dropped out and that the result is unchanged from that in Eq. (17), except that ni=nip+nif.
Now, the expression for the second term on the LHS of Eq. (A9) can be transformed as follows:
tnip(vipvif)b=b(vipvif)nipt+nipvipvifBBt+Bnipt(vipvifB).
(A11)
Using the parallel component of the momentum equation of fluid ions in Eq. (A2), the continuity equation, and the chain rule, we obtain the expression
tnip(vipvif)b=b(vipvif)nipt+nipvipvifB(Ibb)·Bt+nipb(b·(vip·vip)+b·(vif·vif)b··Pipminip+b·Pifminif).
(A12)
The stress tensor T and pressure tensor P can be linked via
bb··T=bb··P+bb··mnvv=bb··P+bb·(mnv·vmvnt).
(A13)
From the combination of Eqs. (A9), (17), and (A12), plus some algebraic manipulations, we can obtain the following equivalent expressions:
mieJift=mie·(Jifvif)Pif+nifni(J×B/cPef(·Tep)menefv,ef2κ+eniEb)+(Ibb)·[nipni(Pif+mee(Jif·)vif)nifni(·Pip+mee(Jip·)vip)minipnifniv,ipv,ifBBt]
(A14)
and
mi(nivif)t=mi·(nifvifvif)Pif+J×B/cPef(·Tep)menefvef,2κ+eniEb+bmivif(nipnif·nifvifb·nipvipbnifvi,·nipnif)minipvipvifB(Ibb)·Btnipnifbb··Tif(·Tip),
(A15)
the second of which is Eq. (26) of Sec. II.
Energy conservation is most simply explored using the equations of motion for each of the four species in Eqs. (A2)–(A5) rather than Eq. (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 to a set of equations that conserve energy.27 Each of the equations has the generic form
(mjnjvj)t+·(mjnjvjvj)=qjnjE+Jj×B/c·Pj,
(B1)
where the subscript j denotes the jth species. Taking the dot product with vj and using the continuity equation, we obtain
t(mjnjvj2/2)+·(mjnjvj2vj/2)=Jj·Evj·(·Pj).
(B2)
For a generic gyrotropic pressure tensor Pj, Eq. (B2) takes the form
Wjt+·Qj=Jj·E,
(B3)
where
W=12mjnjvj2+12P+P
(B4)
is the energy density and
Q=12mjnjvj2vj+(32Pj+Pj)vjb+(12Pj+2Pj)vj
(B5)
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, Eq. (B3) 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
dxt(Wif+Wip+Wef+Wep)=dxJ·E,
(B6)
where J=Jif+Jip+Jef+Jep is the total current. From Faraday's law,
tB28π+c4π·(E×B)+J·E=0.
(B7)
The final energy conservation law then takes the form
dxt(Wif+Wip+Wef+Wep+B2/8π)=0.
(B8)
The derivation of energy conservation from the four individual momentum equations [see Eq. (B1)] 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 me/mi compared with those retained. Moreover, 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 formulation27 can be obtained by defining a center-of-mass (CM) velocity nivicm=nipvip+nifvif. The equation for the CM velocity is obtained by summing the two ion momentum equations [Eqs. (A2) and (A5)],
t(minivicm)+·(minivicmvicm)=eniE+1cnievicm×B·Picm,
(B9)
where the CM pressure tensor is given by
Picm=mi[nip(vipvicm)(vipvicm)+(nif(vifvicm)(vifvicm)]+Pip+IPif,
(B10)
where the velocity dependent terms produce an effective pressure 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
t(minivicm)+·(minivicmvicm)=1cJ×B·Picm(·Tef)(·Tep)+eniEb.
(B11)
The energy equation in CM coordinates then follows by taking the dot product of this equation with vicm,
t(12minivicm2)+·(12minivicm2vicm)=J·Evicm··Picmvicm·[(·Tef)(·Tep)]+enivicmE.
(B12)
The contribution from the CM pressure tensor is linked to the ion thermal energy in the CM frame, WicmP, as in an MHD system. WicmP is given by
WicmP=12mi[nip(vipvicm)2+nif(vifvicm)2]+12Pip+Pip+32Pif=12mi[nipvip2+nifvif2nivicm2]+12Pip+Pip+32Pif.
(B13)
The time derivative of WicmP then follows from the time derivative of the individual kinetic energies as in Eq. (B2),
WicmPt+·(Qip+Qif)·(12nivicm2vicm)=vicm·Picm.
(B14)
This equation is then combined with Eq. (B12) to eliminate the Picm term,
t(12minivicm2+WicmP)+·(Qip+Qif)=J·EJef·EJep·E,
(B15)
where we have substituted nievicm=J+nefevef+nepevep and used vicm=vep=vef since the Hall terms in Ohm's law have been ordered out. The contributions from vef·(Tef) and vep·(Tep) were calculated from their respective momentum equations with their perpendicular bulk flow inertia discarded. Since
12minivicm2+WicmP=Wip+Wif,
(B16)
the integral of Eq. (B15) yields Eq. (B6) and the final energy conservation law in Eq. (B8) follows as calculated previously.
1.
S.
Tsuneta
, “
Structure and dynamics of magnetic reconnection in a solar flare
,”
Astrophys. J.
456
,
840
(
1996
).
2.
M.
Yamada
,
H.
Ji
,
S.
Hsu
,
T.
Carter
,
R.
Kulsrud
,
N.
Bretz
,
F.
Jobes
,
Y.
Ono
, and
F.
Perkins
, “
Study of driven magnetic reconnection in a laboratory plasma
,”
Phys. Plasmas
4
,
1936
(
1997
).
3.
R. P.
Lin
,
S.
Krucker
,
G. J.
Hurford
,
D.
Smith
,
H.
Hudson
,
G.
Holman
,
R.
Schwartz
,
B.
Dennis
,
G.
Share
,
R.
Murphy
et al, “
RHESSI observations of particle acceleration and energy release in an intense solar gamma-ray line flare
,”
Astrophys. J.
595
,
L69
(
2003
).
4.
J. T.
Gosling
,
R. M.
Skoug
,
D. J.
McComas
, and
C. W.
Smith
, “
Direct evidence for magnetic reconnection in the solar wind near 1 AU
,”
J. Geophys. Res.
110
,
A01107
, https://doi.org/10.1029/2004JA010809 (
2005
).
5.
T. D.
Phan
,
J. T.
Gosling
,
M. S.
Davis
,
R. M.
Skoug
,
M.
Øieroset
,
R. P.
Lin
,
R. P.
Lepping
,
D. J.
McComas
,
C. W.
Smith
,
H.
Reme
, and
A.
Balogh
, “
A magnetic reconnection X-line extending more than 390 Earth radii in the solar wind
,”
Nature
439
,
175
(
2006
).
6.
J.
Burch
,
R.
Torbert
,
T.
Phan
,
L.-J.
Chen
,
T.
Moore
,
R.
Ergun
,
J.
Eastwood
,
D.
Gershman
,
P.
Cassak
,
M.
Argall
et al, “
Electron-scale measurements of magnetic reconnection in space
,”
Science
352
,
aaf2939
(
2016
).
7.
D. E.
Gary
,
B.
Chen
,
B. R.
Dennis
,
G. D.
Fleishman
,
G. J.
Hurford
,
S.
Krucker
,
J. M.
McTiernan
,
G. M.
Nita
,
A. Y.
Shih
,
S. M.
White
et al, “
Microwave and hard x-ray observations of the 2017 September 10 solar limb flare
,”
Astrophys. J.
863
,
83
(
2018
).
8.
E. N.
Parker
, “
Sweet's mechanism for merging magnetic fields in conducting fluids
,”
J. Geophys. Res.
62
,
509
, https://doi.org/10.1029/JZ062i004p00509 (
1957
).
9.
H. E.
Petschek
, “
Magnetic field annihilation
,” in
AAS/NASA Symposium on the Physics of Solar Flares
, edited by
W. N.
Ness
(
NASA
,
Washington, DC
,
1964
), p.
425
.
10.
S.
Krucker
,
H. S.
Hudson
,
S. M.
White
,
S.
Masuda
,
J.-P.
Wuelser
, and
R. P.
Lin
, “
Measurements of the coronal acceleration region of a solar flare
,”
Astrophys. J.
714
,
1108
(
2010
).
11.
M.
Oka
,
S.
Ishikawa
,
P.
Saint-Hilaire
,
S.
Krucker
, and
R.
Lin
, “
Kappa distribution model for hard x-ray coronal sources of solar flares
,”
Astrophys. J.
764
,
6
(
2013
).
12.
M.
Øieroset
,
R. P.
Lin
,
T. D.
Phan
,
D. E.
Larson
, and
S. D.
Bale
, “
Evidence for electron acceleration up to ∼300 keV in the magnetic reconnection diffusion region in the Earth's magnetotail
,”
Phys. Rev. Lett.
89
,
195001
(
2002
).
13.
R.
Ergun
,
N.
Ahmadi
,
L.
Kromyda
,
S.
Schwartz
,
A.
Chasapis
,
S.
Hoilijoki
,
F.
Wilder
,
J.
Stawarz
,
K.
Goodrich
,
D.
Turner
et al, “
Observations of particle acceleration in magnetic reconnection–driven turbulence
,”
Astrophys. J.
898
,
154
(
2020
).
14.
M.
Yamada
,
F.
Levinton
,
N.
Pomphrey
,
R.
Budny
,
J.
Manickam
, and
Y.
Nagayama
, “
Investigation of magnetic reconnection during a sawtooth crash in a high-temperature tokamak plasma
,”
Phys. Plasmas
1
,
3269
3276
(
1994
).
15.
V.
Angelopoulos
,
A.
Artemyev
,
T. D.
Phan
, and
Y.
Miyashita
, “
Near-Earth magnetotail reconnection powers space storms
,”
Nat. Phys.
16
,
317
321
(
2020
).
16.
D.
Biskamp
, “
Magnetic reconnection via current sheets
,”
Phys. Fluids
29
,
1520
(
1986
).
17.
J. T.
Karpen
,
S. K.
Antiochos
, and
C. R.
DeVore
, “
The mechanisms for the onset and explosive eruption of coronal mass ejections and eruptive flares
,”
Astrophys. J.
760
,
81
(
2012
).
18.
J. T.
Dahlin
,
S. K.
Antiochos
,
J.
Qiu
, and
C. R.
DeVore
, “
Variability of the reconnection guide field in solar flares
,”
Astrophys. J.
932
,
94
(
2022
).
19.
M. A.
Shay
,
J. F.
Drake
,
R. E.
Denton
, and
D.
Biskamp
, “
Structure of the dissipation region during collisionless magnetic reconnection
,”
J. Geophys. Res.
25
,
9165
, https://doi.org/10.1029/97JA03528 (
1998
).
20.
M.
Hesse
,
K.
Schindler
,
J.
Birn
, and
M.
Kuznetsova
, “
The diffusion region in collisionless magnetic reconnection
,”
Phys. Plasmas
6
,
1781
(
1999
).
21.
M. A.
Shay
,
J. F.
Drake
, and
M.
Swisdak
, “
Two-scale structure of the electron dissipation region during collisionless magnetic reconnection
,”
Phys. Rev. Lett.
99
,
155002
(
2007
).
22.
W.
Daughton
,
V.
Roytershteyn
,
H.
Karimbadi
,
L.
Yin
,
B. J.
Albright
,
B.
Bergen
, and
K. J.
Bowers
, “
Role of electron physics in the development of turbulent magnetic reconnection in collisionless plasmas
,”
Nat. Phys.
7
,
539
(
2011
).
23.
X.
Li
,
F.
Guo
,
H.
Li
,
A.
Stanier
, and
P.
Kilian
, “
Formation of power-law electron energy spectra in three-dimensional low-β magnetic reconnection
,”
Astrophys. J.
884
,
118
(
2019
).
24.
Q.
Zhang
,
F.
Guo
,
W.
Daughton
,
H.
Li
, and
X.
Li
, “
Efficient nonthermal ion and electron acceleration enabled by the flux-rope kink instability in 3D nonrelativistic magnetic reconnection
,”
Phys. Rev. Lett.
127
,
185101
(
2021
).
25.
M.
Onofri
,
H.
Isliker
, and
L.
Vlahos
, “
Stochastic acceleration in turbulent electric fields generated by 3D reconnection
,”
Phys. Rev. Lett.
96
,
151102
(
2006
).
26.
J. F.
Drake
,
H.
Arnold
,
M.
Swisdak
, and
J. T.
Dahlin
, “
A computational model for exploring particle acceleration during reconnection in macroscale systems
,”
Phys. Plasmas
26
,
012901
(
2019
).
27.
H.
Arnold
,
J.
Drake
,
M.
Swisdak
, and
J.
Dahlin
, “
Large-scale parallel electric fields and return currents in a global simulation model
,”
Phys. Plasmas
26
,
102903
(
2019
).
28.
H.
Arnold
,
J. F.
Drake
,
M.
Swisdak
,
F.
Guo
,
J. T.
Dahlin
,
B.
Chen
,
G.
Fleishman
,
L.
Glesener
,
E.
Kontar
,
T.
Phan
et al, “
Electron acceleration during macroscale magnetic reconnection
,”
Phys. Rev. Lett.
126
,
135101
(
2021
).
29.
A.
Emslie
,
B.
Dennis
,
A.
Shih
,
P.
Chamberlin
,
R.
Mewaldt
,
C.
Moore
,
G.
Share
,
A.
Vourlidas
, and
B.
Welsch
, “
Global energetics of thirty-eight large solar eruptive events
,”
Astrophys. J.
759
,
71
(
2012
).
30.
T. G.
Northrop
, “
Adiabatic charged-particle motion
,”
Rev. Geophys. Space Phys.
1
,
283
304
, https://doi.org/10.1029/RG001i003p00283 (
1963
).
31.
J. F.
Drake
,
M.
Swisdak
,
H.
Che
, and
M. A.
Shay
, “
Electron acceleration from contracting magnetic islands during reconnection
,”
Nature
443
,
553
(
2006
).
32.
J. F.
Drake
and
M.
Swisdak
, “
Ion heating and acceleration during magnetic reconnection relevant to the Corona
,”
Space Sci. Rev.
172
,
227
240
(
2012
).
33.
G. F.
Chew
,
M. L.
Goldberger
, and
F. E.
Low
, “
The Boltzmann equation and the one-fluid hydrodynamic equations in the absence of particle collisions
,”
Proc. R. Soc. A
236
,
112
(
1956
).
34.
C.
Tronci
,
E.
Tassi
,
E.
Camporeale
, and
P. J.
Morrison
, “
Hybrid vlasov-MHD models: Hamiltonian vs. non-Hamiltonian
,”
Plasma Phys. Controlled Fusion
56
,
095008
(
2014
).