Higher-order nonequilibrium term: Effective power density quantifying evolution towards or away from local thermodynamic equilibrium

A common approach to assess the nature of energy conversion in a classical fluid or plasma is to compare power densities of the various possible energy conversion mechanisms. A forefront research area is quantifying energy conversion for systems that are not in local thermodynamic equilibrium (LTE), as is common in a number of fluid and plasma systems. Here, we introduce the ``higher-order non-equilibrium term'' (HORNET) effective power density that quantifies the rate of change of departure of a phase space density from LTE. It has dimensions of power density, which allows for quantitative comparisons with standard power densities. We employ particle-in-cell simulations to calculate HORNET during two processes, namely magnetic reconnection and decaying kinetic turbulence in collisionless magnetized plasmas, that inherently produce non-LTE effects. We investigate the spatial variation of HORNET and the time evolution of its spatial average. By comparing HORNET with power densities describing changes to the internal energy (pressure dilatation, $\rm{Pi-D}$, and divergence of the vector heat flux density), we find that HORNET can be a significant fraction of these other measures (8\% and 35\% for electrons and ions, respectively, for reconnection; up to 67\% for both electrons and ions for turbulence), meaning evolution of the system towards or away from LTE can be dynamically important. Applications to numerous plasma phenomena are discussed.


I. INTRODUCTION
The classical laws of thermodynamics quantify the conversion of energy for processes occurring in systems in local thermodynamic equilibrium (LTE).However, many systems of interest are not in LTE.The approaches by Chapman and Enskog [1], Onsager [2], Grad [3] and its extensions [4], and the gyrokinetic approach [5,6] have revealed many important insights about energy conversion for systems not in LTE but the models assume the departures from LTE are small.Many systems are routinely far from LTE, such as plasmas in space and astrophysical settings where densities are low and temperatures are high, leading to low collisionality.This indicates there is a necessity for a non-perturbative approach to quantify non-LTE energy conversion.
Another recent approach that allows for the identification of different non-LTE energy conversion mechanisms in a wide variety of plasma systems is the so-called "fieldparticle correlation" approach [23].It directly identifies the energy conversion between electromagnetic fields and charged particles in a plasma as a function of the phase space density.It has been applied to many plasma systems in simulations and observations [24][25][26][27][28][29][30][31].
In instances where a system deviates significantly from LTE, each of the infinite number of moments of the phase space density can evolve and become dynamically significant.As internal energy is related to only a single moment of the phase space density, a description of the evolution of the other moments can become pertinent.However, a closed set of evolution equations for each of the infinite moments is unattainable due to the closure problem.One potential strategy to describing every moment is to employ the Boltzmann (or Vlasov) transport equation for the phase space density itself, but this is challenging since the phase space density is a sixdimensional (plus time) field.Notably, recent research [32] has demonstrated that the evolution of all moments associated with a non-LTE system can be described by the relative entropy (that is a three-dimensional plus time field).It is a first-principle derivation from the Boltzmann equation using an entropic approach [33].
Mechanisms that can change the internal energy (pressure dilatation, Pi − D, and divergence of vector heat flux density), are typically quantified as power densities, i.e., rates at which the energy density changes.In addition, J σ • E describes the rate of conversion of energy density between electric fields and charged particles of species σ, where J σ is the current density of particle species σ and E is the electric field.Previous studies [10,20] have demonstrated that J σ • E does not directly result in a change in internal energy and can only change the bulk kinetic energy.While questions persist regarding the precise mechanism through which electric fields can induce changes in internal energy, we provide only a cursory treatment of this energy conversion channel because it does not directly go to internal energy and has been compared to other internal energy measures previously.
Here, we use the result of Ref. [32] to define a quantity P σ,HORNET describing the rate of change of the departure of a phase space density from LTE that has dimensions of power density, where HORNET stands for "higher-order non-equilibrium term."We call it an effective power density because it is not a power density in the sense that it does not describe a rate of change of energy density.In this study, we calculate P σ,HORNET in particle-in-cell (PIC) simulations of magnetic reconnection and decaying plasma turbulence.By locally comparing P σ,HORNET to other power densities, we can ascertain the relative prevalence of the evolution of a phase space density towards or away from LTE compared to the energy conversion taking place.We discover that the spatially averaged HORNET power density can be a substantial percentage of the net power density changing the internal energy at the times of interest; for a simulation of reconnection, we find that this percentage in the electron diffusion region is 8% and 35% for electrons and ions, respectively; in a simulation of turbulence, this percentage for the entire domain is 67% for both electrons and ions.We expect P σ,HORNET will be useful for many non-LTE processes in plasma physics and other fields of science.
The layout of this paper is as follows.Sec.II contains a review of the relevant results from Ref. [32] and introduces P σ,HORNET .Sec. III describes the numerical simu-lation setup.Secs.IV and V give the simulation results for a simulation of magnetic reconnection and decaying plasma turbulence, respectively.Sec.VI summarizes the results and discusses applications and implications.

A. Review
The non-relativistic Boltzmann transport equation [34] for species σ is given by where f σ is the phase space density, F σ is the net body force, m σ is the mass of the constituents, r, v, and t are the position-space, velocity-space, and time coordinates of phase space, ∇ and ∇ v are the position-and velocityspace gradient operators, respectively, and C[f ] is a collision operator that may be a sum of intra-species and inter-species collision operators given by C and integrating over all velocity space, where v ′ σ = v−u σ is the random velocity, u σ = (1/n σ ) f σ vd 3 v is the bulk flow velocity and n σ = f σ d 3 v is the number density, one obtains the time evolution equation of internal energy density Ēσ where P σ = n σ k B T σ is the pressure tensor, T σ is the temperature tensor with elements is the volumetric heating rate per particle due to interspecies collisions.
The pressure-strain interaction is defined as Note, Ēσ,int = 3P σ /2 = (3/2)n σ k B T σ , where we call P σ = tr(P σ )/3 the effective pressure and T σ = tr(T σ )/3 the effective temperature, and tr is the trace.The lack of the J σ • E term in this equation is used as evidence that the term does not directly change the internal energy density.
In Ref. [32], it was shown using an entropic approach that for systems not in LTE with a fixed total number of particles N σ of each species, the so-called first law of kinetic theory is given by where dE σ,gen and dQ σ,gen are increments in the so-called generalized energy per particle and generalized heat per particle, respectively, and Equation (5a) describes the rate of change of compressional work per particle dW σ /dt done by the system.The two terms in Eq. (5b) describe the rate of change of internal energy per particle dE σ,int /dt and the rate of change of effective energy per particle (the relative energy per particle) describing the local evolution towards or away from LTE dE σ,rel /dt, respectively.The two terms in Eq. (5c) describe the heating rate per particle associated with changes to the internal energy )] and an effective rate of heat per particle associated with the local phase space density not being in LTE dQ σ,rel /dt, respectively.Equation (5d) describes the rate of energy change per particle associated with collisions.In Eq. (5d), ∆ 3 r σ ∆ 3 v σ is the six-dimensional phase space volume element [33,35,36].Quantitatively, the two relative terms associated with the departure of the phase space density from LTE are given by where the relative entropy per particle [37] and the thermal relative entropy density flux divergence per particle [32] are defined as Here, f σM is the "Maxwellianized" phase space density associated with f σ given by [37] with n σ , u σ , and T σ obtained from f σ .Because s σv,rel /n σ identically vanishes if f σ = f σM [37], so relative entropy per particle contains information about the non-Maxwellianity of a local phase space density, and its Lagrangian time derivative describes how the Maxwellianity locally changes in time.

B. The Higher Order Non-Equilibrium Term
Power Density Pσ,HORNET Here, we construct a quantity describing the evolution of a phase space density towards or away from LTE with dimensions of power density to be on the same footing as the pressure-strain interaction and the divergence of the vector heat flux density.We note from Eq. ( 3) that these quantities give the time rate of change in the Lagrangian frame of the internal energy per particle scaled by the number density n σ .We, therefore, argue the analogous effective power density for describing the evolution of a phase space density towards or away from LTE is −n σ dE σ,rel /dt.We call this the HORNET effective power density P σ,HORNET : Since a Maxwellian phase space density is the maximum entropy state for a fixed number of particles and internal energy [34], a positive P σ,HORNET indicates that the phase space density is evolving away from Maxwellianity in the next instant of time in the Lagrangian reference frame and a negative P σ,HORNET indicates that the phase space density is evolving towards Maxwellianity in the next instant of time, i.e., it thermalizes.The magnitude quantifies the rate that the shape of the phase space density changes scaled by the effective pressure P σ .
To isolate the portion of Eq. ( 4) associated with the departure of the phase space density from LTE, we derive an equation for the time evolution of dE σ,rel (see also Ref. [38]).Substituting Eqs.(5a) -(5c) into Eq.( 4) and using Eq.(3) gives Assuming collisions conserve particle number and are elastic, the collision operator satisfies the relations 0 from conservation of particle number, m σ vC σ d 3 v = 0 from conservation of momentum for intra-species collisions, and (1/2)m σ v 2 C σ d 3 v = 0 from conservation of energy for intra-species collisions.A brief derivation reveals that Qσ,coll,inter is equivalent to where f σM is defined in Eq. (8).Then, Eq. (10) becomes where Eq. ( 12) is a special case of Eq. (3.13) from Ref. [38].
It may seem as though Eq. ( 12) implies that the relative terms are decoupled from the terms related to the internal energy [see Eq. ( 3)], but this is not the case [32] because the thermodynamic heat dQ σ /dt terms depend on non-Maxwellian features of the phase space density, namely Π σ and q σ , and their evolution is included in dE σ,rel /dt.Thus, Eq. (12) shows that the HORNET effective power density describes evolution towards or away from LTE that occurs as a result of relative heat (non-LTE features that change the shape of the phase space density) and/or collisions.
To better show that the relative entropy is associated with changes in shape of the phase space density, we rewrite s σv,rel /n σ from Eq. (7a) as where we define the distribution function as f ′ σ = f σ /n σ .Here, f ′ σ is normalized to 1, while f σ is normalized to n σ .By definition, f ′ σ and f ′ σM are independent of n σ , i.e., they are both unchanged if f σ changes in such a way that its shape does not change even if the density does change.Thus s σv,rel /n σ is independent of changes solely in n σ .Since the HORNET effective power density is proportional to d(s σv,rel /n σ )/dt, it is zero for any process that only changes the density.

III. NUMERICAL SIMULATIONS
Numerical simulations are performed using p3d [39], a massively parallel PIC code.All simulations are 3D in velocity-space, and 2.5D in position-space, i.e., vectors have three components and there is one invariant spatial dimension.Macro-particles are stepped forward using a relativistic Boris particle stepper [40] and electromagnetic fields are stepped using the trapezoidal leapfrog method [41].To enforce Poisson's equation, p3d utilizes the multigrid method [42] to clean the electric field.We employ periodic boundary conditions in both spatial directions for all simulations.
Simulated quantities are presented in normalized units.Lengths are normalized to the ion inertial scale d i0 = c/ω pi0 , where c is the speed of light, ω pi0 = (4πn 0 q 2 i /m i ) 1/2 is the ion plasma frequency based on a reference number density n 0 and q i and m i are the ion charge and mass, respectively.Magnetic fields are normalized to a reference magnetic field B 0 .Velocities are normalized to the Alfvén speed c A0 = B 0 /(4πm i n 0 ) 1/2 .Times are normalized to the inverse ion cyclotron frequency Ω −1 ci0 = (q i B 0 /m i c) −1 .Temperatures are normalized to m i c 2 A0 /k B , vector heat flux densities are in units of (B 2 0 /4π)c A0 , and power densities are in units of (B 2 0 /4π)Ω ci0 .Unrealistic values of the speed of light c = 15 and electron-to-ion mass ratio m e /m i = 0.04 are employed for numerical expedience, but we expect these values do not qualitatively change the results presented here.To calculate kinetic entropy, we employ the implementation from Ref. [36] and optimize the velocity-space grid scale by checking the agreement between the simulated kinetic entropy density s σ for various ∆v σ (where σ = e or i for electrons or ions) and the theoretical value at t = 0 [43].
To reduce the effect of PIC noise in the 2D plots that follow, we implement a recursive smoothing procedure on the raw simulation data.We smooth over a width of four cells four times, followed by taking spatial or temporal derivatives, and then recursively smooth over four cells four times again.We determine the appropriateness of this smoothing procedure by testing various smoothing widths and iterations.Our selection is based on evaluating how spatial variations change with different smoothing parameters, ensuring that clearer, sharper structures are observed in 2D plots and 1D cuts through the domain without the occurrence of oversmoothing (not shown).We ascertain that using six or more iterations results in oversmoothing for the simulations in the present study, particularly in the vicinity of the X-line at y = y 0 in the reconnection simulation.

A. Magnetic Reconnection Simulations
The reconnection simulation system size is L x × L y = 12.8 × 6.4 with 1024 × 512 grid cells that are initialized with 6,400 weighted particles per grid (PPG).The use of weighted particles improves load imbalance by making regions of high density have a similar number of particles for the processors to advance in time as in low density regions.Our choice of the simulation system size ensures that boundary effects do not interfere with the physics near the electron diffusion region (EDR) where we focus our interest.
We initialize the simulation using the standard double tanh magnetic field profile B x (y) = tanh [(y − L y /4)/w 0 ] − tanh [(y − 3L y /4)/w 0 ] − 1, with no initial out-of-plane (guide) magnetic field, where w 0 = 0.5 is the initial half-thickness of the current sheet.Therefore, the reference magnetic field B 0 is the initial asymptotic magnetic field strength.At initialization, the electron and ion density profiles are n(y) = {1/[2(T e + T i )]}{sech 2 [(y − L y /4)/w 0 ] + sech 2 [(y − 3L y /4)/w 0 ]} + n up , where the initial asymptotic upstream plasma density n up has the value of 0.2, and the difference between the peak current sheet number density and the upstream number density is the reference density n 0 .The electron temperature T e and ion temperature T i are uniform and T e = 1/12, T i = 5T e at initialization, which is inspired by the temperature ratio in Earth's magnetotail.Both species are loaded as drifting Maxwellian velocity distribution functions with electron and ion out-of-plane drift speeds satisfying u e,z /T e = −u i,z /T i initially.We seed an Xand O-line pair in each of the two current sheets to initiate reconnection using a magnetic field perturbation of the form δB x = −B pert sin (2πx/L x ) sin (4πy/L y ) and with B pert = 0.05.The rationale behind selecting this perturbation amplitude is that it reduces the computational cost of the simulation by quickening the transition of the reconnecting system into the non-linear growth phase which is when substantial deviations from Maxwellian behavior are expected to occur.This initial perturbation amplitude is not expected to modify the nonlinear phase of reconnection.
The electric field is cleaned every 10 particle time steps to ensure Poisson's equation is satisfied and improve energy conservation.The smallest length scale of the system is the electron Debye length λ De = 0.0176; we choose a grid-length ∆ = 0.0125.The smallest time scale of the system is the inverse of electron plasma frequency ω −1 pe = 0.012 and we choose the particle timestep ∆t = 0.001, with the field time step of half of this.These choices result in excellent total energy conservation, which we need due to the small signals of energy conversion as we discuss in Sec.IV; we see an increase in total energy by only 0.032% by t = 16.The plots we show are from the lower current sheet at time t = 13, during non-linear growth when the rate of reconnection is increasing most rapidly in time.
For the calculation of entropies, the range of the velocity space in each direction is restricted to [−v lim,σ , v lim,σ ] [33], where v lim,σ is chosen to be less than c but contain all the particles.We use v lim,e = 14.2, and v lim,i = 7.8.Each velocity space direction is discretized using n bin,σ bins, where n bin,e = 26, and n bin,i = 34.This gives the optimal velocity space grid scales ∆v σ , which are ∆v e = 1.092 and ∆v i = 0.459 using the procedure described in Refs.[33,35,36,43].The agreement with the theoretical values at t = 0 are within ±0.7% and ±1.5% in the far upstream region for electrons and ions, respectively, and within ±2% and ±3% at the center of the current sheet for electrons and ions, respectively.By t = 16, the total entropy is conserved to 2.08% and 2.29% for electrons and ions, respectively.
In Appendices A and B, we plot electron reduced phase space densities (with one velocity dimension integrated out).The data are from a different simulation with the same parameters used in this study, except for PPG=25,600 instead of 6,400 [32].We use a spatial domain of size 0.0625 × 0.0625 centered around the location for which phase space densities are plotted and bin the particles with a velocity space bin of size 0.1 in all three velocity directions.Reduced phase space densities are in units of n 0 /c 2 A0 .

B. Turbulence Simulations
The turbulence simulation uses a periodic domain of size L x × L y = 37.3912 × 37.3912 with 1024 x 1024 grid cells initialized with 6,400 unweighted PPG.Using weighted particles is unnecessary for this simulation because the density is initially uniform.There is an initial uniform magnetic field that has a strength given by the reference value B 0 = 1 and points along ẑ.The velocity and magnetic field fluctuations are Alfvènic with random phase and are excited in a band of wave numbers k ∈ [2, 4] × 2π/L x with a flat spectrum.The initial root mean square fluctuation amplitude for both velocity and magnetic field are set to 0.25.The initial cross helicity is very small.These initial conditions are analogous to those in a larger domain in Ref. [44].
The initial number density is uniform with a value given by the reference density n 0 = 1.The electron temperature T e and ion temperature T i are initially uniform and equal to 0.3.Both electrons and ions are initially drifting Maxwellian velocity distribution functions at the initial temperature and drifting with the local bulk flow speed.The electric field is cleaned every 40 particle time steps.The cadence of electric field cleaning is determined by performing four simulations with the cadence set to 20, 40, 50, and 100 particle time steps, respectively, and selecting the value that produced the best energy conservation.The smallest length scale of the system is the electron Debye length λ De ≃ 0.03651 and we choose a gridlength ∆ = λ De .The smallest time scale of the system is the inverse of electron plasma frequency ω −1 pe ≃ 0.0133; we choose the particle time-step to be ∆t = 0.005 with the field time step ∆t/3.This results in excellent total energy conservation of 0.013% by the final time of t = 50.The non-linear time, the time scale over which non-linear interactions between different wave modes become significant, is For calculating entropies, we obtain an optimal electron velocity space grid scale of ∆v e = 2.028 using v lim,e = 14.2, and n bin,e = 14.At t = 0, the choice for ∆v e gives agreement with the theoretical values within ±0.16% on average along horizontal and vertical cuts through the center of the simulation domain.For ions, using v lim,i = 5.7, and n bin,i = 50, we obtain ∆v i = 0.228 which gives agreement within ±0.12% of theoretical values on average along the same cuts at t = 0.By the final time of t = 50, the total electron and ion entropies are conserved to 0.28% and 0.55%, respectively.

C. Relative Entropy Implementation
In a previous study [32], the calculation of (d/dt)(s σv,rel /n σ ) was carried out by post-processing output data of s σ and other fluid moments that had been rounded.For the present study, we implement the quantities natively within p3d to reduce rounding errors.We find the overall 2D structures in (d/dt)(s σv,rel /n σ ) are similar when calculated in the two different ways (not shown).
To assess the importance of numerical error, we note that f σ at t = 0 in both the reconnection and turbulence simulations are loaded as drifting Maxwellians.Theoretically, from Eq. (7a), we would expect s σv,rel /n σ to vanish initially everywhere in the simulation domain.However, these initial Maxwellians have PIC noise due to a finite number of particles per grid, so s σv,rel /n σ is non-zero.We confirm that the error of non-zero s σv,rel /n σ values in our reconnection simulation is due to finite PPG by comparing electron data in simulations with different PPG.For the reconnection simulation with PPG = 6,400 used in this manuscript, the average value of s ev,rel /n e over the whole domain at t = 0 is ≃ −0.013.In a test simulation with PPG = 400, we find the average value of s ev,rel /n e at t = 0 is ≃ −0.08.Noise in PIC simulations scales with 1/ √ PPG, so we expect this decrease in PPG by a factor of 16 should increase the noise by a factor of 4. This is close to the factor by which we observe the decrease, suggesting that the values of relative entropy per particle calculated within our code is a result of using finite PPG.This provides evidence that our entropy implementation is valid at t = 0.This lays out an approach for future studies of relative entropy with PIC codes by measuring the desired quantity with a particular value of PPG and then increasing PPG to reduce the noise according to 1/ √ PPG until the signal is higher than the noise.

IV. RESULTS: MAGNETIC RECONNECTION
A. Spatial distribution and amplitude of HORNET Simulation results of P σ,HORNET at t = 13 are shown in Fig. 1, with electrons in (a) and ions in (e).Representative magnetic field line projections are in black, and we shift the coordinate system to be relative to the Xline location at (x 0 , y 0 ) = (9.59375,1.59375).The EDR is approximately located at We begin with a discussion of the spatial structure of P σ,HORNET .For the electrons, we note that Ref. [32] included a plot in their Fig.3(e) of dE e,rel /dt for a reconnection simulation that only differed from the present simulation by its PPG.Since P σ,HORNET = −n σ dE e,rel /dt from Eq. ( 9), the structure of P e,HORNET is largely similar to −dE e,rel /dt from Ref. [32].Consequently, we provide only an abbreviated analysis of the electrons here.
Upstream of the EDR, electrons become trapped in the bent magnetic field lines, leading to phase space densities that are elongated in velocity space in the direction parallel to the local magnetic field [45].This implies the shape of the phase space density is becoming less Maxwellian in the Lagrangian reference frame, which is associated with positive P e,HORNET ≃ 0.003.Near the X-line, triangular striated phase space densities develop due to electron Speiser orbits [46][47][48], a further evolution away from Maxwellianity, therefore associated with positive P e,HORNET ≃ 0.007.Interestingly, sandwiched between these two regions, there is a slim boundary extending from the two inflow regions (0.2 < |y − y 0 | < 0.3) to the two outflow regions (0.7 < |x − x 0 | < 1.2) inside the EDR in which P e,HORNET ≃ −0.0015 is weakly negative.Physically, this happens because the upstream phase space densities are elongated in the parallel direction [45], and positions close enough to the X-line have another population due to the Speiser orbits associated with motion perpendicular to the magnetic field.The resultant phase space density in the (v x , v y )-plane is closer to being a Maxwellian.Therefore, an electron fluid element entering the EDR has a negative P e,HORNET .Phase space densities in relevant regions are plotted in Appendix A. Thermalization of the outflowing electron jets downstream of the X-line [49][50][51][52] is observed at the EDR outflow edge (1.6 < |x − x 0 | < 2.0), associated with negative P e,HORNET ≃ −0.002.
For the ions, the ion diffusion region (IDR) is expected to have a vertical extent of |y − y 0 | < 2.24, but the simulation domain is too small to allow full ion coupling to the reconnection process; the simulation domain would have to be ≥ 40 × 2.24 ≃ 90 for the ions to fully couple [53].Upstream of the EDR, Fig. 1(e) demonstrates that P i,HORNET is weakly positive and gradually increases approaching the X-line, where it has a maximum value of ≃ 0.013.This indicates that ions become more non-Maxwellian (in the Lagrangian frame) approaching the X-line.This is consistent with expectations since ions demagnetize when their gyroradii exceed the distance from their guiding center to the magnetic field reversal, so only the most energetic ions are demagnetized further from the X-line and more ions are demagnetized closer to the X-line.
The ions in the exhaust just downstream of the EDR have positive P i,HORNET with characteristic value 0.006, indicating an evolution away from Maxwellianity in these regions.This is likely from demagnetizing ions crossing the separatrix rather than crossing through the EDR.We expect for a larger system and at later times that ions would thermalize further downstream from the EDR, which would be associated with P i,HORNET < 0, but we do not observe it in our simulation since the system is not in the steady-state and the simulation domain is small.Interestingly, starting at approximately |x − x 0 | ≃ 1.6 along the separatrices outside of the EDR, P i,HORNET is weakly negative with a value of ≃ −0.001.This implies the ion phase space density becomes more Maxwellian (in the Lagrangian frame) there.

B. Comparison to other power densities
We first perform a comparison of P σ,HORNET with the power densities describing changes to the internal energy per particle described in Eq. (3) at t = 13 in order to determine the relative importance of each as a function of position during the non-linear growth phase.For electrons, the pressure dilatation −P σ (∇ • u σ ) is shown in Fig. 1 Looking at ions in the vicinity of the X-line, we observe a strongly negative pressure dilatation (panel (f)), reaching values of approximately −0.06 which by itself would decrease the ion internal energy through expansion near the X-line.In contrast, Pi − D i is positive throughout the EDR and most of its surroundings, with a characteristic value near 0.02.Near the X-line, we find −∇ • q i is positive (the vector heat flux density is converging) [54] with a value ≃ 0.06 (panel (h)), which by itself would increase the ion internal energy.Panel (e) shows that P i,HORNET is positive in the vicinity of the X-line with a value of approximately 0.013.These results imply that, like the electrons, there is a local net increase near the X-line of the ion internal energy, and the sum of the three terms gives a net power of approximately 0.02, and the ion phase space density evolves away from Maxwellianity.For ions, the effective power going into changing the shape of the phase space density is approximately 65% of the net power due to pressure dilatation, Pi − D i , and vector heat flux divergence near the X-line.Near the downstream edge of the EDR, pressure dilatation is instead positive with value near 0.05, whereas Pi − D i is negligible.We also see that −∇ • q i ≃ −0.06, with P i,HORNET ≃ 0.006.Thus, contrary to what is observed for electrons, there is a local decrease in the ion internal energy and the ion phase space density is becoming more non-Maxwellian.
We note that −P i (∇ • u i ) and Pi − D i were plotted in a previous simulation study [21] using a different code and system parameters (see their Fig. 1 for comparison).Their data was taken during the steady-state phase while ours is during the non-linear growth phase.We observe some differences: (i) in our simulation, −P i (∇ • u i ) is notably more localized around the X-line in the x direction [see Fig. 1(f)], presumably because our system size is smaller (theirs is 25 ×25); and (ii) Pi − D i has a single positive peak signature at the X-line in our simulation, while it has a double-peaked structure in Ref. [21], presumably because theirs is in the steady state.
We also note the spatial distribution of −∇ • q σ for electrons and ions reveal structural similarities, except the ion structures are on larger scales.This is consistent with previous work that shows −∇ • q σ can be non-zero and dynamically significant locally [16,17], even though it vanishes when integrated over the whole periodic computational domain [11].Physically, a local non-zero vector heat flux density emerges from velocity-space asymmetries of phase space densities.We provide an explanation of the vector heat flux density structure and its divergence in Appendix B.
We now quantitatively compare the signed spatial averages of each of the previously discussed power densities, for electrons and ions separately, in order to compare their importance over a more extended region and to see how they evolve as a function of time in the simulation.We fix a box centered at the X-line of size 4 d i0 in length and 0.8 d i0 in width (representing the EDR size at t = 13), and then calculate the signed spatial average (denoted by ⟨⟩) of each power density in this region as a function of time.The result is shown in Fig. 2(a) and (b) for electrons and ions, respectively.In each, solid green denotes ⟨P σ,HORNET ⟩, dotted orange denotes ⟨−P σ (∇ • u σ )⟩, dashed purple denotes ⟨Pi − D σ ⟩, and ⟨−∇ • q σ ⟩ is represented by dash-dotted magenta.
At early time in both panels (t ≲ 3), there is an initial positive spike of pressure dilatation, followed by an overshoot where it goes negative.This occurs because the initial conditions of the simulation contain a single drifting Maxwellian population, instead of two populations as is needed for the Harris sheet which is an exact kinetic equilibrium.The system self-consistently adjusts as the simulation is evolved.After this, most of the power densities are small until t ≃ 12 when the reconnection rate begins to increase rapidly in time.Consequently, we fo-cus mostly on these later times.
For the electrons, panel (a) reveals that ⟨Pi − D e ⟩ and ⟨−P i (∇•u e )⟩ are the two largest terms, nearly equal early in the non-linear growth phase with ⟨Pi − D e ⟩ becoming larger in the late non-linear growth phase.⟨P e,HORNET ⟩ remains relatively small throughout the evolution.We know from Fig. 1(a) that P e,HORNET is non-zero locally, so this result suggests that the net positive and negative effective power density contributions within the EDR associated with evolution towards or away from LTE mostly cancels out.Focusing on the net power densities at t = 13 for electrons, we find ⟨Pi − D e ⟩ ≃ 0.008, ⟨−P e (∇ • u e )⟩ ≃ 0.006 and ⟨−∇ • q e ⟩ ≃ −0.002, giving a net power density to change the electron internal energy of 0.012.Since ⟨P e,HORNET ⟩ ≃ 0.001, we find ⟨P e,HORNET ⟩ ≃ 0.13 ⟨Pi − D e ⟩ ≃ 0.17 ⟨−P e (∇ • u e )⟩.Because of the smallness of the average values and potential uncertainties in the measurements of P σ,HORNET , these fractions could have significant uncertainties.The raw values indicate that ⟨P e,HORNET ⟩ represents a small but non-negligible percentage of about 8% of the rate of change of energy density going into changing the electron internal energy.
Interestingly, Fig. 2 shows that ⟨P e,HORNET ⟩ becomes negative at t = 15, suggesting that the net effect in the EDR is for the electron phase space density to evolve to be more Maxwellian at this time, in contrast to t = 13 where it becomes less Maxwellian.While this could be impacted by the level of uncertainty in the measurements of P e,HORNET , there is a physical reason to believe this could be valid.As seen in Fig. 1(a), the strongest negative P e,HORNET regions are at the downstream edges of the EDR where electrons thermalize.The electrons start to thermalize over a larger area as reconnection proceeds through its non-linear growth phase, so these patches of P e,HORNET < 0 grow with time (not shown), leading to the overall change in sign of ⟨P e,HORNET ⟩.
For the ions shown in panel (b) the dominant positive energy conversion channel is ⟨−∇ • q i ⟩, which is nearly balanced by a negative ⟨−P i (∇ • u i )⟩.There is a small net positive ⟨Pi − D i ⟩ late in the non-linear growth phase (t ≥ 12), that leads to a very weak increase in ion internal energy.The net effective energy conversion associated with the evolution towards or away from LTE, i.e., ⟨P i,HORNET ⟩, though positive, is small compared to the other ion terms.This is expected from the 2D plot of P i,HORNET (see Fig. 1(e)) that shows it is non-negative locally inside the EDR leading to a net evolution of ions towards LTE.These results are very sensible physically, as the system size is not large enough for the ions to fully couple to the reconnection process in the simulation, so their net energy conversion is quite weak.At t = 13, we find ⟨−∇ • q i ⟩ ≃ 0.015, ⟨Pi − D i ⟩ ≃ 0.01, and ⟨−P i (∇ • u i )⟩ ≃ −0.005, for a total of 0.02 going into ion internal energy.In comparison, ⟨P i,HORNET ⟩ ≃ 0.007 ≃ 70% ⟨Pi − D i ⟩ ≃ 50% ⟨−∇ • q i ⟩.Thus, in contrast to what we see for electrons, a larger percentage (35%) of the energy going into ion internal energy goes to the evolution of the ions towards LTE.
We briefly examine the comparison between P σ,HORNET and J σ • E.
At t = 13, we observe that the peak amplitude of J σ • E exceeds that of the largest energy conversion metric (i.e., −∇ • q σ ) by more than a factor of two (not shown).Moreover, it surpasses P σ,HORNET by more than an order of magnitude.This is expected because, in reconnection, a significant portion of the electromagnetic energy is converted into bulk kinetic energy of the particles, so the portion associated with internal energy and P σ,HORNET are smaller.In summary, the results of this section suggest P σ,HORNET is a useful tool for assessing the relative rates at which the phase space density evolves towards or away from LTE and for comparing the rate to standard power densities during the reconnection process.

V. RESULTS: DECAYING TURBULENCE A. Spatial distribution and amplitude of HORNET
We now discuss the structure of P σ,HORNET in the decaying turbulence simulation using 2D position-space plots of the various power densities in Fig. 3 at t = 30 ∼ 2 τ nl , with representative in-plane magnetic field line projections as dashed-black lines.This time is chosen as it is after the time at which mean square current density ⟨J 2 ⟩ peaks (plot not shown), The panels are of the same quantities and are in the same order as those plotted in Fig. 1 for the reconnection simulation.
We first consider electrons.In Fig. 3(a), we see that P e,HORNET has substantially non-zero signatures co-located with intermittent structures (current sheets).The regions with strongest P e,HORNET up to a magnitude of about 0.016 are centered around (x, y) ≃ (18,36), (21,2), (29,5), (15,18), (15,4), and (30,22).P e,HORNET can have either sign in and around these structures.We observe one location at (x, y) ≃ (15, 7) around a current sheet where the structure of P e,HORNET somewhat resembles that of the reconnection X-line shown in Fig. 1(a), with a weak positive value surrounded by a negative value.We do not see this feature at other current sheets in our turbulence simulation.The absence is possibly due to the presence of asymmetries [55] or flow shear arising at reconnection sites in the decaying turbulent system or the presence of the mean (guide) field in the turbulence simulation but not in the reconnection simulation that may modify the structure of P e,HORNET .Interestingly, we observe numerous bands with non-zero P e,HORNET coinciding with extended current sheets that are multiple ion-scales in length (∼ 10 d i ) and electronscale in thickness.
Panel (a) makes it appear that P e,HORNET is small away from the current sheets, but this is an artifact of the colorbar because the largest values are significantly higher than values in the rest of the domain.In Appendix C, we reproduce the same plot but with the color- bar saturated at the maximum value used for P i,HORNET to better observe the structures in P e,HORNET away from the current sheets.This plot makes it clear that there are regions with non-zero of P e,HORNET away from the current sheets.There are noticeable features in P e,HORNET within ion-scale eddies of size 5-10 d i .For example, in eddies centered around (x, y) ≃ (12,12), (29,16), and (15,32), we see that the edges exhibit a weak but nonzero signature of P e,HORNET , the centers of eddies have nearly exclusively negligible P e,HORNET .
Turning to the ions, Fig. 3(e) shows that almost the whole simulation domain has noticeable P i,HORNET , implying the phase space densities evolve towards or away from LTE in regions away from the intermittent structures.P i,HORNET is found to be non-zero in the ion-scale eddies, e.g., at regions centered around (x, y) ≃ (12, 12), (30,27), and (29,16).This makes sense physically, as ions are prone to get demagnetized in ion-scale magnetic islands since their characteristic ion gyroradius in the mean field is 0.77.
The peak amplitude of P e,HORNET is higher by a factor of ≃ 4 than P i,HORNET at this time in the simulation, implying the peak effective power density of evolution towards or away from LTE is higher for electrons than ions.Since the density and temperature for the ions and electrons are initially equal, Eq. ( 9) suggests the electrons become more non-Maxwellian than the ions.

B. Comparison to other power densities
We next compare P σ,HORNET with power densities describing changes to internal energy.Pressure dilatation and Pi − D σ have been the subject of many recent 2D turbulence simulation studies [9,10,20,21,56], with a greater emphasis on Pi − D σ .Here, we discuss the overall 2D structures in −P σ (∇ • u σ ) and Pi − D σ and compare their structures with −∇ • q σ and P σ,HORNET .
We first analyze electron power densities at t = 30, shown in Fig. 3.We find that −P e (∇ • u e ) (panel b)) has localized regions of compression or expansion on the order of 1 d i (= 5 d e ) in size and up to an amplitude of 0.16, unlike the structures seen in P e,HORNET .This is presumably due to electron plasma waves propagating through the system.Looking at Pi − D e (panel (c)), we only see a few localized regions with prominent amplitudes up to an amplitude of about 0.06 that are a few (∼ 5) d i in length and electron-scale in thickness, which is similar to the structures seen in P e,HORNET .Finally, the 2D structures seen in −∇ • q e (panel (d)) show localized regions with significant amplitude up to 0.160 that are co-located with current sheets, many of which coincide with regions of strong P e,HORNET and Pi − D e .However, −∇ • q e is appreciably non-zero with positive and negative values throughout the simulation domain, with structure sizes on electron scales.Furthermore, at (x, y) ≃ (15, 7), the structure in −∇ • q e is not coherent while it is in the reconnection simulation, even though P e,HORNET at the same location looks similar in the two.Amongst the four electron power densities, −P e (∇ • u e ) and −∇ • q e are the terms with the largest local magnitudes, but there are more regions with the peak values in −P e (∇ • u e ) than −∇ • q e as −∇ • q e is strongest where reconnection is taking place while −P e (∇ • u e ) is strongest due to plasma waves pervasive in the domain.
We turn to ion power densities.The 2D structures seen for −P i (∇ • u i ) (panel (f)) and Pi − D i (panel (g)) are of the order of a few d i size, similar to what is seen in P i,HORNET .Both have their highest magnitude near intermittent structures (≃ 0.055 for pressure dilatation, ≃ 0.009 for Pi − D i ), with smaller magnitudes elsewhere in the domain, and both quantities have various regions with positive and negative values.−∇ • q i (panel (h)) shows structures similar to the other power densities, with an amplitude up to 0.064, with one notable difference that the structures are slightly smaller in size.Analogous to the electrons, the largest local magnitudes in ion power density are in −∇ • q i and −P i (∇ • u i ), with the former having slightly higher peak magnitude.
Next, we compare the signed spatial average of the power densities considered here calculated over the whole simulation domain as a function of time in Fig. 4, with the same color scheme as in Fig. 2. Panel (a) is for electrons, and panel (b) is for ions with the inset in the latter focusing from t = 10 − 50.A similar analysis of ⟨−P σ (∇ • u σ )⟩ and ⟨Pi − D σ ⟩ was performed in Ref. [11] (see their Figs.2(c) and (d)), which utilized a similar simulation as used in the present study, except their simulation employed a lower PPG of 3,200, had a bigger simulation domain of 150 2 d 2 i and was evolved for more than 300 Ω −1 ci .We note that ⟨−∇ • q σ ⟩ ≃ 0 for both species, as it should, since we are calculating the average over the entire computational domain.
We see in panel (b) that, up until about 1 τ nl , there are oscillations in ⟨Pi − D i ⟩ and ⟨P i,HORNET ⟩ that dwarf all other power densities.We attribute these oscillations to Alfvénic exchange between different modes from the initiation of the simulation, and therefore focus on times after t ≳ 15 ≃ τ nl .
For both electrons and ions, ⟨Pi − D σ ⟩ is the dominant contribution to increasing the internal energy for this simulation, and there are oscillations around zero for ⟨−P e (∇ • u e )⟩, consistent with Ref. [11].The value for ⟨Pi − D σ ⟩ is comparable for the two species, starting at ≃ 0.0003 at t ≃ 1 τ nl and slowly decreasing throughout the rest of the evolution as the turbulence decays.We find ⟨P σ,HORNET ⟩ is ≃ 0.0002 at t ≃ τ nl , which is ≃ 67% of ⟨Pi − D σ ⟩ at that time.While acknowledging that the smallness of the averages increases the potential impact of uncertainties, our turbulence simulation suggests the effective power density associated with the phase space density evolving towards or away from LTE is a significant percentage of the overall energy budget at this time.
We observe that ⟨P σ,HORNET ⟩ decays in time faster than ⟨Pi − D σ ⟩.By t = 30 ≃ 2 τ nl , we find that ⟨P e,HORNET ⟩ ≃ 31% ⟨Pi − D e ⟩ (down from 67% at t ≃ 1 τ nl ) and ⟨P i,HORNET ⟩ ≃ 50% ⟨Pi − D i ⟩ (down from 67% at t ≃ 1 τ nl ).Physically, this implies that the plasma thermalizes more rapidly than it cools as the turbulence decays, whether through physical or numerical effects.Interestingly, we find that ⟨P e,HORNET ⟩ peaks at t = 14, earlier than the mean square current density ⟨J 2 ⟩ (dominated by electrons) which peaks at t = 21 (plots not shown).At t = 14, the domain has a single intermittent structure with strongly positive P e,HORNET , with its spatial distribution more closely resembling that seen in reconnection simulation in Fig. 1(a) than the example discussed at t = 30 in Fig. 3(a).Consequently, the peak effective power density associated with electrons evolving away from LTE occurs prior to current sheets reaching their maximum current as eddies collide.We attribute this to the fact that electron phase space densities start to become non-Maxwellian before electrons demagnetize in the current sheets due to electron trapping in mirror magnetic fields present near the current sheets [45].
Finally, we briefly compare P σ,HORNET and J σ • E. in the decaying turbulence simulation at t = 30.As with the pattern seen in −P e (∇ • u e ) in Fig. 3, J e • E exhibits localized regions of positive and negative amplitudes (not shown).We posit that the structure may be attributed to propagating electron plasma waves within the simulation domain.This hypothesis is supported by the proximity of the peak amplitude of J e • E (0.157) to that of −P e (∇ • u e ), which is an order of magnitude higher than P e,HORNET .For the ions, the structures in J i • E bear resemblance to those observed in −P i (∇ • u i ) and Pi − D i , with peak amplitudes approximately two times larger than the other channels of energy conversion illustrated in Fig. 3 (not shown).This higher value of J σ • E occurs because it quantifies the rate of energy conversion from the fields to the particles, while the power densities considered here contain information about changes to internal energy but not bulk kinetic energy.

VI. DISCUSSION AND CONCLUSION
In this work, we define the "higher order nonequilibrium term" (HORNET) effective power density P σ,HORNET [see Eq. ( 9)].It is local in space and time and quantifies the rate at which a species evolves towards or away from local thermodynamic equilibrium.Positive P σ,HORNET implies the phase space density is locally becoming less Maxwellian in time, while negative P σ,HORNET implies it is locally becoming more Maxwellian (thermalizing).The motivation for defining P σ,HORNET is that it is an effective power density, on the same footing as well-studied power densities describing changes to internal energy, namely pressure dilatation −P σ (∇ • u σ ) describing compression, Pi − D σ describing incompressible effects, and the divergence of the vector heat flux density −∇ • q σ [see Eq. ( 3)], as well as J σ • E. This allows the terms to be compared to determine the relative importance of HORNET compared to the power densities.
To exemplify the utility of P σ,HORNET , we use numerical simulations of plasmas that are not in LTE to study its spatial and temporal variation for both electrons and ions.We use fully kinetic PIC simulations of magnetic reconnection at the time when the reconnection rate is most rapidly increasing.We first look at the spatial variation of P σ,HORNET in and around the EDR.We find that P σ,HORNET identifies regions where particle phase space densities locally evolve toward or away from Maxwellianity in the comoving (Lagrangian) reference frame.We also find that P σ,HORNET is locally significant at and near the X-line for both electrons and ions in our simulation.When compared to power densities describing changes to the internal energy, we find that P e,HORNET ≃ 23% of the sum of −P e (∇ • u e ), Pi − D e , and −∇ • q e , and P i,HORNET ≃ 65% of the sum of −P i (∇ • u i ), Pi − D i , and −∇ • q i .We further examine the spatial average of P σ,HORNET inside the EDR and compare it with the other power densities as a function of time.Though smaller in We also study PIC simulations of decaying turbulence.We investigate the 2D spatial distribution of P σ,HORNET and compare it with other power density measures.Quantitatively, we find P σ,HORNET is comparable to other power densities.Similar to the reconnection simulation, we look at the time evolution of the spatial average of these energy conversion metrics, but averaged over the whole simulation domain.We find ⟨P σ,HORNET ⟩ is a substantial fraction of the only other significant power density ⟨Pi − D σ ⟩; at one non-linear time into the evolution of the system it is 67% for both electrons and for ions.
It is tempting to directly compare the prevalence of non-LTE effects in the reconnection and decaying turbulence simulations in this study.However, we carry out spatial averaging for P σ,HORNET in the reconnection simulation over a localized domain, approximately the EDR, but for the whole domain in the decaying turbulence simulation, so the results for the percentage of P σ,HORNET compared to the power densities are not describing the same physical quantity.Moreover, we find that the signed average of P σ,HORNET is lower in the decaying turbulence simulation than in the reconnection simulation.This difference can be attributed to the fact that the regions undergoing reconnection, embedded within the turbulence simulation, are relatively small compared to the overall simulation domain.Furthermore, the turbulence simulation contains regions with both positive and negative P σ,HORNET , and these oppose each other when averaged over the entire domain.Thus, it is essential to exercise caution when comparing quantities across systems.
We now address the implications of the present work.We discuss how one might expect the magnitudes of electron and ion HORNET effective power density to compare to each other.Using Eq. ( 9), their ratio is For a singly ionized two-component quasi-neutral plasma, the densities cancel.Thus, the electron-to-ion HOR-NET ratio is the electron-to-ion temperature ratio times the ratio of the rate of change of the relative entropy per particle.Calculating this ratio provides information about which phase space density is more rapidly changing its shape.For example, averaging over the EDR at t = 13 in the reconnection simulation, we find ⟨P e,HORNET ⟩/⟨P i,HORNET ⟩ ≃ 0.14.Since the initial temperature ratio T e /T i = 0.2, we conclude the two species are non-Maxwellianizing at a comparable rate, with the ions being at a slightly faster rate.Similarly, averaging over the whole domain at t = 30 in the turbulence simulation, ⟨P e,HORNET ⟩/⟨P i,HORNET ⟩ ≃ 1.1.Since the initial temperature ratio is 1, we conclude the two species are non-Maxwellianizing at a comparable rate, with the electrons being slightly faster.
A second important implication applies for collisionless or weakly collisional plasmas, such as some in space.If the collisional terms in Eq. ( 10) are negligible, then P σ,HORNET is equivalent to This form could be useful because it is easier to calculate in simulations, and potentially in spacecraft or labora-tory experiments, because it does not require taking a temporal derivative.We now discuss possible applications of this study.In turbulent plasmas, energy injected at large scales cascades down to smaller scales [57], leading to intermittent coherent structures [58].It has been argued [9,10,59] that all channels of energy conversion are present in these structures, which convert between bulk flow, thermal and electromagnetic energies.The present results provide an approach to quantify the relative importance of non-LTE effects, which we have shown can be a significant contribution.Magnetic reconnection in thin current sheets, both within turbulent plasmas [55,[60][61][62][63][64][65][66] and in other solar, magnetospheric, planetary, astrophysical, and laboratory settings [67], have long been investigated as locations of energy conversion.Understanding the contribution of non-LTE physics near the diffusion region and beyond is important to understand this conversion.Collisionless shocks, such as Earth's bow shock, convert bulk flow energy into magnetic energy and thermal and nonthermal energy at a relatively sharp boundary layer [68].Plasmas in and around the boundary layer can be very strongly non-Maxwellian, so quantifying the contribution is expected to be very useful.Beyond plasma physics, we expect potential applications in sciences where non-LTE systems are routinely studied such as condensed matter physics [69], open quantum systems [70], and molecular dynamics in chemistry and biology [71,72].
We expect P σ,HORNET to be useful in a variety of analysis techniques for non-LTE systems for which f σ can be measured such as numerical simulations, laboratory experiments, and spacecraft observations of space plasma.In particular, phase space densities are readily obtainable in particle-in-cell and Vlasov/Boltzmann simulations, thus calculating P σ,HORNET is possible.Slices of phase space densities can be measured in plasma experiments [73][74][75][76][77], although full three-dimensional measurements are currently beyond the present technology.For spacecraft data, it has been shown [36,78,79] that kinetic entropy can reliably be measured using NASA's Magnetospheric Multiscale (MMS) satellites [80].One direct avenue of future work is to measure P σ,HORNET using MMS data during turbulence, reconnection, and collisionless shocks.
Another avenue for future work is to ascertain the impact of the system size on the 2D structure of P i,HORNET during magnetic reconnection.A limitation of the present reconnection study is the small system size; larger simulations are necessary to determine how P i,HORNET will change in a system in which ions fully couple to large-scale physics [53].In a simulation with a larger domain, we expect ions to recouple to the reconnected field near the downstream edges of the IDR and thermalize, which will appear as a P i,HORNET < 0 signature, analogous to the negative P e,HORNET seen for electrons near the downstream edges of the EDR in the present study.Other avenues for future work, for both reconnection and turbulence, include utilizing a more realistic mass ratio, performing 3D simulations, studying the parametric dependence on P σ,HORNET with ambient plasma parameters, and studying the dependence on out-of-plane mean (guide) magnetic field, upstream asymmetries and/or upstream bulk flow.
In the absence of physical collisions as in the PIC simulations discussed here, energy and entropy are conserved reasonably well but are not conserved perfectly due to numerical dissipation.This implies that the physical energy conversion to internal energy by the discussed energy conversion measures, i.e., pressure dilatation, Pi − D σ , and the divergence of the vector heat flux density, is not equal to the rate of change of internal energy density (in the Lagrangian frame) as suggested by Eq. ( 2).Similarly, the rate of change in relative energy in Eq. ( 12) is affected by numerical dissipation.Entropy metrics can be useful for quantifying numerical dissipation (e.g., [43]).The effects of physical collisions and entropic approaches to quantifying numerical dissipation will be topics of future study.We acknowledge basing our color table for 2D plots on the "bam" color table [81] in an effort to improve color-vision deficiency accessibility [82,83].The colorbar extrema are green [(r g , g g , b g ) = (2, 73, 29) = C g ] and purple [(r p , g p , b p ) = (68, 3, 79) = C p ], where (r, g, b) are the red, green, and blue components of each color, given by the subscript g for green and p for purple.Bam varies linearly from purple to white to green, but we stretch the color table to provide higher contrast for small values.Letting C(k) = (r(k), g(k), b(k)) be the color table with index k from 0 to 255, we use a color table of velocity u e,z , the phase space density weighted by v ′2 is more significant at v z < u e,z than v z > u e,z , so there is a negative heat flux density q e,z at the X-line (panel (f)).As one goes away from the X-line in either outflow direction, the phase space density is rotated by the reconnected magnetic field B y along the y − y 0 = 0 line, as has been discussed previously [12,47,48,84].To the left of the X-line, there is a negative u e,x (panel (a)).The part of the phase space density below the bulk velocity u e,z rotates in the positive v x direction, so this rotation creates an asymmetry in f e that gives a positive q e,x , consistent with panel (d).Similarly, to the right of the X-line, q e,x is negative, also consistent with panel (d).This is the kinetic manifestation of the spatial structure in q e,x .This term is the dominant contribution to the converging vec-tor heat flux −∇ • q e at and near the X-line as seen by the green color near the X-line in Fig. 1(d).
Appendix C: Pe,HORNET in the turbulence simulation In Sec.V A, we display a 2D plot of P e,HORNET in a PIC simulation of decaying turbulence at t = 30 and discuss its structure.To highlight the structure beyond the regions of its peak value, we show the same quantity shown in Fig. 3(a) in Fig. 7, but with the colorbar saturated to have the same range as that for P i,HORNET in Fig. 3(e).The regions with non-zero P e,HORNET outside the intermittent structures are more clearly visible here.

FIG. 1 .
FIG. 1. Particle-in-cell simulation results of power densities for electrons and ions during the non-linear growth phase of reconnection (t = 13).(a) and (e) higher order non-equilibrium terms effective power density Pe,HORNET and Pi,HORNET, (b) and (f) pressure dilatation −Pe(∇ • ue) and −Pi(∇ • ui), (c) and (g) Pi − D e and Pi − D i , and (d) and (h) vector heat flux density divergence −∇ • qe and −∇ • qi.
(b), Pi − D σ is shown in Fig.1(c), and Fig.1(d) exhibits the negative divergence of the vector heat flux density −∇ • q σ .These quantities are repeated in Fig.1(f) -(h) for ions.We first discuss electrons.In a recent study with a sim-ulation with the same initial conditions except with PPG = 25,600 instead of 6,400[14], −P e (∇ • u e ) and Pi − D e were extensively discussed, but the vector heat flux density divergence was not analyzed.Near the X-line, we find −∇•q e is positive (the vector heat flux is converging) with a value near 0.04 (panel (d)), which by itself would increase the electron internal energy.The pressure dilatation is negative with a value near −0.015 at the X-line (panel (b)), which by itself would decrease the electron internal energy.Pi − D e is weakly positive at 0.005 near the X-line (panel (c)).The net change to the internal energy is positive at this time, and the sum of the pressure dilatation, Pi − D e , and negative vector heat flux divergence is approximately 0.03.Concomitantly, panel (a) reveals that P e,HORNET is positive at the X-line with a value near 0.007.Thus, the electron phase space densities are becoming more non-Maxwellian with an effective power density approximately 23% of that going to changing the electron internal energy.In contrast, near the outflow edges of the EDR (|x−x 0 | ≃ 1.7), −P e (∇•u e ) ≃ 0.04 and Pi − D e ≃ 0.02, whereas −∇ • q e ≃ −0.01 with P e,HORNET ≃ 0.002.This indicates that the electron internal energy is increasing, while the electron phase space densities are becoming more Maxwellian at this location.

FIG. 4 .
FIG. 4. Instantaneous spatially-averaged power densities in the turbulence simulation domain as a function of time for (a) electrons and (b) ions.The inset in (b) is a plot from time t = 10 − 50.The same color scheme as Fig. 2 is used.In this simulation, τ nl ≃ 17.

FIG. 5 .
FIG.5.Reduced electron phase space densities fe(vx, vy) from a particle-in-cell magnetic reconnection simulation (using the simulation with PPG = 25,600 from Ref.[14]).(a) Inside the EDR, close to its upstream edge, (b) further inside the EDR, and (c) near the X-line.

P
e,HORNET P i,HORNET = n e T e d dt (s ev,rel /n e ) n i T i d dt (s iv,rel /n i ) .
ACKNOWLEDGMENTS P.A.C. gratefully acknowledges primary support from NSF Grant PHY-1804428 and supplemental support from NSF Grant AGS-1602769, NASA Grant 80NSSC19M0146 and DOE Grant DE-SC0020294.M.R.A. acknowledges support by NASA grant 80NSSC23K0409.This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0027083.Data used for the figures are available publicly at https://doi.org/10.5281/zenodo.8147803.