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

M. Hasan Barbhuiya<EMAIL_ADDRESS>Paul A. Cassak  Subash Adhikari
Department of Physics and Astronomy and the Center for KINETIC Plasma Physics,
West Virginia University, Morgantown, WV 26506, USA  Tulasi N. Parashar
School of Chemical and Physical Sciences, Victoria University of Wellington,
Gate 7, Kelburn Parade, Wellington, 6012, New Zealand  Haoming Liang
Department of Astronomy, University of Maryland College Park, College Park, MD
20742, USA and NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Matthew R. Argall  Space Science Center, Institute for the Study of Earth,
Oceans, and Space
University of New Hampshire, Durham, NH 03824, USA

###### Abstract

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.

††preprint: APS/123-QED

##  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 Chapman and Cowling (1970), Onsager Onsager (1931), Grad Grad
(1949) and its extensions Jou  _et al._ (2010), and the gyrokinetic approach
Howes  _et al._ (2006); Schekochihin  _et al._ (2009) 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.

One non-perturbative approach is to use the non-LTE internal (thermal) energy
density evolution equation Braginskii (1965), which follows directly from the
Boltzmann transport equation. As will be reviewed in the next section, four
physical mechanisms can cause the internal energy density to change in the
comoving (Lagrangian) reference frame. They are a) pressure dilatation, which
describes compressible heating; b) incompressible heating, including
deformation and shear, captured by the term called ${\rm Pi-D}$ Del Sarto  _et
al._ (2016); Yang  _et al._ (2017a, b, 2022); Cassak and Barbhuiya (2022);
Cassak  _et al._ (2022); Barbhuiya and Cassak (2022); c) the divergence of the
vector heat flux density Zhou  _et al._ (2021); Song  _et al._ (2020); Du  _et
al._ (2020); and d) inter-species collisions. In LTE, $\mathrm{Pi-D}$, vector
heat flux density, and the collision operator vanish, so these terms describe
non-LTE energy conversion. They have been extensively used to study energy
conversion in space and astrophysical plasmas, especially during magnetic
reconnection and turbulence Del Sarto  _et al._ (2016); Yang  _et al._
(2017b); Chasapis  _et al._ (2018); Yang  _et al._ (2019); Pezzi  _et al._
(2019); Song  _et al._ (2020); Du  _et al._ (2020); Pezzi  _et al._ (2021);
Bandyopadhyay  _et al._ (2021); Zhou  _et al._ (2021); Yang  _et al._ (2022);
Cassak and Barbhuiya (2022); Cassak  _et al._ (2022); Barbhuiya and Cassak
(2022).

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 “field-particle correlation” approach Klein and Howes (2016). 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 Howes
(2017); Klein  _et al._ (2017); Chen  _et al._ (2019); Klein  _et al._ (2020);
Juno  _et al._ (2021); McCubbin  _et al._ (2022); Montag and Howes (2022);
Montag  _et al._ (2023).

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 six-dimensional (plus time)
field. Notably, recent research Cassak  _et al._ (2023) 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 Liang  _et al._ (2019).

Mechanisms that can change the internal energy (pressure dilatation, ${\rm
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, $\mathbf{J}_{\sigma}\cdot\mathbf{E}$ describes the rate of
conversion of energy density between electric fields and charged particles of
species $\sigma$, where $\mathbf{J}_{\sigma}$ is the current density of
particle species $\sigma$ and $\mathbf{E}$ is the electric field. Previous
studies Yang  _et al._ (2017b); Pezzi  _et al._ (2019) have demonstrated that
$\mathbf{J}_{\sigma}\cdot\mathbf{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. Cassak  _et al._ (2023) to define a quantity
$P_{\sigma,{\rm 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_{\sigma,{\rm HORNET}}$ in particle-in-cell (PIC) simulations of magnetic
reconnection and decaying plasma turbulence. By locally comparing
$P_{\sigma,{\rm 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_{\sigma,{\rm 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. Cassak  _et al._ (2023) and introduces
$P_{\sigma,{\rm HORNET}}$. Sec. III describes the numerical simulation 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.

##  II Theory

###  II.1 Review

The non-relativistic Boltzmann transport equation Boltzmann (1872) for species
$\sigma$ is given by

$\frac{\partial f_{\sigma}}{\partial t}+{\bf v}\cdot\nabla
f_{\sigma}+\frac{{\bf F}_{\sigma}}{m_{\sigma}}\cdot\nabla_{v}f_{\sigma}=C[f],$
(1)

where $f_{\sigma}$ is the phase space density, ${\bf F}_{\sigma}$ is the net
body force, $m_{\sigma}$ is the mass of the constituents, ${\bf r}$, ${\bf
v}$, and $t$ are the position-space, velocity-space, and time coordinates of
phase space, $\nabla$ and $\nabla_{v}$ are the position- and velocity-space
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[f]=C_{\sigma}[f_{\sigma}]+\sum_{\sigma^{\prime}\neq\sigma}C_{\sigma\sigma^{\prime}}[f_{\sigma},f_{\sigma^{\prime}}]$.
Multiplying Eq. (1) by $(1/2)m_{\sigma}v_{\sigma}^{\prime 2}$ and integrating
over all velocity space, where ${\bf v}_{\sigma}^{\prime}={\bf v}-{\bf
u}_{\sigma}$ is the random velocity, ${\bf u}_{\sigma}=(1/n_{\sigma})\int
f_{\sigma}{\bf v}d^{3}v$ is the bulk flow velocity and $n_{\sigma}=\int
f_{\sigma}d^{3}v$ is the number density, one obtains the time evolution
equation of internal energy density $\bar{\mathcal{E}}_{\sigma,{\rm
int}}=\int(1/2)m_{\sigma}v_{\sigma}^{\prime 2}f_{\sigma}d^{3}v$ given by
Braginskii (1965)

$\displaystyle\frac{\partial\bar{\mathcal{E}}_{\sigma,{\rm int}}}{\partial t}$
$\displaystyle+$ $\displaystyle\bm{\nabla}\cdot(\bar{\mathcal{E}}_{\sigma,{\rm
int}}{\bf u}_{\sigma})=$ (2) $\displaystyle-$ $\displaystyle({\bf
P}_{\sigma}\cdot\bm{\nabla})\cdot{\bf u}_{\sigma}-\bm{\nabla}\cdot{\bf
q}_{\sigma}+n_{\sigma}\dot{Q}_{\sigma,{\rm coll,inter}},$

where ${\bf P}_{\sigma}=n_{\sigma}k_{B}{\bf T}_{\sigma}$ is the pressure
tensor, ${\bf T}_{\sigma}$ is the temperature tensor with elements
$T_{\sigma,jk}=m_{\sigma}/(n_{\sigma}k_{B})\int v_{\sigma j}^{\prime}v_{\sigma
k}^{\prime}f_{\sigma}d^{3}v$, $k_{B}$ is Boltzmann’s constant, ${\bf
q}_{\sigma}=\int(1/2)m_{\sigma}v_{\sigma}^{\prime 2}{\bf
v}_{\sigma}^{\prime}f_{\sigma}d^{3}v$ is the vector heat flux density, and
$\dot{Q}_{\sigma,{\rm
coll,inter}}=(1/n_{\sigma})\int(1/2)m_{\sigma}v_{\sigma}^{\prime
2}\sum_{\sigma^{\prime}\neq\sigma}C_{\sigma\sigma^{\prime}}[f_{\sigma},f_{\sigma^{\prime}}]d^{3}v$
is the volumetric heating rate per particle due to inter-species collisions.
The pressure-strain interaction is defined as $-({\bf
P}_{\sigma}\cdot\bm{\nabla})\cdot{\bf u}_{\sigma}$ Del Sarto  _et al._ (2016);
Yang  _et al._ (2022); Cassak and Barbhuiya (2022); Yang  _et al._ (2017b, a).
Note, $\bar{\mathcal{E}}_{\sigma,{\rm int}}=3{\cal
P}_{\sigma}/2=(3/2)n_{\sigma}k_{B}\mathcal{T}_{\sigma}$, where we call ${\cal
P}_{\sigma}={\rm tr}({\bf P}_{\sigma})/3$ the effective pressure and ${\cal
T}_{\sigma}={\rm tr}({\bf T}_{\sigma})/3$ the effective temperature, and tr is
the trace. The lack of the ${\bf J}_{\sigma}\cdot{\bf E}$ term in this
equation is used as evidence that the term does not directly change the
internal energy density.

It is advantageous to write Eq. (2) in the comoving (Lagrangian) reference
frame using the convective derivative $d/dt=\partial/\partial t+{\bf
u}_{\sigma}\cdot\nabla$, which gives Braginskii (1965)

$\frac{d{\cal E}_{\sigma,{\rm int}}}{dt}=\frac{-{\cal
P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})-\Pi_{\sigma,jk}{\cal
D}_{\sigma,jk}-\bm{\nabla}\cdot{\bf
q}_{\sigma}}{n_{\sigma}}+\dot{Q}_{\sigma,{\rm coll,inter}},$ (3)

where ${\cal E}_{\sigma,{\rm int}}=(3/2)k_{B}{\cal T}_{\sigma}$ (without the
overbar) is the internal energy per particle. In writing Eq. (3), we decompose
the pressure-strain interaction into a sum of pressure dilatation $-{\cal
P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})$ and ${\rm
Pi-D}_{\sigma}=-\Pi_{\sigma,jk}{\cal D}_{\sigma,jk}$ that describe
compressible and incompressible energy conversion between bulk kinetic energy
and internal energy, respectively, and we use the Einstein summation
convention for repeated dummy indices. Here,
$\Pi_{\sigma,jk}=P_{\sigma,jk}-{\cal P}_{\sigma}\delta_{jk}$ is the $jk$
element of the deviatoric pressure tensor, ${\cal
D}_{\sigma,jk}=(1/2)\left(\partial u_{\sigma j}/\partial r_{k}+\partial
u_{\sigma k}/\partial r_{j}\right)-(1/3)\delta_{jk}(\bm{\nabla}\cdot{\bf
u}_{\sigma})$ is the $jk$ element of the traceless symmetric strain rate
tensor, and $\delta_{jk}$ is the Kronecker delta Yang  _et al._ (2017b).

In Ref. Cassak  _et al._ (2023), it was shown using an entropic approach that
for systems not in LTE with a fixed total number of particles $N_{\sigma}$ of
each species, the so-called first law of kinetic theory is given by

$\frac{d{\cal W}_{\sigma}}{dt}+\frac{d{\cal E}_{\sigma,{\rm
gen}}}{dt}=\frac{d{\cal Q}_{\sigma,{\rm gen}}}{dt}+\dot{{\cal Q}}_{\sigma,\rm
coll},$ (4)

where $d{\cal E}_{\sigma,{\rm gen}}$ and $d{\cal Q}_{\sigma,{\rm gen}}$ are
increments in the so-called generalized energy per particle and generalized
heat per particle, respectively, and

$\displaystyle\frac{d{\cal W}_{\sigma}}{dt}$ $\displaystyle=$
$\displaystyle{\cal P}_{\sigma}\frac{d(1/n_{\sigma})}{dt},$ (5a)
$\displaystyle\frac{d{\cal E}_{\sigma,{\rm gen}}}{dt}$ $\displaystyle=$
$\displaystyle\frac{d{\cal E}_{\sigma,{\rm int}}}{dt}+\frac{d{\cal
E}_{\sigma,{\rm rel}}}{dt},$ (5b) $\displaystyle\frac{d{\cal Q}_{\sigma,{\rm
gen}}}{dt}$ $\displaystyle=$ $\displaystyle\frac{d{\cal
Q}_{\sigma}}{dt}+\frac{d{\cal Q}_{\sigma,{\rm rel}}}{dt},$ (5c)
$\displaystyle\dot{{\cal Q}}_{\sigma,{\rm coll}}$ $\displaystyle=$
$\displaystyle-\frac{k_{B}\mathcal{T}_{\sigma}}{n_{\sigma}}\int
C[f]\ln\left(\frac{f_{\sigma}\Delta^{3}r_{\sigma}\Delta^{3}v_{\sigma}}{N_{\sigma}}\right)d^{3}v.$
(5d)

Equation (5a) describes the rate of change of compressional work per particle
$d{\cal W}_{\sigma}/dt$ done by the system. The two terms in Eq. (5b) describe
the rate of change of internal energy per particle $d{\cal E}_{\sigma,{\rm
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 $d{\cal E}_{\sigma,{\rm rel}}/dt$, respectively. The two terms in Eq.
(5c) describe the heating rate per particle associated with changes to the
internal energy $d{\cal Q}_{\sigma}/dt=(-\nabla\cdot{\bf
q}_{\sigma}-\Pi_{\sigma,jk}{\cal D}_{\sigma,jk})/n_{\sigma}$ [see Eq. (3)] and
an effective rate of heat per particle associated with the local phase space
density not being in LTE $d{\cal Q}_{\sigma,{\rm rel}}/dt$, respectively.
Equation (5d) describes the rate of energy change per particle associated with
collisions. In Eq. (5d), $\Delta^{3}r_{\sigma}\Delta^{3}v_{\sigma}$ is the
six-dimensional phase space volume element Liang  _et al._ (2019, 2020a);
Argall  _et al._ (2022).

Quantitatively, the two relative terms associated with the departure of the
phase space density from LTE are given by

$\displaystyle\frac{d{\cal E}_{\sigma,{\rm rel}}}{dt}$ $\displaystyle=$
$\displaystyle{\cal T}_{\sigma}\frac{d(s_{\sigma v,{\rm
rel}}/n_{\sigma})}{dt},$ (6a) $\displaystyle\frac{d{\cal Q}_{\sigma,{\rm
rel}}}{dt}$ $\displaystyle=$ $\displaystyle-{\cal
T}_{\sigma}\frac{(\nabla\cdot\bm{\mathcal{J}}_{\sigma,{\rm th}})_{{\rm
rel}}}{n_{\sigma}},$ (6b)

where the relative entropy per particle Grad (1965) and the thermal relative
entropy density flux divergence per particle Cassak  _et al._ (2023) are
defined as

$\displaystyle\frac{s_{\sigma v,{\rm rel}}}{n_{\sigma}}$ $\displaystyle=$
$\displaystyle-\frac{k_{B}}{n_{\sigma}}\int
f_{\sigma}\ln\left(\frac{f_{\sigma}}{f_{\sigma M}}\right)d^{3}v,$ (7a)
$\displaystyle\frac{(\bm{\nabla}\cdot\bm{\mathcal{J}}_{\sigma,{\rm th}})_{{\rm
rel}}}{n_{\sigma}}$ $\displaystyle=$
$\displaystyle-\frac{k_{B}}{n_{\sigma}}\int\left[\frac{}{}\bm{\nabla}\cdot({\bf
v}_{\sigma}^{\prime}f_{\sigma})\right]\ln\left(\frac{f_{\sigma}}{f_{\sigma
M}}\right)d^{3}v.$ (7b)

Here, $f_{\sigma M}$ is the “Maxwellianized” phase space density associated
with $f_{\sigma}$ given by Grad (1965)

$f_{\sigma M}=n_{\sigma}\left(\frac{m_{\sigma}}{2\pi k_{B}{\cal
T}_{\sigma}}\right)^{3/2}e^{{-m_{\sigma}({\bf v}-{\bf
u}_{\sigma})^{2}/2k_{B}{\cal T}_{\sigma}}},$ (8)

with $n_{\sigma}$, ${\bf u}_{\sigma}$, and ${\cal T}_{\sigma}$ obtained from
$f_{\sigma}$.

Because $s_{\sigma v,{\rm rel}}/n_{\sigma}$ identically vanishes if
$f_{\sigma}=f_{\sigma M}$ Grad (1965), 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.

###  II.2 The Higher Order Non-Equilibrium Term Power Density $P_{\sigma,{\rm
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_{\sigma}$. We, therefore, argue the
analogous effective power density for describing the evolution of a phase
space density towards or away from LTE is $-n_{\sigma}d{\cal E}_{\sigma,{\rm
rel}}/dt$. We call this the HORNET effective power density $P_{\sigma,{\rm
HORNET}}$:

$P_{\sigma,{\rm HORNET}}=-n_{\sigma}\frac{d{\cal E}_{\sigma,{\rm
rel}}}{dt}=-n_{\sigma}\mathcal{T}_{\sigma}\frac{d}{dt}\left(\frac{s_{\sigma
v,{\rm rel}}}{n_{\sigma}}\right).$ (9)

Since a Maxwellian phase space density is the maximum entropy state for a
fixed number of particles and internal energy Boltzmann (1872), a positive
$P_{\sigma,{\rm 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_{\sigma,{\rm 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 ${\cal
P}_{\sigma}$.

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
$d{\cal E}_{\sigma,{\rm rel}}$ (see also Ref. Eu (1995)). Substituting Eqs.
(5a) - (5c) into Eq. (4) and using Eq. (3) gives

$\frac{d{\cal E}_{\sigma,{\rm rel}}}{dt}=\frac{d{\cal Q}_{\sigma,{\rm
rel}}}{dt}+\dot{Q}_{\sigma,{\rm coll}}-\dot{Q}_{\sigma,{\rm coll,inter}}.$
(10)

Assuming collisions conserve particle number and are elastic, the collision
operator satisfies the relations $\int C_{\sigma}d^{3}v=\int
C_{\sigma\sigma^{\prime}}d^{3}v=0$ from conservation of particle number, $\int
m_{\sigma}{\bf v}C_{\sigma}d^{3}v=0$ from conservation of momentum for intra-
species collisions, and $\int(1/2)m_{\sigma}v^{2}C_{\sigma}d^{3}v=0$ from
conservation of energy for intra-species collisions. A brief derivation
reveals that $\dot{Q}_{\sigma,{\rm coll,inter}}$ is equivalent to

$\dot{Q}_{\sigma,{\rm
coll,inter}}=-\frac{k_{B}\mathcal{T}_{\sigma}}{n_{\sigma}}\int
C[f]\ln\left(\frac{f_{\sigma
M}\Delta^{3}r_{\sigma}\Delta^{3}v_{\sigma}}{N_{\sigma}}\right)d^{3}v,$ (11)

where $f_{\sigma M}$ is defined in Eq. (8). Then, Eq. (10) becomes

$\frac{d{\cal E}_{\sigma,{\rm rel}}}{dt}=\frac{d{\cal Q}_{\sigma,{\rm
rel}}}{dt}+\dot{Q}_{\sigma,{\rm coll,rel}},$ (12)

where

$\dot{Q}_{\sigma,{\rm
coll,rel}}=-\frac{k_{B}\mathcal{T}_{\sigma}}{n_{\sigma}}\int
C[f]\ln\left(\frac{f_{\sigma}}{f_{\sigma M}}\right)d^{3}v.$ (13)

Eq. (12) is a special case of Eq. (3.13) from Ref. Eu (1995). 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 Cassak
_et al._ (2023) because the thermodynamic heat $d{\cal Q}_{\sigma}/dt$ terms
depend on non-Maxwellian features of the phase space density, namely
${\bm{\Pi}}_{\sigma}$ and ${\bf q}_{\sigma}$, and their evolution is included
in $d{\cal E}_{\sigma,{\rm 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_{\sigma v,{\rm rel}}/n_{\sigma}$
from Eq. (7a) as

$\displaystyle\frac{s_{\sigma v,{\rm rel}}}{n_{\sigma}}$ $\displaystyle=$
$\displaystyle-k_{B}\int
f^{\prime}_{\sigma}\ln\left(\frac{f^{\prime}_{\sigma}}{f^{\prime}_{\sigma
M}}\right)d^{3}v,$ (14)

where we define the distribution function as
$f_{\sigma}^{\prime}=f_{\sigma}/n_{\sigma}$. Here, $f_{\sigma}^{\prime}$ is
normalized to 1, while $f_{\sigma}$ is normalized to $n_{\sigma}$. By
definition, $f_{\sigma}^{\prime}$ and $f_{\sigma M}^{\prime}$ are independent
of $n_{\sigma}$, i.e., they are both unchanged if $f_{\sigma}$ changes in such
a way that its shape does not change even if the density does change. Thus
$s_{\sigma v,{\rm rel}}/n_{\sigma}$ is independent of changes solely in
$n_{\sigma}$. Since the HORNET effective power density is proportional to
$d(s_{\sigma v,{\rm rel}}/n_{\sigma})/dt$, it is zero for any process that
only changes the density.

##  III Numerical Simulations

Numerical simulations are performed using p3d Zeiler  _et al._ (2002), 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 Birdsall and Langdon (1991) and
electromagnetic fields are stepped using the trapezoidal leapfrog method
Guzdar  _et al._ (1993). To enforce Poisson’s equation, p3d utilizes the
multigrid method Trottenberg  _et al._ (2000) 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/\omega_{pi0}$, where $c$ is the speed of
light, $\omega_{pi0}=(4\pi n_{0}q_{i}^{2}/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\pi m_{i}n_{0})^{1/2}$. Times are normalized to the
inverse ion cyclotron frequency $\Omega_{ci0}^{-1}=(q_{i}B_{0}/m_{i}c)^{-1}$.
Temperatures are normalized to $m_{i}c_{A0}^{2}/k_{B}$, vector heat flux
densities are in units of $(B_{0}^{2}/4\pi)c_{A0}$, and power densities are in
units of $(B_{0}^{2}/4\pi)\Omega_{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. Argall  _et al._
(2022) and optimize the velocity-space grid scale by checking the agreement
between the simulated kinetic entropy density $s_{\sigma}$ for various $\Delta
v_{\sigma}$ (where $\sigma=e$ or $i$ for electrons or ions) and the
theoretical value at $t=0$ Liang  _et al._ (2020b).

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.

###  III.1 Magnetic Reconnection Simulations

The reconnection simulation system size is $L_{x}\times L_{y}=12.8\times 6.4$
with $1024\times 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{\left[(y-L_{y}/4)/w_{0}\right]}-\tanh{\left[(y-3L_{y}/4)/w_{0}\right]}-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})]\\}\\{{\rm
sech}^{2}\left[(y-L_{y}/4)/w_{0}\right]+{\rm
sech}^{2}\left[(y-3L_{y}/4)/w_{0}\right]\\}+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 X- and
O-line pair in each of the two current sheets to initiate reconnection using a
magnetic field perturbation of the form $\delta B_{x}=-B_{pert}\sin\left(2\pi
x/L_{x}\right)\sin\left(4\pi y/L_{y}\right)$ and $\delta
B_{y}=B_{pert}\left[L_{y}/(2L_{x})\right]\cos\left(2\pi
x/L_{x}\right)\left[1-\cos\left(4\pi y/L_{y}\right)\right]$ 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 $\lambda_{De}=0.0176$; we
choose a grid-length $\Delta=0.0125$. The smallest time scale of the system is
the inverse of electron plasma frequency $\omega^{-1}_{pe}=0.012$ and we
choose the particle time-step $\Delta 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,\sigma},v_{lim,\sigma}]$ Liang  _et al._
(2019), where $v_{lim,\sigma}$ 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,\sigma}$ bins, where
$n_{bin,e}=26$, and $n_{bin,i}=34$. This gives the optimal velocity space grid
scales $\Delta v_{\sigma}$, which are $\Delta v_{e}=1.092$ and $\Delta
v_{i}=0.459$ using the procedure described in Refs. Liang  _et al._ (2019,
2020a, 2020b); Argall  _et al._ (2022). The agreement with the theoretical
values at $t=0$ are within $\pm 0.7\%$ and $\pm 1.5\%$ in the far upstream
region for electrons and ions, respectively, and within $\pm 2\%$ and $\pm
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 Cassak  _et al._ (2023). We use a spatial domain of size
$0.0625\times 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}$.

###  III.2 Turbulence Simulations

The turbulence simulation uses a periodic domain of size $L_{x}\times
L_{y}=37.3912\times 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 $\hat{z}$. The velocity and magnetic field
fluctuations are Alfvènic with random phase and are excited in a band of wave
numbers $k\in[2,4]\times 2\pi/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. Parashar  _et al._ (2018).

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 $\lambda_{De}\simeq 0.03651$ and we choose a
grid-length $\Delta=\lambda_{De}$. The smallest time scale of the system is
the inverse of electron plasma frequency $\omega^{-1}_{pe}\simeq 0.0133$; we
choose the particle time-step to be $\Delta t=0.005$ with the field time step
$\Delta 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
$\tau_{nl}=L_{x}/[2\pi(\delta b_{rms}^{2}+\delta v_{rms}^{2})^{1/2}]\simeq
17$.

For calculating entropies, we obtain an optimal electron velocity space grid
scale of $\Delta v_{e}=2.028$ using $v_{lim,e}=14.2$, and $n_{bin,e}=14$. At
$t=0$, the choice for $\Delta v_{e}$ gives agreement with the theoretical
values within $\pm 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 $\Delta v_{i}=0.228$ which gives agreement
within $\pm 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.

###  III.3 Relative Entropy Implementation

In a previous study Cassak  _et al._ (2023), the calculation of
$(d/dt)(s_{\sigma v,{\rm rel}}/n_{\sigma})$ was carried out by post-processing
output data of $s_{\sigma}$ 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_{\sigma v,{\rm
rel}}/n_{\sigma})$ are similar when calculated in the two different ways (not
shown).

To assess the importance of numerical error, we note that $f_{\sigma}$ at
$t=0$ in both the reconnection and turbulence simulations are loaded as
drifting Maxwellians. Theoretically, from Eq. (7a), we would expect $s_{\sigma
v,{\rm rel}}/n_{\sigma}$ 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_{\sigma v,{\rm rel}}/n_{\sigma}$ is non-
zero. We confirm that the error of non-zero $s_{\sigma v,{\rm
rel}}/n_{\sigma}$ 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,{\rm rel}}/n_{e}$ over the whole domain at $t=0$ is
$\simeq-0.013$. In a test simulation with PPG = 400, we find the average value
of $s_{ev,{\rm rel}}/n_{e}$ at $t=0$ is $\simeq-0.08$. Noise in PIC
simulations scales with $1/\sqrt{\rm 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/\sqrt{{\rm PPG}}$
until the signal is higher than the noise.

##  IV Results: Magnetic Reconnection

###  IV.1 Spatial distribution and amplitude of HORNET

Figure 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 $P_{e,{\rm
HORNET}}$ and $P_{i,{\rm HORNET}}$, (b) and (f) pressure dilatation $-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ and $-{\cal P}_{i}(\bm{\nabla}\cdot{\bf
u}_{i})$, (c) and (g) ${\rm Pi-D}_{e}$ and ${\rm Pi-D}_{i}$, and (d) and (h)
vector heat flux density divergence $-\bm{\nabla}\cdot{\bf q}_{e}$ and
$-\bm{\nabla}\cdot{\bf q}_{i}$.

Simulation results of $P_{\sigma,{\rm 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 X-line location at $(x_{0},y_{0})=(9.59375,1.59375)$. The EDR is
approximately located at $|x-x_{0}|\lesssim 2,|y-y_{0}|\lesssim 0.4$.

We begin with a discussion of the spatial structure of $P_{\sigma,{\rm
HORNET}}$. For the electrons, we note that Ref. Cassak  _et al._ (2023)
included a plot in their Fig. 3(e) of $d{\cal E}_{e,{\rm rel}}/dt$ for a
reconnection simulation that only differed from the present simulation by its
PPG. Since $P_{\sigma,{\rm HORNET}}=-n_{\sigma}d{\cal E}_{e,{\rm rel}}/dt$
from Eq. (9), the structure of $P_{e,{\rm HORNET}}$ is largely similar to
$-d{\cal E}_{e,{\rm rel}}/dt$ from Ref. Cassak  _et al._ (2023). 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 Egedal  _et al._ (2013).
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,{\rm HORNET}}\simeq 0.003$. Near the X-line, triangular striated phase
space densities develop due to electron Speiser orbits Speiser (1965); Ng  _et
al._ (2011); Shuster  _et al._ (2015), a further evolution away from
Maxwellianity, therefore associated with positive $P_{e,{\rm HORNET}}\simeq
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,{\rm
HORNET}}\simeq-0.0015$ is weakly negative. Physically, this happens because
the upstream phase space densities are elongated in the parallel direction
Egedal  _et al._ (2013), 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,{\rm HORNET}}$. Phase
space densities in relevant regions are plotted in Appendix A. Thermalization
of the outflowing electron jets downstream of the X-line Egedal  _et al._
(2012); Lavraud  _et al._ (2016); Wang  _et al._ (2016); Norgren  _et al._
(2020) is observed at the EDR outflow edge ($1.6<|x-x_{0}|<2.0$), associated
with negative $P_{e,{\rm HORNET}}\simeq-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 $\geq 40\times 2.24\simeq 90$ for the ions to fully couple Sharma
Pyakurel  _et al._ (2019). Upstream of the EDR, Fig. 1(e) demonstrates that
$P_{i,{\rm HORNET}}$ is weakly positive and gradually increases approaching
the X-line, where it has a maximum value of $\simeq 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,{\rm
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,{\rm 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}|\simeq 1.6$ along the separatrices outside of the EDR, $P_{i,{\rm
HORNET}}$ is weakly negative with a value of $\simeq-0.001$. This implies the
ion phase space density becomes more Maxwellian (in the Lagrangian frame)
there.

###  IV.2 Comparison to other power densities

We first perform a comparison of $P_{\sigma,{\rm 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 $-{\cal P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})$ is
shown in Fig. 1(b), ${\rm Pi-D}_{\sigma}$ is shown in Fig. 1(c), and Fig. 1(d)
exhibits the negative divergence of the vector heat flux density
$-\bm{\nabla}\cdot{\bf q}_{\sigma}$. These quantities are repeated in Fig.
1(f) - (h) for ions.

We first discuss electrons. In a recent study with a simulation with the same
initial conditions except with PPG = 25,600 instead of 6,400 Barbhuiya and
Cassak (2022), $-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ and ${\rm
Pi-D}_{e}$ were extensively discussed, but the vector heat flux density
divergence was not analyzed. Near the X-line, we find $-\bm{\nabla}\cdot{\bf
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.
${\rm 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, ${\rm Pi-D}_{e}$, and negative vector heat flux
divergence is approximately 0.03. Concomitantly, panel (a) reveals that
$P_{e,{\rm 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}|\simeq 1.7$), $-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})\simeq
0.04$ and ${\rm Pi-D}_{e}\simeq 0.02$, whereas $-\bm{\nabla}\cdot{\bf
q}_{e}\simeq-0.01$ with $P_{e,{\rm HORNET}}\simeq 0.002$. This indicates that
the electron internal energy is increasing, while the electron phase space
densities are becoming more Maxwellian at this location.

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, ${\rm 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 $-\bm{\nabla}\cdot{\bf q}_{i}$ is positive (the vector heat
flux density is converging) Eastwood  _et al._ (2020) with a value $\simeq
0.06$ (panel (h)), which by itself would increase the ion internal energy.
Panel (e) shows that $P_{i,{\rm 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,
${\rm 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 ${\rm Pi-D}_{i}$ is negligible. We also see that
$-\bm{\nabla}\cdot{\bf q}_{i}\simeq-0.06$, with $P_{i,{\rm HORNET}}\simeq
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 $-{\cal P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})$ and ${\rm Pi-D}_{i}$
were plotted in a previous simulation study Pezzi  _et al._ (2021) 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,
$-{\cal P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})$ is notably more localized around
the X-line in the $\hat{x}$ direction [see Fig. 1(f)], presumably because our
system size is smaller (theirs is 25 $\times 25$); and (ii) ${\rm 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. Pezzi  _et al._ (2021), presumably
because theirs is in the steady state.

We also note the spatial distribution of $-\bm{\nabla}\cdot{\bf q}_{\sigma}$
for electrons and ions reveal structural similarities, except the ion
structures are on larger scales. This is consistent with previous work that
shows $-\bm{\nabla}\cdot{\bf q}_{\sigma}$ can be non-zero and dynamically
significant locally Du  _et al._ (2020); Song  _et al._ (2020), even though it
vanishes when integrated over the whole periodic computational domain Yang
_et al._ (2022). 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 $\langle\rangle$) 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 $\langle P_{\sigma,{\rm
HORNET}}\rangle$, dotted orange denotes $\langle-{\cal
P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})\rangle$, dashed purple denotes
$\langle{\rm Pi-D}_{\sigma}\rangle$, and $\langle-\bm{\nabla}\cdot{\bf
q}_{\sigma}\rangle$ is represented by dash-dotted magenta.

Figure 2: Instantaneous spatially-averaged power densities in the reconnection
electron diffusion region as a function of time for (a) electrons and (b)
ions. $\langle P_{\sigma,{\rm HORNET}}\rangle$ are the solid green lines,
$\langle-{\cal P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})\rangle$ are the
dotted orange lines, $\langle{\rm Pi-D}_{\sigma}\rangle$ are the dashed purple
lines, and $\langle-\bm{\nabla}\cdot{\bf q}_{\sigma}\rangle$ are dash-dotted
magenta lines.

At early time in both panels ($t\lesssim 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\simeq 12$ when the reconnection rate begins
to increase rapidly in time. Consequently, we focus mostly on these later
times.

For the electrons, panel (a) reveals that $\langle{\rm Pi-D}_{e}\rangle$ and
$\langle-{\cal P}_{i}(\bm{\nabla}\cdot{\bf u}_{e})\rangle$ are the two largest
terms, nearly equal early in the non-linear growth phase with $\langle{\rm
Pi-D}_{e}\rangle$ becoming larger in the late non-linear growth phase.
$\langle P_{e,{\rm HORNET}}\rangle$ remains relatively small throughout the
evolution. We know from Fig. 1(a) that $P_{e,{\rm 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 $\langle{\rm Pi-D}_{e}\rangle\simeq 0.008$,
$\langle-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})\rangle\simeq 0.006$ and
$\langle-\bm{\nabla}\cdot{\bf q}_{e}\rangle\simeq-0.002$, giving a net power
density to change the electron internal energy of 0.012. Since $\langle
P_{e,{\rm HORNET}}\rangle\simeq 0.001$, we find $\langle P_{e,{\rm
HORNET}}\rangle\simeq 0.13\ \langle{\rm Pi-D}_{e}\rangle\simeq 0.17\
\langle-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})\rangle$. Because of the
smallness of the average values and potential uncertainties in the
measurements of $P_{\sigma,{\rm HORNET}}$, these fractions could have
significant uncertainties. The raw values indicate that $\langle P_{e,{\rm
HORNET}}\rangle$ 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 $\langle P_{e,{\rm HORNET}}\rangle$ 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,{\rm
HORNET}}$, there is a physical reason to believe this could be valid. As seen
in Fig. 1(a), the strongest negative $P_{e,{\rm 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,{\rm HORNET}}<0$ grow with time (not
shown), leading to the overall change in sign of $\langle P_{e,{\rm
HORNET}}\rangle$.

For the ions shown in panel (b) the dominant positive energy conversion
channel is $\langle-\bm{\nabla}\cdot{\bf q}_{i}\rangle$, which is nearly
balanced by a negative $\langle-{\cal P}_{i}(\bm{\nabla}\cdot{\bf
u}_{i})\rangle$. There is a small net positive $\langle{\rm Pi-D}_{i}\rangle$
late in the non-linear growth phase ($t\geq 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., $\langle
P_{i,{\rm HORNET}}\rangle$, though positive, is small compared to the other
ion terms. This is expected from the 2D plot of $P_{i,{\rm 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 $\langle-\bm{\nabla}\cdot{\bf
q}_{i}\rangle\simeq 0.015$, $\langle{\rm Pi-D}_{i}\rangle\simeq 0.01$, and
$\langle-{\cal P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})\rangle\simeq-0.005$, for a
total of 0.02 going into ion internal energy. In comparison, $\langle
P_{i,{\rm HORNET}}\rangle\simeq 0.007\simeq 70\%\ \langle{\rm
Pi-D}_{i}\rangle\simeq 50\%\ \langle-\bm{\nabla}\cdot{\bf q}_{i}\rangle$.
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_{\sigma,{\rm HORNET}}$ and
$\mathbf{J}_{\sigma}\cdot\mathbf{E}$. At $t=13$, we observe that the peak
amplitude of $\mathbf{J}_{\sigma}\cdot\mathbf{E}$ exceeds that of the largest
energy conversion metric (i.e., $-\bm{\nabla}\cdot{\bf q}_{\sigma}$) by more
than a factor of two (not shown). Moreover, it surpasses $P_{\sigma,{\rm
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_{\sigma,{\rm HORNET}}$ are smaller. In summary, the
results of this section suggest $P_{\sigma,{\rm 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

###  V.1 Spatial distribution and amplitude of HORNET

We now discuss the structure of $P_{\sigma,{\rm HORNET}}$ in the decaying
turbulence simulation using 2D position-space plots of the various power
densities in Fig. 3 at $t=30\sim 2\ \tau_{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 $\langle
J^{2}\rangle$ 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.

Figure 3: Particle-in-cell simulation results comparing power densities for
electrons and ions during decaying turbulence at $t=30$. (a) and (e) higher
order non-equilibrium terms effective power density $P_{e,{\rm HORNET}}$ and
$P_{i,{\rm HORNET}}$, (b) and (f) pressure dilatation $-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ and $-{\cal P}_{i}(\bm{\nabla}\cdot{\bf
u}_{i})$, (c) and (g) ${\rm Pi-D}_{e}$ and ${\rm Pi-D}_{i}$, (d) and (h)
vector heat flux density divergence $-\bm{\nabla}\cdot{\bf q}_{e}$ and
$-\bm{\nabla}\cdot{\bf q}_{i}$.

We first consider electrons. In Fig. 3(a), we see that $P_{e,{\rm HORNET}}$
has substantially non-zero signatures co-located with intermittent structures
(current sheets). The regions with strongest $P_{e,{\rm HORNET}}$ up to a
magnitude of about 0.016 are centered around
$(x,y)\simeq(18,36),~{}(21,2),~{}(29,5),~{}(15,18)$, $(15,4)$, and $(30,22)$.
$P_{e,{\rm HORNET}}$ can have either sign in and around these structures. We
observe one location at $(x,y)\simeq(15,7)$ around a current sheet where the
structure of $P_{e,{\rm 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
Servidio  _et al._ (2009) 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,{\rm HORNET}}$. Interestingly, we observe numerous
bands with non-zero $P_{e,{\rm HORNET}}$ coinciding with extended current
sheets that are multiple ion-scales in length ($\sim 10\ d_{i}$) and electron-
scale in thickness.

Panel (a) makes it appear that $P_{e,{\rm 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 colorbar saturated at the
maximum value used for $P_{i,{\rm HORNET}}$ to better observe the structures
in $P_{e,{\rm HORNET}}$ away from the current sheets. This plot makes it clear
that there are regions with non-zero of $P_{e,{\rm HORNET}}$ away from the
current sheets. There are noticeable features in $P_{e,{\rm HORNET}}$ within
ion-scale eddies of size 5-10 $d_{i}$. For example, in eddies centered around
$(x,y)\simeq(12,12),~{}(29,16)$, and $(15,32)$, we see that the edges exhibit
a weak but non-zero signature of $P_{e,{\rm HORNET}}$, while the centers of
eddies have nearly exclusively negligible $P_{e,{\rm HORNET}}$.

Turning to the ions, Fig. 3(e) shows that almost the whole simulation domain
has noticeable $P_{i,{\rm HORNET}}$, implying the phase space densities evolve
towards or away from LTE in regions away from the intermittent structures.
$P_{i,{\rm HORNET}}$ is found to be non-zero in the ion-scale eddies, e.g., at
regions centered around $(x,y)\simeq(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,{\rm HORNET}}$ is higher by a factor of $\simeq$ 4
than $P_{i,{\rm 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.

###  V.2 Comparison to other power densities

We next compare $P_{\sigma,{\rm HORNET}}$ with power densities describing
changes to internal energy. Pressure dilatation and ${\rm Pi-D}_{\sigma}$ have
been the subject of many recent 2D turbulence simulation studies Yang  _et
al._ (2017a, b); Pezzi  _et al._ (2019, 2021); Agudelo Rueda  _et al._ (2022),
with a greater emphasis on ${\rm Pi-D}_{\sigma}$. Here, we discuss the overall
2D structures in $-{\cal P}_{\sigma}(\bm{\nabla}\cdot{\bf u}_{\sigma})$ and
${\rm Pi-D}_{\sigma}$ and compare their structures with $-\bm{\nabla}\cdot{\bf
q}_{\sigma}$ and $P_{\sigma,{\rm HORNET}}$.

We first analyze electron power densities at $t=30$, shown in Fig. 3. We find
that $-{\cal P}_{e}(\bm{\nabla}\cdot{\bf 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,{\rm HORNET}}$. This is presumably due to electron plasma waves
propagating through the system. Looking at ${\rm 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 ($\sim 5$) $d_{i}$ in length and electron-scale
in thickness, which is similar to the structures seen in $P_{e,{\rm HORNET}}$.
Finally, the 2D structures seen in $-\bm{\nabla}\cdot{\bf 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,{\rm HORNET}}$ and ${\rm Pi-D}_{e}$. However, $-\bm{\nabla}\cdot{\bf
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)\simeq(15,7)$, the structure in $-\bm{\nabla}\cdot{\bf q}_{e}$ is not
coherent while it is in the reconnection simulation, even though $P_{e,{\rm
HORNET}}$ at the same location looks similar in the two. Amongst the four
electron power densities, $-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ and
$-\bm{\nabla}\cdot{\bf q}_{e}$ are the terms with the largest local
magnitudes, but there are more regions with the peak values in $-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ than $-\bm{\nabla}\cdot{\bf q}_{e}$ as
$-\bm{\nabla}\cdot{\bf q}_{e}$ is strongest where reconnection is taking place
while $-{\cal P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$ is strongest due to plasma
waves pervasive in the domain.

We turn to ion power densities. The 2D structures seen for $-{\cal
P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})$ (panel (f)) and ${\rm Pi-D}_{i}$ (panel
(g)) are of the order of a few $d_{i}$ size, similar to what is seen in
$P_{i,{\rm HORNET}}$. Both have their highest magnitude near intermittent
structures ($\simeq 0.055$ for pressure dilatation, $\simeq 0.009$ for ${\rm
Pi-D}_{i}$), with smaller magnitudes elsewhere in the domain, and both
quantities have various regions with positive and negative values.
$-\bm{\nabla}\cdot{\bf 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
$-\bm{\nabla}\cdot{\bf q}_{i}$ and $-{\cal P}_{i}(\bm{\nabla}\cdot{\bf
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 $\langle-{\cal P}_{\sigma}(\bm{\nabla}\cdot{\bf
u}_{\sigma})\rangle$ and $\langle{\rm Pi-D}_{\sigma}\rangle$ was performed in
Ref. Yang  _et al._ (2022) (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_{i}^{2}$ and was evolved for more than $300~{}\Omega^{-1}_{ci}$.
We note that $\langle-\bm{\nabla}\cdot{\bf q}_{\sigma}\rangle\simeq 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 $\tau_{nl}$, there are oscillations
in $\langle{\rm Pi-D}_{i}\rangle$ and $\langle P_{i,{\rm HORNET}}\rangle$ 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\gtrsim 15\simeq\tau_{nl}$.

Figure 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, $\tau_{nl}\simeq 17$.

For both electrons and ions, $\langle{\rm Pi-D}_{\sigma}\rangle$ is the
dominant contribution to increasing the internal energy for this simulation,
and there are oscillations around zero for $\langle-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})\rangle$, consistent with Ref. Yang  _et
al._ (2022). The value for $\langle{\rm Pi-D}_{\sigma}\rangle$ is comparable
for the two species, starting at $\simeq 0.0003$ at $t\simeq 1\ \tau_{nl}$ and
slowly decreasing throughout the rest of the evolution as the turbulence
decays. We find $\langle P_{\sigma,{\rm HORNET}}\rangle$ is $\simeq 0.0002$ at
$t\simeq\tau_{nl}$, which is $\simeq 67\%$ of $\langle{\rm
Pi-D}_{\sigma}\rangle$ 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 $\langle P_{\sigma,{\rm HORNET}}\rangle$ decays in time faster
than $\langle{\rm Pi-D}_{\sigma}\rangle$. By $t=30\simeq 2\ \tau_{nl}$, we
find that $\langle P_{e,{\rm HORNET}}\rangle\simeq 31\%$ $\langle{\rm
Pi-D}_{e}\rangle$ (down from 67% at $t\simeq 1\ \tau_{nl}$) and $\langle
P_{i,{\rm HORNET}}\rangle\simeq 50\%$ $\langle{\rm Pi-D}_{i}\rangle$ (down
from 67% at $t\simeq 1\ \tau_{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 $\langle P_{e,{\rm HORNET}}\rangle$ peaks at
$t=14$, earlier than the mean square current density $\langle J^{2}\rangle$
(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,{\rm 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 Egedal  _et al._
(2013).

Figure 5: Reduced electron phase space densities $f_{e}(v_{x},v_{y})$ from a
particle-in-cell magnetic reconnection simulation (using the simulation with
PPG = 25,600 from Ref. Barbhuiya and Cassak (2022)). (a) Inside the EDR, close
to its upstream edge, (b) further inside the EDR, and (c) near the X-line.

Finally, we briefly compare $P_{\sigma,{\rm HORNET}}$ and
$\mathbf{J}_{\sigma}\cdot\mathbf{E}$. in the decaying turbulence simulation at
$t=30$. As with the pattern seen in
$-\mathcal{P}_{e}(\bm{\nabla}\cdot\mathbf{u}_{e})$ in Fig. 3, ${\bf
J}_{e}\cdot{\bf 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 ${\bf
J}_{e}\cdot{\bf E}$ (0.157) to that of
$-\mathcal{P}_{e}(\bm{\nabla}\cdot\mathbf{u}_{e})$, which is an order of
magnitude higher than $P_{e,{\rm HORNET}}$. For the ions, the structures in
${\bf J}_{i}\cdot{\bf E}$ bear resemblance to those observed in
$-\mathcal{P}_{i}(\bm{\nabla}\cdot\mathbf{u}_{i})$ and ${\rm 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
${\bf J}_{\sigma}\cdot{\bf 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 non-equilibrium term” (HORNET)
effective power density $P_{\sigma,{\rm 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_{\sigma,{\rm
HORNET}}$ implies the phase space density is locally becoming less Maxwellian
in time, while negative $P_{\sigma,{\rm HORNET}}$ implies it is locally
becoming more Maxwellian (thermalizing). The motivation for defining
$P_{\sigma,{\rm 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 $-{\cal P}_{\sigma}(\bm{\nabla}\cdot{\bf
u}_{\sigma})$ describing compression, ${\rm Pi-D}_{\sigma}$ describing
incompressible effects, and the divergence of the vector heat flux density
$-\bm{\nabla}\cdot{\bf q}_{\sigma}$ [see Eq. (3)], as well as ${\bf
J}_{\sigma}\cdot{\bf 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_{\sigma,{\rm 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_{\sigma,{\rm
HORNET}}$ in and around the EDR. We find that $P_{\sigma,{\rm 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_{\sigma,{\rm 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,{\rm HORNET}}\simeq 23\%$ of the sum of $-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})$, ${\rm Pi-D}_{e}$, and
$-\bm{\nabla}\cdot{\bf q}_{e}$, and $P_{i,{\rm HORNET}}\simeq 65\%$ of the sum
of $-{\cal P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})$, ${\rm Pi-D}_{i}$, and
$-\bm{\nabla}\cdot{\bf q}_{i}$. We further examine the spatial average of
$P_{\sigma,{\rm HORNET}}$ inside the EDR and compare it with the other power
densities as a function of time. Though smaller in magnitude than the other
power densities, $\langle P_{\sigma,{\rm HORNET}}\rangle$ represents the
dynamical importance of the electrons departing from LTE.  During the non-
linear growth phase of reconnection (at $t=13$), $\langle P_{e,{\rm
HORNET}}\rangle\simeq 8.3\%$ of the sum of $\langle-{\cal
P}_{e}(\bm{\nabla}\cdot{\bf u}_{e})\rangle$, $\langle{\rm Pi-D}_{e}\rangle$
and $\langle-\bm{\nabla}\cdot{\bf q}_{e}\rangle$ and $\langle P_{i,{\rm
HORNET}}\rangle\simeq 35\%$ of the sum of $\langle-{\cal
P}_{i}(\bm{\nabla}\cdot{\bf u}_{i})\rangle$, $\langle{\rm Pi-D}_{i}\rangle$
and $\langle-\bm{\nabla}\cdot{\bf q}_{i}\rangle$.

We also study PIC simulations of decaying turbulence. We investigate the 2D
spatial distribution of $P_{\sigma,{\rm HORNET}}$ and compare it with other
power density measures. Quantitatively, we find $P_{\sigma,{\rm 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
$\langle P_{\sigma,{\rm HORNET}}\rangle$ is a substantial fraction of the only
other significant power density $\langle{\rm Pi-D}_{\sigma}\rangle$; 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_{\sigma,{\rm 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_{\sigma,{\rm HORNET}}$ compared to the power densities are
not describing the same physical quantity. Moreover, we find that the signed
average of $P_{\sigma,{\rm 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_{\sigma,{\rm 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

$\frac{P_{e,{\rm HORNET}}}{P_{i,{\rm
HORNET}}}=\frac{n_{e}\mathcal{T}_{e}\frac{d}{dt}\left(s_{ev,{\rm
rel}}/n_{e}\right)}{n_{i}\mathcal{T}_{i}\frac{d}{dt}\left(s_{iv,{\rm
rel}}/n_{i}\right)}.$ (15)

For a singly ionized two-component quasi-neutral plasma, the densities cancel.
Thus, the electron-to-ion HORNET 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 $\langle P_{e,{\rm
HORNET}}\rangle/\langle P_{i,{\rm HORNET}}\rangle\simeq 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, $\langle P_{e,{\rm HORNET}}\rangle/\langle P_{i,{\rm
HORNET}}\rangle\simeq 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_{\sigma,{\rm HORNET}}$ is equivalent to

$P_{\sigma,{\rm HORNET}}={\cal T}_{\sigma}(\bm{\nabla}\cdot{\bf
J}_{\sigma,{\rm th}})_{{\rm rel}}.$ (16)

This form could be useful because it is easier to calculate in simulations,
and potentially in spacecraft or laboratory experiments, because it does not
require taking a temporal derivative.

Figure 6: Particle-in-cell magnetic reconnection simulation results (using the
simulation with PPG = 25,600 from Ref. Barbhuiya and Cassak (2022)) to explain
the vector heat flux density profiles observed. Reduced electron phase space
densities $f_{e}(v_{x},v_{z})$ at the (a) green, (b) magenta, and (c) cyan
“x”’s overlayed on the $x$-component of electron heat flux density
$\textbf{q}_{e}$ shown in panel (d). Panels (e) and (f) show the $y$ and
$z$-components of $\textbf{q}_{e}$. The magenta dashed lines in (a) - (c)
denote the bulk flow velocity components $u_{e,x}$ and $u_{e,z}$.

We now discuss possible applications of this study. In turbulent plasmas,
energy injected at large scales cascades down to smaller scales Matthaeus
(2021), leading to intermittent coherent structures Wan  _et al._ (2012). It
has been argued Yang  _et al._ (2017b, a); Matthaeus  _et al._ (2020) 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 Matthaeus and Lamkin
(1986); Servidio  _et al._ (2009, 2010); Haggerty  _et al._ (2017); Retinò
_et al._ (2007); Sundkvist  _et al._ (2007); Phan  _et al._ (2018); Stawarz
_et al._ (2019) and in other solar, magnetospheric, planetary, astrophysical,
and laboratory settings Zweibel and Yamada (2009), 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 non-thermal energy at a relatively
sharp boundary layer Marcowith  _et al._ (2016). 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 Jaeger and Liu (2010), open quantum systems
Schaller (2014), and molecular dynamics in chemistry and biology Evans and
Hoover (1986); Karplus and Petsko (1990).

We expect $P_{\sigma,{\rm HORNET}}$ to be useful in a variety of analysis
techniques for non-LTE systems for which $f_{\sigma}$ 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_{\sigma,{\rm HORNET}}$ is possible. Slices of phase space densities can be
measured in plasma experiments Raymond  _et al._ (2018); Milder  _et al._
(2021a, b); Shi  _et al._ (2022a, b), although full three-dimensional
measurements are currently beyond the present technology. For spacecraft data,
it has been shown Argall  _et al._ (2022); Lindberg  _et al._ (2022); Agapitov
_et al._ (2023) that kinetic entropy can reliably be measured using NASA’s
Magnetospheric Multiscale (MMS) satellites Burch  _et al._ (2016). One direct
avenue of future work is to measure $P_{\sigma,{\rm 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,{\rm 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,{\rm HORNET}}$ will change in
a system in which ions fully couple to large-scale physics Sharma Pyakurel
_et al._ (2019). 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,{\rm HORNET}}<0$ signature, analogous
to the negative $P_{e,{\rm 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_{\sigma,{\rm 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, $\mathrm{{Pi-D}}_{\sigma}$, 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., Liang  _et al._ (2020b)). The effects of physical
collisions and entropic approaches to quantifying numerical dissipation will
be topics of future study.

###### Acknowledgements.

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. We acknowledge basing our color table
for 2D plots on the “bam” color table Crameri (2020) in an effort to improve
color-vision deficiency accessibility Crameri  _et al._ (2020); Brewer (2020).
The colorbar extrema are green [$(r_{g},g_{g},b_{g})=(2,73,29)={\cal C}_{g}$]
and purple [$(r_{p},g_{p},b_{p})=(68,3,79)={\cal 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 ${\cal 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 ${\cal C}(k)=255-(255-{\cal
C}_{p})F(k)$ from $k=0$ to 127 (purple) and ${\cal C}(k)=255-(255-{\cal
C}_{g})F(-k)$ from $k=127$ to 254 (green), with $r(255)=g(255)=b(255)=0$
(black) and $F(k)=\tanh{\left[3(127-k)/127\right]}/\tanh{3}$.

##  Appendix A Why $P_{e,{\rm HORNET}}$ is negative inside the EDR near its
upstream edge during reconnection non-linear growth phase

In Sec. IV.1, we provided a physical reason for the negative $P_{e,{\rm
HORNET}}$ inside the EDR between the upstream edge and the region near the
X-line ($0.15\lesssim|y-y_{0}|\lesssim 0.25$) during the non-linear growth
phase of anti-parallel magnetic reconnection. Here, we provide further
details. We show reduced electron phase space densities $f_{e}(v_{x},v_{y})$
(with $v_{z}$ integrated out) in Fig. 5 at three locations from the
reconnection simulation from Ref. Barbhuiya and Cassak (2022) that is
identical to the simulation presented in this study except PPG = 25,600. Panel
(a) is at $(x-x_{0},y-y_{0})=(0,-0.3)$ inside the EDR, close to its upstream
edge. Here, $f_{e}$ is elongated in the parallel ($\approx\hat{x}$) direction,
as explained by Egedal and colleagues Egedal  _et al._ (2013). Panel (b) is at
$(x-x_{0},y-y_{0})=(0,-0.2),$ closer to the X-line inside the EDR where
$P_{e,{\rm HORNET}}<0$. Panel (c) is at $(x-x_{0},y-y_{0})=(0,-0.1),$ even
closer to the X-line where electrons undergo Speiser orbits. Panel (b) reveals
that $f_{e}$ contains particles from the distributions in both (a) and (c),
which reduces the anisotropy in $f_{e}$ going from panel (a) to panel (b).
Therefore, in the comoving reference frame for a plasma parcel convecting from
outside the EDR to inside the EDR close to its upstream edge, the phase space
density becomes more Maxwellian, so $P_{e,{\rm HORNET}}$ is negative, in
agreement with the simulations.

##  Appendix B Physical Mechanism for the Vector Heat Flux Density Structure
in the Reconnection Diffusion Region

Figure 7: Plot of $P_{e,{\rm HORNET}}$ identical to Fig. 3(a), except with
colorbar saturated at values identical to that of $P_{i,{\rm HORNET}}$ in Fig.
3(e).

At the X-line, the reduced phase space density is symmetric around $v_{x}=0$,
so $q_{e,x}=0$ (panel (d)). Similarly, $f_{e}$ is symmetric in $v_{y}$ (not
shown), so $q_{e,y}=0$ (panel (e)). However, the striated phase space density
$f_{e}$ at the X-line is not symmetric in $v_{z}$. Relative to the bulk flow
velocity $u_{e,z}$, the phase space density weighted by $v^{\prime 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 Ng  _et al._ (2011); Shuster  _et al._ (2015); Cassak and
Barbhuiya (2022); Egedal  _et al._ (2023). 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 vector heat flux $-\bm{\nabla}\cdot{\bf q}_{e}$ at and near the
X-line as seen by the green color near the X-line in Fig. 1(d).

##  Appendix C $P_{e,{\rm HORNET}}$ in the turbulence simulation

In Sec. V.1, we display a 2D plot of $P_{e,{\rm 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,{\rm HORNET}}$ in Fig. 3(e). The regions with non-zero
$P_{e,{\rm HORNET}}$ outside the intermittent structures are more clearly
visible here.

## References

  * Chapman and Cowling (1970) S. Chapman and T. G. Cowling,  _The mathematical theory of non-uniform gases. an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_ , 3rd ed. (Cambridge University Press, 1970).
  * Onsager (1931) L. Onsager, Phys. Rev. 37, 405 (1931).
  * Grad (1949) H. Grad, Communications on Pure and Applied Mathematics 2, 331 (1949).
  * Jou  _et al._ (2010) D. Jou, G. Lebon, and J. Casas-Vázquez,  _Extended Irreversible Thermodynamics_ , 4th ed. (Springer, Dordrecht, 2010).
  * Howes  _et al._ (2006) G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin, The Astrophysical Journal 651, 590 (2006).
  * Schekochihin  _et al._ (2009) A. Schekochihin, S. Cowley, W. Dorland, G. Hammett, G. G. Howes, E. Quataert, and T. Tatsuno, The Astrophysical Journal Supplement Series 182, 310 (2009).
  * Braginskii (1965) S. I. Braginskii, Reviews of Plasma Physics 1, 205 (1965).
  * Del Sarto  _et al._ (2016) D. Del Sarto, F. Pegoraro, and F. Califano, Physical Review E 93, 053203 (2016).
  * Yang  _et al._ (2017a) Y. Yang, W. H. Matthaeus, T. N. Parashar, P. Wu, M. Wan, Y. Shi, S. Chen, V. Roytershteyn, and W. Daughton, Physical Review E 95, 061201 (2017a).
  * Yang  _et al._ (2017b) Y. Yang, W. H. Matthaeus, T. N. Parashar, C. C. Haggerty, V. Roytershteyn, W. Daughton, M. Wan, Y. Shi, and S. Chen, Physics of Plasmas 24, 072306 (2017b), https://doi.org/10.1063/1.4990421 .
  * Yang  _et al._ (2022) Y. Yang, W. H. Matthaeus, S. Roy, V. Roytershteyn, T. N. Parashar, R. Bandyopadhyay, and M. Wan, The Astrophysical Journal 929, 142 (2022).
  * Cassak and Barbhuiya (2022) P. A. Cassak and M. H. Barbhuiya, Physics of Plasmas 29, 122306 (2022).
  * Cassak  _et al._ (2022) P. A. Cassak, M. H. Barbhuiya, and H. A. Weldon, Physics of Plasmas 29, 122307 (2022).
  * Barbhuiya and Cassak (2022) M. H. Barbhuiya and P. A. Cassak, Physics of Plasmas 29, 122308 (2022).
  * Zhou  _et al._ (2021) M. Zhou, H. Man, Y. Yang, Z. Zhong, and X. Deng, Geophysical Research Letters 48, e2021GL096372 (2021).
  * Song  _et al._ (2020) L. Song, M. Zhou, Y. Yi, X. Deng, Z. Zhong, and H. Man, Journal of Geophysical Research: Space Physics 125, e2020JA028278 (2020).
  * Du  _et al._ (2020) S. Du, G. P. Zank, X. Li, and F. Guo, Phys. Rev. E 101, 033208 (2020).
  * Chasapis  _et al._ (2018) A. Chasapis, Y. Yang, W. H. Matthaeus, T. N. Parashar, C. C. Haggerty, J. L. Burch, T. E. Moore, C. J. Pollock, J. Dorelli, D. J. Gershman, R. B. Torbert, and C. T. Russell, The Astrophysical Journal 862, 32 (2018).
  * Yang  _et al._ (2019) Y. Yang, M. Wan, W. H. Matthaeus, L. Sorriso-Valvo, T. N. Parashar, Q. Lu, Y. Shi, and S. Chen, Monthly Notices of the Royal Astronomical Society 482, 4933 (2019).
  * Pezzi  _et al._ (2019) O. Pezzi, Y. Yang, F. Valentini, S. Servidio, A. Chasapis, W. H. Matthaeus, and P. Veltri, Physics of Plasmas 26, 072301 (2019), https://doi.org/10.1063/1.5100125 .
  * Pezzi  _et al._ (2021) O. Pezzi, H. Liang, J. L. Juno, P. A. Cassak, C. L. Vásconez, L. Sorriso-Valvo, D. Perrone, S. Servidio, V. Roytershteyn, J. M. TenBarge, and W. H. Matthaeus, MNRAS 505, 4857 (2021).
  * Bandyopadhyay  _et al._ (2021) R. Bandyopadhyay, A. Chasapis, W. H. Matthaeus, T. N. Parashar, C. C. Haggerty, M. A. Shay, D. J. Gershman, B. L. Giles, and J. L. Burch, Physics of Plasmas 28, 112305 (2021).
  * Klein and Howes (2016) K. G. Klein and G. G. Howes, The Astrophysical Journal Letters 826, L30 (2016).
  * Howes (2017) G. G. Howes, Physics of Plasmas 24, 055907 (2017), https://doi.org/10.1063/1.4983993 .
  * Klein  _et al._ (2017) K. G. Klein, G. G. Howes, and J. M. TenBarge, Journal of Plasma Physics 83 (2017).
  * Chen  _et al._ (2019) C. H. K. Chen, K. G. Klein, and G. G. Howes, Nature Communications 10, 740 (2019).
  * Klein  _et al._ (2020) K. G. Klein, G. G. Howes, J. M. TenBarge, and F. Valentini, Journal of Plasma Physics 86 (2020).
  * Juno  _et al._ (2021) J. Juno, G. G. Howes, J. M. TenBarge, L. B. Wilson, A. Spitkovsky, D. Caprioli, K. G. Klein, and A. Hakim, Journal of Plasma Physics 87 (2021).
  * McCubbin  _et al._ (2022) A. J. McCubbin, G. G. Howes, and J. M. TenBarge, Physics of plasmas 29 (2022).
  * Montag and Howes (2022) P. Montag and G. G. Howes, Physics of Plasmas 29, 032901 (2022).
  * Montag  _et al._ (2023) P. Montag, G. Howes, D. McGinnis, A. Afshari, M. Starkey, and M. Desai, “Mms observations of the velocity-space signature of shock-drift acceleration,” (2023), arXiv:2306.09061 .
  * Cassak  _et al._ (2023) P. A. Cassak, M. H. Barbhuiya, H. Liang, and M. R. Argall, Physical Review Letters 130, 085201 (2023).
  * Liang  _et al._ (2019) H. Liang, P. A. Cassak, S. Servidio, M. A. Shay, J. F. Drake, M. Swisdak, M. R. Argall, J. C. Dorelli, E. E. Scime, W. H. Matthaeus, V. Roytershteyn, and G. L. Delzanno, Phys. Plasmas 26, 082903 (2019).
  * Boltzmann (1872) L. Boltzmann, Wiener Berichte 66, 275 (1872), in (Boltzmann 1909) Vol. I, paper 23.
  * Liang  _et al._ (2020a) H. Liang, M. H. Barbhuiya, P. A. Cassak, O. Pezzi, S. Servidio, F. Valentini, and G. P. Zank, Journal of Plasma Physics 86, 825860502 (2020a).
  * Argall  _et al._ (2022) M. R. Argall, M. H. Barbhuiya, P. A. Cassak, S. Wang, J. Shuster, H. Liang, D. J. Gershman, R. B. Torbert, and J. L. Burch, Physics of Plasmas 29, 022902 (2022).
  * Grad (1965) H. Grad, Journal of the Society for Industrial and Applied Mathematics 13, 259 (1965).
  * Eu (1995) B. C. Eu, The Journal of Chemical Physics 102, 7169 (1995).
  * Zeiler  _et al._ (2002) A. Zeiler, D. Biskamp, J. F. Drake, B. N. Rogers, M. A. Shay, and M. Scholer, J. Geophys. Res. 107, 1230 (2002).
  * Birdsall and Langdon (1991) C. K. Birdsall and A. B. Langdon, “Plasma physics via computer simulation,” (Institute of Physics Publishing, Philadelphia, 1991) Chap. 15.
  * Guzdar  _et al._ (1993) P. N. Guzdar, J. F. Drake, D. McCarthy, A. B. Hassam, and C. S. Liu, Phys. Fluids B 5, 3712 (1993).
  * Trottenberg  _et al._ (2000) U. Trottenberg, C. W. Oosterlee, and A. Schuller,  _Multigrid_ (Academic Press, San Diego, 2000).
  * Liang  _et al._ (2020b) H. Liang, P. A. Cassak, M. Swisdak, and S. Servidio, Journal of Physics: Conference Series 1620, 012009 (2020b).
  * Parashar  _et al._ (2018) T. N. Parashar, W. H. Matthaeus, and M. A. Shay, The Astrophysical Journal Letters 864, L21 (2018).
  * Egedal  _et al._ (2013) J. Egedal, A. Le, and W. Daughton, Physics of Plasmas 20, 061201 (2013).
  * Speiser (1965) T. W. Speiser, Journal of Geophysical Research (1896-1977) 70, 4219 (1965).
  * Ng  _et al._ (2011) J. Ng, J. Egedal, A. Le, W. Daughton, and L.-J. Chen, Physical Review Letters 106, 065002 (2011).
  * Shuster  _et al._ (2015) J. R. Shuster, L.-J. Chen, M. Hesse, M. R. Argall, W. Daughton, R. B. Torbert, and N. Bessho, Geophysical Research Letters 42, 2586 (2015).
  * Egedal  _et al._ (2012) J. Egedal, W. Daughton, and A. Le, Nature Physics 8 (2012).
  * Lavraud  _et al._ (2016) B. Lavraud, Y. C. Zhang, Y. Vernisse, D. J. Gershman, J. Dorelli, P. A. Cassak, J. Dargent, C. Pollock, B. Giles, N. Aunai, M. Argall, L. Avanov, A. Barrie, J. Burch, M. Chandler, L.-J. Chen, G. Clark, I. Cohen, V. Coffey, J. P. Eastwood, J. Egedal, S. Eriksson, R. Ergun, C. J. Farrugia, S. A. Fuselier, V. Génot, D. Graham, E. Grigorenko, H. Hasegawa, C. Jacquey, I. Kacem, Y. Khotyaintsev, E. MacDonald, W. Magnes, A. Marchaudon, B. Mauk, T. E. Moore, T. Mukai, R. Nakamura, W. Paterson, E. Penou, T. D. Phan, A. Rager, A. Retino, Z. J. Rong, C. T. Russell, Y. Saito, J.-A. Sauvaud, S. J. Schwartz, C. Shen, S. Smith, R. Strangeway, S. Toledo-Redondo, R. Torbert, D. L. Turner, S. Wang, and S. Yokota, Geophysical Research Letters 43, 3042 (2016).
  * Wang  _et al._ (2016) S. Wang, L. Chen, N. Bessho, L. M. Kistler, J. R. Shuster, and R. Guo, Journal of Geophysical Research: Space Physics 121, 2104 (2016).
  * Norgren  _et al._ (2020) C. Norgren, M. Hesse, D. B. Graham, Y. V. Khotyaintsev, P. Tenfjord, A. Vaivads, K. Steinvall, Y. Xu, D. J. Gershman, P.-A. Lindqvist, F. Plaschke, and J. L. Burch, Journal of Geophysical Research: Space Physics 125, e2019JA027440 (2020), e2019JA027440 10.1029/2019JA027440.
  * Sharma Pyakurel  _et al._ (2019) P. Sharma Pyakurel, M. A. Shay, T. D. Phan, W. H. Matthaeus, J. F. Drake, J. M. TenBarge, C. C. Haggerty, K. G. Klein, P. A. Cassak, T. N. Parashar, M. Swisdak, and A. Chasapis, Physics of Plasmas 26, 082307 (2019), 10.1063/1.5090403 .
  * Eastwood  _et al._ (2020) J. P. Eastwood, M. V. Goldman, T. D. Phan, J. E. Stawarz, P. A. Cassak, J. F. Drake, D. Newman, B. Lavraud, M. A. Shay, R. E. Ergun, J. L. Burch, D. J. Gershman, B. L. Giles, P. A. Lindqvist, R. B. Torbert, R. J. Strangeway, and C. T. Russell, Physical Review Letters 125, 265102 (2020).
  * Servidio  _et al._ (2009) S. Servidio, W. H. Matthaeus, M. A. Shay, P. A. Cassak, and P. Dmitruk, Phys. Rev. Lett. 102, 115003 (2009).
  * Agudelo Rueda  _et al._ (2022) J. A. Agudelo Rueda, D. Verscharen, R. T. Wicks, C. J. Owen, G. Nicolaou, K. Germaschewski, A. P. Walsh, I. Zouganelis, and S. V. Domínguez, The Astrophysical Journal 938, 4 (2022).
  * Matthaeus (2021) W. Matthaeus, Physics of Plasmas 28 (2021).
  * Wan  _et al._ (2012) M. Wan, W. H. Matthaeus, H. Karimabadi, V. Roytershteyn, M. Shay, P. Wu, W. Daughton, B. Loring, and S. C. Chapman, Physical Review Letters 109, 195001 (2012).
  * Matthaeus  _et al._ (2020) W. H. Matthaeus, Y. Yang, M. Wan, T. N. Parashar, R. Bandyopadhyay, A. Chasapis, O. Pezzi, and F. Valentini, The Astrophysical Journal 891, 101 (2020).
  * Matthaeus and Lamkin (1986) W. H. Matthaeus and S. L. Lamkin, Physics of Fluids 29, 2513 (1986).
  * Servidio  _et al._ (2010) S. Servidio, W. H. Matthaeus, M. A. Shay, P. Dmitruk, P. A. Cassak, and M. Wan, Physics of Plasmas 17, 032315 (2010), https://doi.org/10.1063/1.3368798 .
  * Haggerty  _et al._ (2017) C. C. Haggerty, T. N. Parashar, W. H. Matthaeus, M. A. Shay, Y. Yang, M. Wan, P. Wu, and S. Servidio, Physics of Plasmas 24, 102308 (2017).
  * Retinò  _et al._ (2007) A. Retinò, D. Sundkvist, A. Vaivads, F. Mozer, M. André, and C. J. Owen, Nature Physics 3, 235 (2007).
  * Sundkvist  _et al._ (2007) D. Sundkvist, A. Retinò, A. Vaivads, and S. D. Bale, Physical Review Letters 99, 025004 (2007).
  * Phan  _et al._ (2018) T. Phan, J. Eastwood, M. Shay, J. Drake, B. Sonnerup, M. Fujimoto, P. Cassak, M. Øieroset, J. Burch, R. Torbert, A. Rager, J. Dorelli, D. Gershman, C. Pollock, P. Pyakurel, C. Haggerty, Y. Khotyaintsev, B. Lavraud, Y. Saito, M. Oka, R. Ergun, A. Retino, O. Le Contel, M. Argall, B. Giles, T. Moore, F. Wilder, R. Strangeway, C. Russell, P. Lindqvist, and W. Magnes, Nature 557, 202 (2018).
  * Stawarz  _et al._ (2019) J. E. Stawarz, J. P. Eastwood, T. D. Phan, I. L. Gingell, M. A. Shay, J. L. Burch, R. E. Ergun, B. L. Giles, D. J. Gershman, O. Le Contel, P.-A. Lindqvist, C. T. Russell, R. J. Strangeway, R. B. Torbert, M. R. Argall, D. Fischer, W. Magnes, and L. Franci, Astrophysical Journal Letters 877, L37 (2019).
  * Zweibel and Yamada (2009) E. G. Zweibel and M. Yamada, Annu. Rev. Astron. Astrophys. 47, 291 (2009).
  * Marcowith  _et al._ (2016) A. Marcowith, A. Bret, A. Bykov, M. E. Dieckman, L. O. Drury, B. Lembège, M. Lemoine, G. Morlino, G. Murphy, G. Pelletier, I. Plotnikov, B. Reville, M. Riquelme, L. Sironi, and A. S. Novo, Rep. Prog. Phys. 79 (2016).
  * Jaeger and Liu (2010) H. M. Jaeger and A. J. Liu, arXiv preprint arXiv:1009.4874 (2010).
  * Schaller (2014) G. Schaller,  _Open Quantum Systems Far from Equilibrium_ (Springer International Publishing, 2014).
  * Evans and Hoover (1986) D. J. Evans and W. G. Hoover, Annual Review of Fluid Mechanics 18, 243 (1986), 10.1146/annurev.fl.18.010186.001331 .
  * Karplus and Petsko (1990) M. Karplus and G. A. Petsko, Nature 347, 631 (1990).
  * Raymond  _et al._ (2018) A. E. Raymond, C. F. Dong, A. McKelvey, C. Zulick, N. Alexander, A. Bhattacharjee, P. T. Campbell, H. Chen, V. Chvykov, E. Del Rio, P. Fitzsimmons, W. Fox, B. Hou, A. Maksimchuk, C. Mileham, J. Nees, P. M. Nilson, C. Stoeckl, A. G. R. Thomas, M. S. Wei, V. Yanovsky, K. Krushelnick, and L. Willingale, Phys. Rev. E 98, 043207 (2018).
  * Milder  _et al._ (2021a) A. L. Milder, J. Katz, R. Boni, J. P. Palastro, M. Sherlock, W. Rozmus, and D. H. Froula, Physics of Plasmas 28, 082102 (2021a).
  * Milder  _et al._ (2021b) A. Milder, J. Katz, R. Boni, J. Palastro, M. Sherlock, W. Rozmus, and D. Froula, Physical Review Letters 127, 015001 (2021b).
  * Shi  _et al._ (2022a) P. Shi, P. Srivastav, M. H. Barbhuiya, P. A. Cassak, E. E. Scime, M. Swisdak, C. Beatty, T. Gilbert, R. John, M. Lazo, R. S. Nirwan, M. Paul, E. E. Scime, K. Stevenson, and T. Steinberger, Physics of Plasmas 29, 032101 (2022a), publisher: American Institute of Physics.
  * Shi  _et al._ (2022b) P. Shi, P. Srivastav, M. H. Barbhuiya, P. A. Cassak, E. E. Scime, and M. Swisdak, Physical Review Letters 128, 025002 (2022b).
  * Lindberg  _et al._ (2022) M. Lindberg, A. Vaivads, S. Raptis, P.-A. Lindqvist, B. L. Giles, and D. J. Gershman, Entropy 24, 745 (2022), number: 6 Publisher: Multidisciplinary Digital Publishing Institute.
  * Agapitov  _et al._ (2023) O. V. Agapitov, V. Krasnoselskikh, M. Balikhin, J. W. Bonnell, F. S. Mozer, and L. Avanov, “Energy Repartition and Entropy Generation across the Earth’s Bow Shock: MMS Observations,”  (2023), arXiv:2306.03194.
  * Burch  _et al._ (2016) J. L. Burch, T. E. Moore, R. B. Torbert, and B. L. Giles, Space Sci. Rev. 199, 5 (2016).
  * Crameri (2020) F. Crameri, “Scientific colour maps,” (2020), http://www.fabiocrameri.ch/colourmaps .
  * Crameri  _et al._ (2020) F. Crameri, G. E. Shephard, and P. J. Heron, Nature Communications 11, 5444 (2020).
  * Brewer (2020) C. A. Brewer, “Colorbrewer colour maps,” (2020), https://colorbrewer2.org/ .
  * Egedal  _et al._ (2023) J. Egedal, H. Gurram, S. Greess, W. Daughton, and A. Lê, Physics of Plasmas 30, 062106 (2023).