††thanks: Corresponding author.
<EMAIL_ADDRESS>
# Intermittency and electron heating in kinetic-Alfvén-wave turbulence

Muni Zhou  Plasma Science and Fusion Center, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Zhuo Liu  Plasma Science and Fusion
Center, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Nuno
F. Loureiro  Plasma Science and Fusion Center, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA

(August 27, 2024)

###### Abstract

We report analytical and numerical investigations of sub-ion-scale turbulence
in low-beta plasmas, focusing on the spectral properties of the fluctuations
and electron heating. In the isothermal limit, the numerical results strongly
support a description of the turbulence as a critically-balanced Kolmogorov-
like cascade of kinetic Alfvén wave fluctuations, as amended by Boldyrev &
Perez (Astrophys. J. Lett. 758, L44 (2012)) to include intermittent effects.
When the constraint of isothermality is removed (i.e., with the inclusion of
electron kinetic physics), the energy spectrum is found to steepen due to
electron Landau damping, which is enabled by the local weakening of advective
nonlinearities around current sheets, and yields significant energy
dissipation via a velocity-space cascade. The use of a Hermite-polynomial
representation to express the velocity-space dependence of the electron
distribution function allows us to obtain an analytical, lowest-order solution
for the Hermite moments of the distribution, which is borne out by numerical
simulations.

##  I Introduction

Plasma turbulence is ubiquitous in space and astrophysical systems as diverse
as the Earth’s magnetosphere [borovsky2003mhd], the solar wind [Goldreich1995,
tu1995mhd], the solar corona [parker1983], accretion disks
[balbus1998instability], and the interstellar and intracluster media
[begelman1990turbulent]. Many of these environments are sufficiently dilute
that particle collisions are rare on the dynamical timescales of interest. As
such, kinetic plasma descriptions are needed to understand turbulence in those
environments. Despite decades of intensive research, formulating a predictive
theoretical framework for the dynamics of turbulence in the phase space of
positions and velocities has been a notoriously difficult problem, especially
in the sub-ion scale range (the kinetic range), where various microphysical
plasma processes are dynamically important.

Recent high-resolution, in situ measurements of electromagnetic fluctuations
and plasma distribution functions from satellites such as the Magnetospheric
Multiscale mission (MMS), Cluster, and the Parker Solar Probe (PSP) have
provided unprecedented opportunities to study the rich plasma dynamics in the
sub-ion range of the turbulence [sahraoui2009evidence,
alexandrova2009universality, alexandrova2012solar, alexandrova2013solar,
chen2019evidence]. The interpretation and understanding of those observations
essentially reduces to answering two questions: (i) What are the spectral
properties of the electromagnetic fluctuations, and more specifically, what
physical processes determine the energy spectra in the kinetic range? (ii)
What is the nature of the dissipative processes, i.e., how are ions and
electrons heated? Answering these two questions requires a comprehensive
understanding of the phase-space dynamics of kinetic plasma turbulence.

Alfvén waves (AWs), an essential building block of long-wavelength
magnetohydrodynamic (MHD) turbulence, convert into (dispersive) kinetic Alfvén
waves (KAWs) in the sub-ion range. The KAW nature of kinetic turbulence has
been supported by a substantial amount of analytical and numerical work
[gary2004kinetic, howes2008model, howes2008kinetic, salem2012identification,
tenbarge2012interpreting, franci2018solar, cerri2019kinetic,
grovselj2019kinetic], as well as observations [leamon1998observational,
sahraoui2009evidence, podesta2010kinetic, chen2013nature, kiyani2012enhanced,
chen2016recent]. Therefore, theoretical models conjecture a critically
balanced [Goldreich1995] cascade of KAWs, arriving at predictions of a $\sim
k_{\perp}^{-7/3}$ magnetic energy spectrum [cho2004anisotropy,
cho2009simulations, howes2008model, schekochihin2009astrophysical]. However, a
large number of observations [leamon1998observational, bale2005measurement,
sahraoui2009evidence, salem2012identification, chen2013nature] and numerical
simulations [tenbarge2012interpreting, servidio2015kinetic,
wan2015intermittent, told2015multiscale, cerri2018dual, arzamasskiy2019hybrid,
grovselj2018fully] yield steeper spectra, typically with a spectral index
close to $-2.8$ down to electron scales (e.g., [alexandrova2009universality,
kiyani2009global, sahraoui2013scaling]).

One important attempt to reconcile theory with observations invokes
intermittency. Indeed, several studies suggest a positive correlation between
intermittency and energy dissipation [Sundkvist2007, karimabadi2011flux,
Osman2014, wan2015, chasapis2015thin, Camporeale2018] in the kinetic range;
and a phenomenological model has been proposed to describe the role of
intermittency in the energy cascade [boldyrev2012spectrum]: assuming that the
fluctuations, and thus their energy cascade, are concentrated in two-
dimensional (2D) structures, a $\propto k_{\perp}^{-8/3}$ magnetic energy
spectrum can be derived, which to some extent explains the steepening of the
spectrum (without invoking electron kinetic effects) with respect to a
standard KAW-cascade model.

A distinct possibility that has been mooted is that of the KAW cascade being
mediated by the tearing mode instability. This hypothesis is inspired by
similar thinking within the scope of MHD turbulence. There, the assumption of
progressive alignment of the turbulent fluctuations towards smaller scales
[boldyrev2006spectrum] has led to the theoretical conjecture (and tentative
numerical confirmation) that, at sufficiently small scales, turbulent eddies
should possess a large aspect ratio and become susceptible to tearing,
yielding a dynamical range where the energy cascade is mediated by that
instability [loureiro2017role, mallet2017statistical,
boldyrev2017magnetohydrodynamic, dong2018role, walker2018influence,
loureiro2020nonlinear]. Extending these ideas to the sub-ion range is non-
trivial because of complications caused by various kinetic effects and the
lack of analytical models that predict the scale-dependent anisotropy of
turbulent fluctuations in this range. However, if tearing mediation were to
exist in the sub-ion range, a $k_{\perp}^{-8/3}$ (or $k_{\perp}^{-3}$,
depending on the choice of local magnetic configuration) magnetic spectrum is
predicted [loureiro2017collisionless], which is also consistent with the
steepening of the spectrum. It is an unfortunate coincidence that this
prediction is numerically the same as that arising from boldyrev2012spectrum
intermittency model: the underlying physics is entirely different.

Apart from understanding the spatial properties of the fluctuations, one of
the ultimate goals in studying kinetic turbulence is to understand how energy
is dissipated in weakly-collisional plasmas and how the electrons and ions are
energized [parashar2015turbulent]. In a collisional plasma, energy can only be
thermalized through the “fluid-channel”, in which energy cascades to small
spatial scales through nonlinear advection and finally dissipates by viscosity
and resistivity [dmitruk2004test, Matthaeus2011who, servidio2011, osman2011,
zhdankin2013statistical], while the particle distribution remains close to
local thermodynamic equilibrium [chapman1990mathematical]. In a weakly
collisional plasma, however, particle free streaming along magnetic field
lines leads to the “phase mixing” process, which smooths out the
electromagnetic fluctuations and develops complex structures in velocity
space. This kinetic effect enables an additional “kinetic-channel” for
dissipation, in which energy transfers to small scales in velocity space
(i.e., high velocity moments of the distribution function) and dissipates
through collisions [howes2008model, schekochihin2008gyrokinetic,
navarro2016structure, grovselj2017fully, servidio2017magnetospheric]. Energy
dissipation in kinetic turbulence occurs via a combination of these two
channels — their relative importance is determined by the complex phase-space
dynamics of kinetic turbulence [schekochihin2016phase, adkins2018solvable,
cerri2018dual, eyink2018cascades].

The kinetic dissipation channel is, in principle, susceptible to the plasma
echo effect [Gould1967echo, malmberg1968plasma], whereby free energy inverse
cascades in velocity space, returning from the high to the low moments of the
distribution function (i.e., phase unmixing). The occurrence of a collection
of the stochastic plasma echoes can cause kinetic turbulence to resemble fluid
turbulence, as the only allowed energy-cascade channel in that case is the
fluid-type cascade towards smaller spatial scales [schekochihin2016phase,
adkins2018solvable]. The significance of this effect has been confirmed by
numerical simulations of collisionless plasma turbulence at scales above the
ion Larmor radius ($\rho_{i}$) with isothermal electrons focusing on
compressive fluctuations [meyrand2019fluidization], and below $\rho_{i}$ with
an electrostatic drift-kinetic model [parker2016suppression]. However, how the
possible occurrence of plasma echoes can be reconciled with observations of
efficient electron and ion heating (e.g., [servidio2017magnetospheric,
chen2019evidence]) at kinetic scales is an open question.

In this work, we provide answers to the aforementioned two questions — what
the physical mechanism that sets the energy spectrum is, and how electron
heating occurs — in the low-$\beta$ limit of sub-$\rho_{i}$ turbulence. The
theoretical framework we adopt is described in Sec. II. In Sec. III, we first
review three leading models for the KAW spectrum in the isothermal limit, and
then, with electron kinetic effects accounted for, derive a (nonlinear)
lowest-order solution of the electron distribution function in the velocity
space. We then test these models using numerical simulations, whose details
are provided in Sec. IV and from which the results presented in Sec. V are
obtained. Sec. VI presents our conclusions and discussion of our results.

##  II Theoretical framework

Our results in this paper are obtained with an analytical framework known as
the Kinetic Reduced Electron Heating Model (KREHM) [zocco2011reduced]. KREHM
is a rigorous asymptotic reduction of gyrokinetics valid in the limit of low
electron plasma-beta, $\beta_{e}\sim m_{e}/m_{i}$, with $\beta_{e}\equiv 8\pi
n_{0e}T_{0e}/B_{0}^{2}$, where $B_{0}$ is the background magnetic (guide)
field strength (assumed constant and directed along $\hat{z}$), and $n_{0e}$,
$T_{0e}$ are the background electron density and temperature, respectively;
$m_{e}$ and $m_{i}$ denote the electron and ion masses. Within this limit, the
parallel streaming and electromagnetic effects are ordered out in the ion
gyrokinetic equation. It follows that the ion flow velocity parallel to the
background field is $u_{zi}=0$, and ions become isothermal and electrostatic.
Ion finite Larmor radius (FLR) effects are retained in the gyrokinetic
Poisson’s law,

$\frac{\delta
n_{e}}{n_{0e}}=\frac{1}{\tau}(\hat{\Gamma}_{0}-1)\frac{e\varphi}{T_{0e}},$ (1)

where $\varphi$ is the electrostatic potential, $\delta n_{e}/n_{0e}$ is the
electron density perturbation normalized to its background value, $\tau\equiv
T_{0i}/T_{0e}$ is the temperature ratio, and $\hat{\Gamma}_{0}$ is a
gyroaveraging operator that has the closed-form expression
$\Gamma_{0}(\alpha)=I_{0}(\alpha)e^{-\alpha}$ in Fourier space; here, $I_{0}$
is the zeroth-order modified Bessel function of the first kind and
$\alpha=k_{\perp}^{2}\rho_{i}^{2}/2$, where $\rho_{i}=v_{{\rm
th}i}/\Omega_{i}$ is the ion Larmor radius, $v_{{\rm
th}i}=\sqrt{2T_{0i}/m_{i}}$ is the ion thermal velocity, and
$\Omega_{i}=eB_{0}/m_{i}c$ is the ion Larmor frequency (we consider single
charge ions in this work, $Z=1$).

The electrons are described by a perturbed distribution function which, to
order $\sqrt{m_{e}/m_{i}}\sim\sqrt{\beta_{e}}$ in the gyrokinetic expansion,
can be written as $\delta f_{e}=g_{e}+(\delta
n_{e}/n_{0e}+2v_{z}u_{ze}/v_{{\rm th}e}^{2})F_{0e}$, where $F_{0e}$ is the
equilibrium Maxwellian defined with the mean electron temperature $T_{0e}$ and
its corresponding thermal speed $v_{{\rm th}e}=\sqrt{2T_{0e}/m_{e}}$, and
$\delta n_{e}$ and $u_{ze}$ (electron flow parallel to the guide field) are
the zeroth and first moments of $\delta f_{e}$, respectively. Since
$u_{zi}=0$, the parallel component of Ampère’s law,
$J_{z}=-(c/4\pi)\nabla_{\perp}^{2}A_{z}$, where $A_{z}$ is the parallel
component of the vector potential, leads to
$u_{ze}=(e/cm_{e})d_{e}^{2}\nabla_{\perp}^{2}A_{z}$, where
$d_{e}=c/\omega_{pe}$ is the electron skin depth. Information about the second
and higher moments of $\delta f_{e}$ is contained in the (reduced)
distribution function $g_{e}$.

The dynamics of this system are described by fluid equations for the first two
moments, coupled to a drift-kinetic equation for $g_{e}$:

$\frac{1}{n_{0e}}\frac{d\delta
n_{e}}{dt}=-\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\frac{e}{cm_{e}}d_{e}^{2}\nabla_{\perp}^{2}A_{z},$
(2)
$\frac{d}{dt}(A_{z}-d_{e}^{2}\nabla_{\perp}^{2}A_{z})=-c\frac{\partial\varphi}{\partial
z}+\frac{cT_{e0}}{e}\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\left(\frac{\delta
n_{e}}{n_{0e}}+\frac{\delta T_{ze}}{T_{0e}}\right),$ (3)
$\begin{split}&\frac{dg_{e}}{dt}+v_{z}\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\left(g_{e}-\frac{\delta
T_{ze}}{T_{0e}}F_{0e}\right)\\\ &=C[g_{e}]+\left(1-\frac{v_{z}^{2}}{v_{{\rm
th}e}^{2}}\right)F_{0e}\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\frac{e}{cm_{e}}d_{e}^{2}\nabla_{\perp}^{2}A_{z}.\end{split}$
(4)

In these equations, $d/dt\equiv\partial/\partial t+c/B_{0}\\{\varphi,...\\}$
denotes the convective time derivative, with the Poisson bracket defined as
$\\{P,Q\\}\equiv\partial_{x}P\partial_{y}Q-\partial_{y}P\partial_{x}Q$. The
parallel (to the total field) gradient operator is
$\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}P\equiv\partial_{z}P-\\{A_{z},P\\}/B_{0}$.
The last term in the generalized Ohm’s law features the (normalized) parallel
electron temperature perturbation, defined in terms of the reduced electron
distribution function as

$\frac{\delta T_{ze}}{T_{0e}}=\frac{1}{n_{0e}}\left(\int
d^{3}{\mbox{\boldmath{$v$}}}\frac{2v_{z}^{2}}{v_{{\rm th}e}^{2}}g_{e}\right);$
(5)

the isothermal limit of these equations corresponds to $\delta T_{ze}=0$, that
is, $g_{e}=0$. Lastly, $C[g_{e}]$ represents the collision operator. Since Eq.
(4) has no explicit dependence on $v_{\perp}$, this coordinate can be
integrated out of the problem if a collision operator is chosen that bares no
such dependence, such as the (modified) Lenard-Bernstein collision operator
[zocco2011reduced]. In that case, the reduced distribution function becomes
four-dimensional, $g_{e}=g_{e}(\mbox{\boldmath{$r$}},v_{z},t)$.

In the absence of collisions, Eqs. (2)-(4) conserve a quadratic invariant that
is usually referred to as the total free energy

$\displaystyle W$
$\displaystyle\equiv\int\frac{d^{3}\mbox{\boldmath{$r$}}}{V}\left[\frac{1}{\tau}(1-\hat{\Gamma}_{0})\frac{e^{2}n_{0e}\varphi^{2}}{2T_{0e}}+\frac{|\nabla_{\perp}A_{z}|^{2}}{8\pi}\right.$
(6)
$\displaystyle+\left.\frac{1}{\tau^{2}}(1-\hat{\Gamma}_{0})^{2}\frac{e^{2}n_{0e}\varphi^{2}}{2T_{0e}}+\frac{d_{e}^{2}|\nabla^{2}_{\perp}A_{z}|^{2}}{8\pi}+\int
d^{3}\mbox{\boldmath{$v$}}\frac{T_{0e}g_{e}^{2}}{2F_{0e}}\right],$

where the different terms on the right-hand side, in the order in which they
appear, correspond to the ion perturbed entropy, the magnetic energy, the
electron density variance, the kinetic energy of the parallel electron flow,
and the electron free energy, respectively. For completeness, we note that the
isothermal ($g_{e}=0$) limit of these equations admits another ideal invariant
— the generalized helicity,

$H\equiv\int\frac{(1-\hat{\Gamma}_{0})\varphi}{\rho_{i}^{2}}(1-d_{e}^{2}\nabla_{\perp}^{2})A_{z}dV,$
(7)

which reduces to $H=\int n_{e}A_{z}dV$ in the range $\rho_{i}\gg
k_{\perp}^{-1}\gg d_{e}$, and to the cross-helicity at MHD scales
($k_{\perp}\rho_{i}\ll 1$).

To better describe the velocity-space dynamics, we follow Ref.
[zocco2011reduced] and expand $g_{e}$ in Hermite polynomials:

$g_{e}(\mbox{\boldmath{$r$}},v_{z},t)=\frac{1}{\sqrt{2^{m}m!}}\sum_{m=0}^{\infty}H_{m}\left(\frac{v_{z}}{v_{{\rm
th}e}}\right)g_{m}(\mbox{\boldmath{$r$}},t)F_{0e}(v_{z}),$ (8)

where

$g_{m}=\frac{1}{n_{0e}}\frac{1}{\sqrt{2^{m}m!}}\int_{\infty}^{\infty}dv_{z}H_{m}\left(\frac{v_{z}}{v_{{\rm
th}e}}\right)g_{e}.$ (9)

Since $g_{e}$ only contains moments higher than $\delta n_{e}$ and $u_{ze}$,
it follows that $g_{0}=g_{1}=0$. Using the Hermite expansion, Eq. (4) can be
decomposed into a series of coupled equations for the Hermite coefficients
$g_{m}$ (for $m\geq 2$):

$\begin{split}\frac{dg_{m}}{dt}&=-v_{{\rm
th}e}\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\left(\sqrt{\frac{m+1}{2}}g_{m+1}+\sqrt{\frac{m}{2}}g_{m-1}-\delta_{m,1}g_{2}\right)\\\
&-\sqrt{2}\delta_{m,2}\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}\frac{e}{cm_{e}}d_{e}^{2}\nabla_{\perp}^{2}A_{z}-C[g_{m}].\end{split}$
(10)

At large $m$, free energy is converted to electron entropy and thence to heat
via the collisional operator $C[g_{m}]$. In summary, the KREHM framework
captures ion FLR effects and electron (drift) kinetic dynamics, including
linear Landau damping. It is accurate at all scales above the electron Larmor
radius, under the low-beta constraint $\beta_{e}\sim m_{e}/m_{i}$. The only
linear mode captured by these equations is the kinetic Alfvén wave (KAW); as
such, it provides the simplest-possible fully kinetic platform to study
turbulence below the ion scales. In the long-wavelength limit
$k_{\perp}\rho_{i}\ll 1$, the KREHM equations reduce to the reduced-
magnetohydrodynamic (RMHD) equations [kadomtsev1973nonlinear,
strauss1976nonlinear, zank1992, schekochihin2009astrophysical].

##  III Theoretical models of Kinetic-Alfvén-wave turbulence

In this section, we first focus on the $\rho_{i}\gg k_{\perp}^{-1}\gg d_{e}$
scale range and discuss the spectral properties of the fluctuations in
position and velocity space. The position-space spectra is discussed in the
limit of isothermal electrons, i.e., the limit when the perturbed electron
distribution vanishes, $g_{e}=0$. In the case of $g_{e}\neq 0$, we obtain a
lowest-order solution for $g_{e}$ in velocity space, and discuss its
implications for electron heating. We then extend our discussion to the sub-
electron scale range, $k_{\perp}d_{e}\gtrsim 1$.

###  III.1 Spectra and spectral anisotropy

At sub-$\rho_{i}$ scales, the turbulent fluctuations are composed of
nonlinearly interacting, strongly anisotropic KAWs propagating along the
“local mean” magnetic field. The total magnetic field is
$\mbox{\boldmath{$B$}}=B_{0}\mbox{\boldmath{$\hat{z}$}}+\delta\mbox{\boldmath{$B$}}_{\perp}$,
where the perpendicular (fluctuating) magnetic field is
$\delta\mbox{\boldmath{$B$}}_{\perp}=-\hat{\mbox{\boldmath{$z$}}}\times\mbox{\boldmath{$\nabla$}}A_{z}$,
and $\delta B_{\perp}/B_{0}\ll 1$. The wave numbers of the fluctuations can be
decomposed into parallel and perpendicular components, $k_{\parallel}$ and
$k_{\perp}$, with respect to $B$. The multi-scale fluctuations can be
characterized as turbulent eddies with the correlation length $\ell\sim
1/k_{\parallel}$ in the parallel (to the local magnetic field) direction, the
width $\xi$ in the direction of the local perpendicular fluctuation, and the
thickness $\lambda\sim 1/k_{\perp}$ in the direction perpendicular to the
other two. In the rest of the paper, we use the subscript $\lambda$ to refer
to the amplitude of fluctuating quantities at scale $\lambda\sim
k_{\perp}^{-1}$.

First, assuming scale-by-scale equipartition between the density and magnetic
energy fluctuations, $(\delta
n_{e}/n_{0e})^{2}n_{0e}T_{0e}\sim|\nabla_{\perp}A_{z}|^{2}/8\pi$, leads to the
relation

$\varphi_{\lambda}\sim(\rho_{i}V_{A}/c)\delta B_{\perp\lambda},$ (11)

with $V_{A}=B_{0}/\sqrt{4\pi\rho_{0}}$ the Alfvén speed pertaining to the
guide field, where $\rho_{0}=n_{0}m_{i}$ is the background mass density, and
we have used the $k_{\perp}\rho_{i}\gg 1$ limit of GK Poisson’s law (Eq. (1)),

$\delta n_{e\lambda}/n_{0e}\sim-e\varphi_{\lambda}/(\tau T_{0e}).$ (12)

Dynamically, the parallel component of fluctuations at scale $\ell\sim
k_{\parallel}^{-1}$ is assumed to be determined by the propagation of KAWs,
and so the characteristic (parallel) time scale is set by the (linear)
frequency of those waves:

$\gamma_{l}\sim\omega_{\rm KAW}\propto
k_{\perp}\rho_{s}k_{\parallel}V_{A}\sim\rho_{s}V_{A}/(\ell\lambda),$ (13)

where $\rho_{s}=\rho_{i}/\sqrt{2\tau}$ is the ion sound Larmor radius. The
perpendicular fluctuations at scale $\lambda\sim k_{\perp}^{-1}$ undergo
nonlinear interactions which lead to a cascade of constant energy flux toward
smaller scales. Dimensionally, this inverse nonlinear time scale, usually
referred to as the “eddy-turnover rate”, is

$\gamma_{nl}\sim\varepsilon/(\rho_{0}v_{A\lambda}^{2}/2)\sim\varepsilon/(\delta
B_{\perp\lambda}^{2}/8\pi),$ (14)

where $\varepsilon$ is the constant energy flux and $v_{A\lambda}=\delta
B_{\perp\lambda}/\sqrt{4\pi\rho_{0}}$. The critical balance [Goldreich1995]
conjecture declares that these two frequencies should be comparable at each
scale, $\gamma_{l}\sim\gamma_{nl}$, setting the energy cascade rate of the
turbulence.

The above arguments are relatively general (not model-specific) as they do not
involve any assumption on the nonlinear physics of the fluctuations. In what
follows, we discuss three existing models for describing the spectrum of KAW
turbulence in the low-$\beta$ regime. Each model makes specific assumptions
about the nonlinear turbulent cascade and derives the corresponding spectra
and spectral anisotropy.

#### KAW cascade model.

Conventional models [cho2004anisotropy, howes2008model,
schekochihin2009astrophysical] of KAW turbulence assume isotropic fluctuations
in the field-perpendicular plane (i.e., $\xi\sim\lambda$) and a Kolmogorov-
like energy cascade; the nonlinear eddy-turnover rate (Eq. (14)) can thus be
written as $\gamma_{nl}\sim u_{\lambda}/\lambda$, where
$u_{\lambda}\sim(c/B_{0})\varphi_{\lambda}/\lambda$ is the velocity of the
$\mbox{\boldmath{$E$}}\times\mbox{\boldmath{$B$}}$ flow. Combined with the
equipartition relation Eq. (11), a $\propto k_{\perp}^{-7/3}$ energy spectrum
for magnetic and density fluctuations is obtained, with a
$k_{\parallel}\propto k_{\perp}^{-1/3}$ spectral anisotropy following from
imposing critical balance.

#### Intermittency model.

As discussed in the Introduction, observations and numerical simulations yield
spectra steeper than the $k_{\perp}^{-7/3}$ predicted by this simple KAW
cascade model. To address this discrepancy, boldyrev2012spectrum conjectured
that intermittent effects might play an important role in KAW turbulence. They
argued that density and magnetic fluctuations, and thus their energy cascade,
are concentrated in highly-intermittent, quasi-two-dimensional (2D) structures
elongated in the $\hat{\mbox{\boldmath{$z$}}}$-direction. Such 2D structures
can be viewed as “energy-containing sheets” (not to be confused with current
sheets, where the gradient of the magnetic fluctuation is large) 111In the
conventional $\beta$-model [frisch1995turbulence] for intermittency widely
used in hydrodynamics, the dimension of intermittent structures is a fractal
dimension and is not necessarily the same as their spatial dimension. However,
in boldyrev2012spectrum and other turbulence simulations, it is found that
intermittent structures do appear to be 2D in position space.. The volume-
filling fraction of these “active” 2D structures at each scale $\lambda$
should scale as $p_{\lambda}\propto\lambda$ [frisch1995turbulence]. It follows
that the energy density of the fluctuations should scale as
$E_{B_{\lambda}}\sim\delta B_{\perp\lambda}^{2}p_{\lambda}$, and the energy
cascade rate (Eq. (14)) thus becomes $\gamma_{nl}\sim\varepsilon/(\delta
B_{\perp\lambda}^{2}\ p_{\lambda})$. The turbulent eddies, however, are
assumed to be isotropic on field-perpendicular planes (consistent with the
numerical results reported below in Sec. V.1). The intermittent “energy-
containing” sheets can be viewed as chains of such cylindrical magnetic
structures. No reduction of nonlinearity is expected for such circular eddies,
and thus the nonlinear eddy-turnover rate should remain unchanged:
$\gamma_{\rm nl}\sim\varphi_{\lambda}/\lambda^{2}$. Combined with the
equipartition relation (Eq. (11)), this leads to the scaling $\delta
B_{\perp\lambda}\propto\lambda^{1/3}$ and energy density
$E_{B_{\lambda}}\propto\lambda^{5/3}$. The perpendicular magnetic energy
spectrum is thus

$\displaystyle E_{B}(k_{\perp})dk_{\perp}\propto k_{\perp}^{-8/3}dk_{\perp}.$
(15)

Balancing the frequency of KAWs and the nonlinear turnover rate, $\omega_{\rm
KAW}\sim\gamma_{\rm nl}$, yields $\ell\propto\lambda^{2/3}$ or, equivalently,

$\displaystyle k_{\parallel}\propto k_{\perp}^{2/3},$ (16)

and the parallel magnetic energy spectrum

$\displaystyle E_{B}(k_{\parallel})dk_{\parallel}\propto
k_{\parallel}^{-7/2}dk_{\parallel}.$ (17)

We denote by $E_{\varphi}(k_{\perp})$ the spectrum of the ion perturbed
entropy $K_{i}$ (the first term on the r.h.s. of the total free energy
invariant Eq. (6)). For $k_{\perp}\rho_{i}\gg 1$, $K_{i}$ reduces to the
energy of density fluctuations, and for $k_{\perp}\rho_{i}\ll 1$ (in RMHD),
$K_{i}$ reduces to the kinetic energy of the
$\mbox{\boldmath{$E$}}\times\mbox{\boldmath{$B$}}$ flow. In both regimes
(i.e., for all the scales above $d_{e}$ ), $K_{i}$ is expected to be in
equipartition with the magnetic energy, and so in the KAW turbulence,
$E_{\varphi}(k_{\perp})dk_{\perp}\propto k_{\perp}^{-8/3}dk_{\perp}$.

#### Tearing-mediation model.

Another fundamental process that could modify the physics of the nonlinear
cascade is the onset of tearing instability in the turbulent eddies. In
analogy with the reasoning underlying the conjecture of tearing-mediation of
MHD turbulence [loureiro2017role, mallet2017statistical,
boldyrev2017magnetohydrodynamic], if dynamical alignment
[boldyrev2006spectrum] were to be present in the sub-$\rho_{i}$ range and lead
to sufficiently elongated eddies, tearing modes could become unstable in those
eddies and modify the spectrum of the turbulence. A critical assumption of
tearing-mediated turbulence models [loureiro2017role, mallet2017disruption,
loureiro2017collisionless] is that the eddy turnover rate, or the nonlinear
frequency, is determined by the time scale of the tearing instability, i.e.,
$\gamma_{nl}\sim\gamma_{t}$, where $\gamma_{t}\sim
v_{A\lambda}\rho_{i}d_{e}^{1/2}\lambda^{-5/2}$ is the growth rate of the
fastest-growing tearing mode in a collisionless plasma for scales
$\lambda>\rho_{i}$, assuming sinusoidal (rather than Harris-sheet-like)
profiles of magnetic fluctuations [loureiro2017collisionless]. This leads to
$v_{A\lambda}\propto\lambda^{5/6}$, and the spectrum of $\delta
B_{\perp\lambda}^{2}$ thus scales as

$\displaystyle E_{B}(k_{\perp})dk_{\perp}$ $\displaystyle\propto
k_{\perp}^{-8/3}dk_{\perp}.$ (18)

Invoking, as before, the critical balance conjecture leads to the spectral
anisotropy $k_{\parallel}\propto k_{\perp}^{2/3}$. Note the fact that these
are exactly the same scalings with $k_{\perp}$ as found in the intermittency
model of boldyrev2012spectrum, even though the physics underlying each of
these models are of a completely different nature.

The derivation above assumes that the magnetic field profile in the eddies is
sinusoidal (i.e., $\delta B_{\lambda}\sim\sin(x/\lambda$)), and that the
frozen flux constraint is broken by electron inertia. Different assumptions of
spatial profiles of the fluctuations, or different choices of flux unfreezing
mechanisms, lead to different scalings of a tearing-mediated spectrum, whereas
the intermittency is not expected to be affected. If a Harris sheet-type
configuration (i.e., $\tanh(x/\lambda)$ profile) of magnetic fluctuations is
assumed instead, the derived spectrum and anisotropy become
$E_{B}(k_{\perp})\propto k_{\perp}^{-3}$ and $k_{\parallel}\propto k_{\perp}$
222For the Harris sheet-type configuration, the growth rate of the fastest
tearing mode becomes $\gamma_{t}\sim v_{A\lambda}d_{e}\rho_{s}/\lambda^{3}$
[loureiro2017collisionless], leading to the scale dependence of magnetic
fluctuation $v_{A\lambda}\sim\varepsilon^{1/3}\lambda(d_{e}\rho_{s})^{-1/3}$,
the magnetic spectrum $E_{B}(k_{\perp})\sim
4\pi\rho_{0}\varepsilon^{2/3}(d_{e}\rho_{s})^{-2/3}k_{\perp}^{-3}dk_{\perp}$,
and the spectral anisotropy
$k_{\parallel}\sim\varepsilon^{1/3}(d_{e}^{2}/\rho_{s})^{1/3}V_{A}^{-1}k_{\perp}$..
In numerical simulations, electron inertia can be replaced by (hyper-)
resistivity as the flux-unfreezing mechanism by taking $d_{e}\rightarrow 0$
(more precisely, by setting the (hyper) resistivity to be such that the
corresponding thickness of the inner layer of the tearing mode is larger than
$d_{e}$). In this case, if the cascade is indeed tearing-mediated, the
expected magnetic spectrum would be

$E_{B}(k_{\perp})dk_{\perp}\propto
k_{\perp}^{-(7n\alpha+2\alpha+2)/(3n\alpha+2)}dk_{\perp},$ (19)

where $\alpha$ is the order of the spatial derivative in the (hyper) resistive
term and $n=1$ ($n=2$) corresponds to the Harris-sheet (sinusoidal) profile [a
detailed derivation of this expression is provided in Sec. .1 of the
Supplemental Materials (SM)]. The flexibility introduced by the use of
different flux-unfreezing methods enables the identification of which of these
turbulence models — intermittency or tearing-mediation — is a better
description of the turbulent cascade at sub-ion scales. This is one of the
main goals of our numerical simulations (see Sec. V.1).

###  III.2 Phase-space cascade and electron heating

The above-described models ignore electron kinetic physics. In reality,
however, electrons are not isothermal and the kinetic physics of electrons
leads to complex dynamics in velocity space, relevant to energy dissipation
and electron heating. It follows that velocity-space dynamics is as important
to the understanding of KAW turbulence as real-space dynamics — and, indeed,
we will show that the two are tightly connected.

In velocity space, the phase-mixing rate can be estimated through terms on the
r.h.s. of Eq. (10) which represent the transfer of free energy from the $m$-th
to $(m+1)$-th moment. In the large-$m$ limit, these terms can be written
approximately as the derivative with respect to $m$, from which follows that
the phase-mixing rate can be approximated as $\sim|k_{\parallel}|v_{{\rm
th}e}/\sqrt{m}\sim(v_{{\rm th}e}/B_{0})\\{A_{z},...\\}/\sqrt{m}$
[zocco2011reduced, schekochihin2016phase]. In position space, on the other
hand, free energy cascades to smaller spatial scales at the nonlinear
advection rate $d/dt\sim(c/B_{0})\\{\varphi,...\\}$ (l.h.s. of Eq. (10)).
These rates are not, in general, the same. At each scale $\lambda$, there is a
critical Hermite mode number, denoted as $m_{\rm cr}$, at which these two
rates balance:

$\frac{c}{B_{0}}\\{\varphi,g_{m_{\rm cr}}\\}\sim\frac{v_{{\rm
th}e}}{B_{0}}\\{A_{z},g_{m_{\rm cr}}\\}/\sqrt{m_{\rm cr}},$ (20)

resulting in the relation $\sqrt{m_{\rm cr}}\sim(v_{{\rm
th}e}/c)A_{z\lambda}/\varphi_{\lambda}$. Combined with the GK Poisson’s law
(Eq. (1)) and the equipartition relation (Eq. (11)), we obtain the spatial
scale dependence of the critical Hermite order:

$m_{\rm cr}(\lambda)\sim(\lambda/d_{e})^{2}/(2\tau^{2}).$ (21)

That is, for $\sqrt{m}k_{\perp}d_{e}\gg 1$, the nonlinear advection is
expected to dominate, and the rapid coupling between modes with different
$k_{\parallel}$ can result in a new mode with negative $k_{\parallel}$ and
with a returning free-energy flux to lower Hermite moments. This classic
phenomenon is know as the “plasma echo” [Gould1967echo, malmberg1968plasma,
schekochihin2016phase].

In many astrophysical environments, collisions are so infrequent that the
collisional cutoff in velocity space will only occur at asymptotically large
$m$ ($\gg m_{\rm cr}$). Therefore, one may expect the free energy to reach
$m_{\rm cr}$ and return to the low Hermite moments due to the restoring
stochastic plasma echo. Kinetic turbulence would thus resemble fluid
turbulence, in the sense that the only available route for energy dissipation
would be the real-space cascade. However, this argument is seemingly at odds
with measurements of strong electron heating in the near-Earth solar wind
[chen2019evidence] and in kinetic simulations [howes2018spatially,
klein2016measuring, howes2017diagnosing]. Consistently, in the numerical
results reported in Sec. V, we observe no evidence for plasma echoes. We argue
that this paradox is resolved by the presence of spontaneously formed current
sheets in kinetic turbulent systems. On the one hand, efficient phase mixing
is expected (as a result of the local suppression of nonlinearity) around
current sheets [loureiro2013fast, numata2015]. On the other hand, it is at the
sites with large gradient of current density (around the edge of current
sheets) that most of the energy in fluid quantities is pumped into Hermite
moments, because the coupling between them is through
$\propto\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}(e/cm_{e})d_{e}^{2}\nabla_{\perp}^{2}A_{z}\propto\hat{\mbox{\boldmath{$b$}}}\cdot\mbox{\boldmath{$\nabla$}}J_{z}$
[last term in Eqs. (4) and (10)]. The combination of these two effects of
current sheets locally regulates the kinetic turbulence so that the free
energy can freely transfer to large Hermite moments, unimpeded by echo, and
collisionally dissipate, leading to electron heating. Therefore, the phase-
mixing-dominated regime in phase-space can be much wider than the estimation
(Eq. (21)) based on the simple time-scale comparison (Eq. (20)).

###  III.3 A lowest-order solution for the Hermite-expansion coefficients of
$g_{e}$

Based on the time-scale comparison discussed in the previous section, we
present a lowest-order solution for the perturbed electron distribution
function $g_{e}$ in velocity space, valid in the phase-mixing-dominated
regime. In Eq. (10), each term can be ordered based on its inherent frequency.
In the phase-mixing-dominated regime, as discussed in Sec. III.2, the l.h.s.
(proportional to the nonlinear-advection rate) of Eq. (10) is much smaller
than its r.h.s. (proportional to the phase-mixing rate). To lowest order, the
two terms on the r.h.s. should therefore balance, and each Hermite moment
should satisfy the relation

$\displaystyle g_{m+1}=-\sqrt{m/(m+1)}g_{m-1}$ $\displaystyle\text{for}\ m\geq
3,\ \text{and}$ (22) $\displaystyle g_{3}=-2/\sqrt{3}\ e/(cm_{e}v_{{\rm
th}e})d_{e}^{2}\nabla_{\perp}^{2}A_{z}.$

That is, to lowest order, all the odd Hermite moments share the same spatial
configuration as the current, and all the even ones that of the temperature
fluctuations, $g_{2}$.

The Hermite spectrum of $g_{e}$, $E_{m}\equiv\braket{}{g_{m}}{{}^{2}}/2$ where
$\braket{...}$ represents the volume average, can be derived from Eq.
(LABEL:eq:gm_solution). With the approximation
$E_{m}\approx\braket{}{g_{m}}{g_{m-1}|}/2$ at large $m$, and using the
relation $\lvert g_{m+1}/g_{m-1}\rvert=\sqrt{m/(m+1)}$, we obtain
$E_{m+1}/E_{m-1}\approx g_{m+1}g_{m}/(g_{m-1}g_{m-2})=\sqrt{(m-1)/(m+1)}$.
This recursive relation leads to the solution $E_{m}\propto m^{-1/2}$ — the
same as the expected spectrum in a linear system with a Landau-damped kinetic
field, corresponding to a constant flux of free energy transferring from small
to large velocity moments [zocco2011reduced, kanekar2015fluctuation].

###  III.4 Turbulence at scales below the electron skin depth, $\lambda\ll
d_{e}$

Lastly, one may also use the KREHM formalism to investigate the turbulent
cascade at scales below $d_{e}$, i.e., in the range $\lambda\ll
d_{e}\lesssim\rho_{i}$, where fluctuations become predominantly electrostatic.
In this range, equipartition between the density fluctuations and kinetic
energy of the parallel electron flows, $(\delta
n_{e}/n_{0e})^{2}n_{0e}T_{0e}\sim d_{e}^{2}|\nabla_{\perp}^{2}A_{z}|^{2}/8\pi$
is expected, leading to the relation
$\varphi_{\lambda}\sim(\rho_{i}V_{A}/c)d_{e}\delta B_{\perp\lambda}/\lambda$
(invoking Eq. (12)). Since electromagnetic effects are subdominant in this
scale range, neither the intermittency corrections nor the tearing mediation
that we discussed earlier are expected to be important here. One might thus
expect a standard Kolmogorov-type cascade model to be a reasonable description
of turbulence at these scales. The energy flux of the cascade is
$\varepsilon\sim\gamma_{\rm nl}e^{2}n_{0e}\varphi^{2}/T_{0e}$, where the eddy
turn-over rate is $\gamma_{\rm nl}\sim\varphi_{\lambda}/\lambda^{2}$. Together
with the equipartition relation, we obtain the kinetic and magnetic spectra

$\displaystyle E_{\varphi}(k_{\perp})\propto k_{\perp}^{-7/3},\quad
E_{B}(k_{\perp})\propto k_{\perp}^{-13/3}.$ (23)

Imposing critical balance of these fluctuations, that is, $\omega_{\rm
KAW}\sim\gamma_{nl}$ (at $k_{\perp}d_{e}\gg 1$, $\omega_{\rm KAW}\sim
k_{\parallel}V_{A}\rho_{s}/d_{e}$), we obtain $k_{\parallel}\propto
k_{\perp}^{4/3}$.

As before, let us now estimate the critical Hermite moment, $m_{\rm cr}$,
where the rate of free energy transfer between consecutive Hermite moments
through linear phase mixing balances the nonlinear advection rate, following
Eq. (20). Combined with the equipartition relation, we obtain $m_{\rm
cr}\sim(\lambda/d_{e})^{4}/(2\tau^{2})$. That is, at scales $k_{\perp}d_{e}\gg
1$, the nonlinear advection of free energy is always faster than its
channeling to larger Hermite moments. Linear phase mixing should therefore be
subdominant as an energy dissipation channel, and the isothermal closure
should be a good approximation to the dynamics at these scales.

##  IV Numerical setup

We perform direct numerical simulations of the KREHM equations (see Section
II) using the pseudo-spectral code Viriato [loureiro2016viriato]. In what
follows, quantities are given in dimensionless form. The domain is a triply
periodic cubic box with sides of length $L=2\pi$. Hyper-diffusion terms of the
form $\nu_{H}\mbox{\boldmath{$\nabla$}}_{\perp}^{6}$ are included in the
r.h.s. of Eq. (2) (electron hyper-viscosity), Eq. (3) (hyper-resistivity), and
Eq. (10), with $\nu_{H}$ dynamically adjusted to absorb energy at the grid
scale [loureiro2016viriato]. A hyper-collision operator of the form $-\nu_{\rm
coll}m^{6}$ is added to the r.h.s. of Eq. (10) for $m\geq 2$, with $\nu_{\rm
coll}$ likewise set to remove energy at the smallest velocity scale,
determined by $M$, the order of the highest Hermite polynomial kept in
simulations. The simulations are driven with a white-noise forcing term added
to Eq. (2), which injects energy into the largest scales in the simulation box
(perpendicular and parallel wavenumbers ranging from 1 to 2). Once a steady
state is reached, the energy injection is balanced by its dissipation either
through hyper-collisions at large $m$ or through hyper-diffusion at large
$k_{\perp}$, self-consistently determined by the nature of kinetic turbulence.
As mentioned in Section I, it is one of the goals of this study to identify
which of these dissipation channels is privileged by the turbulence.

We perform one set of simulations with isothermal electrons (i.e., $g_{e}=0$),
and another set containing electron kinetic effects ($g_{e}\neq 0$, meaning
that Eq. (10) is solved for the Hermite moments of electron perturbed
distribution). The key simulation parameters are summarized in Table 1 and
Table 2, respectively. For the isothermal simulations, we focus on the sub-ion
range ($\rho_{i}>k_{\perp}^{-1}>d_{e}$) and vary the dominant mechanisms that
unfreeze the magnetic flux (such as to be able to distinguish intermittency
from tearing in regards to the physical mechanism that determines the energy
spectrum; see discussion in Section III). For the non-isothermal (i.e.,
kinetic) simulations, we investigate three dynamical ranges of interest: (1)
the inertial range $L>k_{\perp}^{-1}>\rho_{i}$, (2) the sub-ion range
$\rho_{i}>k_{\perp}^{-1}>d_{e}$, and (3) the electrostatic range
$d_{e}>k_{\perp}^{-1}$. In run K1, we include both the inertial and
sub-$\rho_{i}$ ranges, enabling us to study the transition across the
$\rho_{i}$ scale. However, the resolution of the $d_{e}$ scale is sacrificed
and the effect of electron inertia on turbulence is not properly taken into
account. This is compensated by runs K2(a,b) where the whole domain is in the
sub-ion range and $d_{e}$ is well resolved. Run K3 is designed to focus on the
sub-$d_{e}$ range only. Runs K4(a-e) consist of a scan of the parameter
$\rho_{i}/d_{e}$ with fixed $d_{e}$ and varying $\rho_{i}$ and are used to
study the dependence of the energy dissipation on the scale separation between
ion and electron scales.

Runs |  $\rho_{i}/L$ |  $d_{e}/L$ |  $\eta$ |  $\nu_{H}$ |  $N^{3}$ |  flux-unfreezing
---|---|---|---|---|---|---
Iso-1 |  2 |  0.03 |  0 |  0.11e-9 |  $512^{3}$ |  electron inertia
Iso-2 |  2 |  0.001 |  0.01 |  0.39e-9 |  $512^{3}$ |  resistivity
Iso-3 |  2 |  0.001 |  0 |  0.39e-9 |  $512^{3}$ |  hyper-resistivity
Table 1: Summary of isothermal ($g_{e}=0$) simulations’ parameters. The last column shows the dominant mechanism that breaks the frozen-flux constraint in the simulation. Runs |  $\rho_{i}/L$ |  $\rho_{i}/d_{e}$ |  $k_{\perp{\rm max}}d_{e}$ |  $M$ |  $N^{3}$
---|---|---|---|---|---
K1 |  0.1 |  20 |  0.9 |  30 |  $256^{3}$
K2a |  2 |  40 |  9 |  30 |  $256^{3}$
K2b |  2 |  20 |  9 |  60 |  $128^{3}$
K2c |  2 |  66 |  5.4 |  30 |  $256^{3}$
K3 |  2 |  1 |  180 |  30 |  $128^{3}$
K4a-e |  3,2,1,0.3,0.1 |  30,20,10,3,1 |  9 |  30 |  $128^{3}$
Table 2: Summary of kinetic ($g_{e}\neq 0$) simulations’ parameters.
$k_{\perp{\rm max}}$ is the maximum value of the perpendicular wavenumber,
representing the resolution of simulations.

##  V Numerical results

In Sec. V.1, we present the results from runs Iso-(1-3) (with isothermal
electrons) on the spatial energy spectra, structure functions of the
fluctuations, and intermittency. With such understanding of the isothermal
limit of sub-$\rho_{i}$ turbulence, in Sec. V.2 we proceed to take into
account the electron kinetic effects and study the spectra in position and
velocity space, the energy dissipation and electron heating, and the
correlation of Hermite moments.

###  V.1 Energy spectra set by intermittency

Figure 1: Comparison of magnetic spectra from isothermal runs (see Table 1)
with different flux-freezing-breaking mechanisms (top). Structure functions
(middle) and spectral anisotropy (bottom) from the run Iso-1 (the other Iso
runs show the similar results).

In Sec. III.1 we discussed three models for energy spectra based on distinct
physics that might govern the nonlinear cascade. As we noted, in spite of the
different underlying physical mechanisms, the spectra and spectral anisotropy
predicted by the intermittency model and the tearing-mediation model (with
electron inertia effects) are the same: $E_{B}(k_{\perp})\propto
k_{\perp}^{-8/3}$, $E_{B}(k_{\parallel})\propto k_{\parallel}^{-7/2}$, and
$k_{\parallel}\propto k_{\perp}^{2/3}$. This prediction is confirmed by run
Iso-1 in which $d_{e}$ is resolved (and thus with electron inertia effects
accounted for). The normalized (to its maximum value) magnetic spectrum
exhibits a $\propto k_{\perp}^{-8/3}$ scaling, as shown in Fig. 1, top panel.

To obtain a more detailed description of the turbulent fluctuations in these
simulations we compute the three-dimensional structure functions of the
magnetic fluctuations $S_{q}(\delta r)\equiv\braket{}{\Delta
B(\mbox{\boldmath{$r$}},\delta\mbox{\boldmath{$r$}})}{{}^{q}}_{\mbox{\boldmath{$r$}}}$,
where $\delta r$ is the incremental scale, $\Delta
B(\mbox{\boldmath{$r$}},\delta\mbox{\boldmath{$r$}})$ is the local-field
increment calculated with a five-point stencil, and $q$ is the order of the
structure function. The detailed description of this diagnostic and of its
numerical implementation can be found in Sec. 3 in SM. We use the second-order
structure function $S_{2}(\delta r)$ to study the parallel spectrum and
spectral anisotropy by projecting it into local coordinates with respect to
the local magnetic field. We denote by $S_{2}(\lambda)$ and $S_{2}(\ell)$ the
field-perpendicular and field-parallel structure functions, which are mapped
to the perpendicular and parallel magnetic energy spectrum through Parseval’s
theorem: a $\sim k^{-\alpha}$ energy spectrum corresponds to the structure
function $S_{2}(\delta r)\propto\delta r^{\alpha-1}$ with $\delta r\sim
k^{-1}$. Therefore, the $E_{B}(k_{\perp})\propto k_{\perp}^{-8/3}$ (Eq. (15))
and $E_{B}(k_{\parallel})\propto k_{\parallel}^{-7/2}$ (Eq. (17)) spectra
should correspond to the scalings $S_{2}(\lambda)\propto\lambda^{5/3}$ and
$S_{2}(\ell)\propto\ell^{5/2}$, respectively. These predicted scalings are
consistent with the measurements shown in the middle panel of Fig. 1. The
predicted spectral anisotropy based on the critical balance of the turbulent
fluctuations, $k_{\parallel}\propto k_{\perp}^{2/3}$, is also confirmed by the
measured $\ell\propto\lambda^{2/3}$ scaling shown in the bottom panel.
Finally, we find that the fluctuations exhibit a $\xi\propto\lambda$ scaling
(where $\xi$ is the coherence length in the fluctuation direction, as
explained in section III.1) on perpendicular (to the local magnetic field)
planes. This scaling indicates the isotropic morphology of eddies in
perpendicular planes. That is, the electromagnetic fluctuations in our
simulations do not become progressively more aligned with decreasing scale (in
contrast with what is conjectured to happen in MHD turbulence
[boldyrev2006spectrum, mallet2016, schekochihin2020mhd]), which is the
necessary premise for the hypothesis of tearing-mediation of the energy
cascade.

It is worth noting that the three-dimensional eddy anisotropies at sub-ion
scales have recently been investigated through in situ measurements by the MMS
in the Earth’s magnetosheath [wang2020observational] and by the PSP in the
inner heliosphere [zhang2022three]; the scalings between $\ell$, $\xi$, and
$\lambda$ that these investigations report are quantitatively similar to our
measurements shown in the bottom panel of Fig. 1.

As explained earlier, the above results for the magnetic spectrum, by
themselves, do not clarify whether it is intermittency or tearing-mediation
that determines the spectrum. To shed light on this issue, we perform two
additional runs whose only difference from the one we just discussed (Iso-1)
is the flux-unfreezing mechanism that is used: instead of electron inertia, it
is Laplacian resistivity for run Iso-2, and hyper-resistivity for run Iso-3.
In the top panel of Fig. 1, we compare the perpendicular magnetic spectra from
all three isothermal runs. If the spectra were to be set by tearing-mediation,
we would expect the scalings $E_{B}(k_{\perp})\propto k_{\perp}^{-2.4}$ and
$E_{B}(k_{\perp})\propto k_{\perp}^{-2.6}$ for runs Iso-2 and Iso-3,
respectively, according to Eq. (19) for the $n=2$ sinusoidal profile. Instead,
we observe a clear overlap of the spectra for all three simulations (in the
inertial range). This finding, combined with the measured $\xi\sim\lambda$
scaling (Fig. 1, bottom panel), rules out tearing-mediation as the physical
mechanism underpinning the KAW turbulent cascade in our simulations, and
leaves intermittency as the only known explanation for the $k_{\perp}^{-8/3}$
magnetic spectrum in the low-$\beta$, sub-$\rho_{i}$ regime.

Figure 2: The exponents of $q$-th-order structure functions in the
perpendicular plane for isothermal runs. The deviation of high-order exponents
from the linear scaling $\zeta_{1}q$ (black dashed line) indicates the
presence of intermittency.

To characterize the intermittency of the turbulent fluctuations in our
simulations we investigate the higher order ($q>2$) structure functions. In
Fig. 2, we present the dependence on the order $q$ of the measured index
$\zeta_{q}$ for the field-perpendicular structure function
$S_{q}(\lambda)\propto\lambda^{\zeta_{q}}$ from the three isothermal runs. In
a turbulent system with purely self-similar fluctuations (i.e., no
intermittency), a linear scaling $\zeta_{q}=\zeta_{1}q$ is expected
[frisch1995turbulence]. The deviation from this linear scaling as $q$
increases confirms the presence of intermittency in our system. Furthermore,
the measured scaling for $\zeta_{q}$ is consistent with the $\beta$-model
[frisch1995turbulence] and the main assumption made in boldyrev2012spectrum
that the intermittent structures are essentially 2D. There are two
contributions to $S_{q}(\lambda)$: the fluctuations $\delta B_{\perp\lambda}$
from the “active” intermittent regions, and their volume-filling fraction
$p_{\lambda}$, leading to the relation $S_{q}(\lambda)\equiv\braket{}{\Delta
B(\mbox{\boldmath{$r$}},\lambda)}{{}^{q}}_{\mbox{\boldmath{$r$}}}\sim\delta
B_{\perp\lambda}^{q}p_{\lambda}\sim\lambda^{q/3+1}$. This $\zeta_{q}=q/3+1$
scaling for intermittency is indeed measured in the our simulations (with the
small deviation of run Iso-1 caused by the relatively narrow inertial range
compared to the other runs.).

###  V.2 Electron heating and the absence of plasma echo

Figure 3: Runs K1 (left column), K2a (middle column), and K3 (right column).
Top: visualization of $J_{z}$ (color map) and contours of the electron heating
rate that is twice its volume-averaged value (magenta lines). Middle: Magnetic
energy and density spectra. Bottom: Energy dissipation spectra in
$m$-$k_{\perp}\rho_{i}$ phase space. The vertical dotted lines indicate the
$d_{e}$ scale.

With the understanding of the isothermal limit in place, we proceed to
investigate the more realistic case where electron kinetic physics — and, in
particular, electron Landau damping — is included in the simulations (i.e.,
$g_{e}\neq 0$). As discussed above, this enables the cascade of energy in
velocity (or Hermite), as well as in real (or Fourier), space.

In Fig. 3 we present results from the three main non-isothermal runs (K1, K2a,
and K3, corresponding to the left, middle, and right panels) that altogether
cover the spatial scales ranging from MHD to sub-$d_{e}$. Run K1 captures the
transition at $\rho_{i}$ and run K2 captures that at $d_{e}$. For each run,
the results are plotted at an arbitrarily chosen time during the saturated
stage of the turbulence.

#### Visuals.

The top row shows the current densities normalized to their root-mean-square
values on arbitrarily chosen $xy$ planes. The visualization of run K1 (left
panel) is dominated by the features of long-wavelength MHD turbulence, with
indications of plasmoid instability [loureiro2007instability] in some of the
large-scale current sheets. At scales below $\rho_{i}$ (middle panel), the
current still exhibits sheet-type structures, but there are no visual signs of
plasmoid generation. At scales below $d_{e}$ (right panel), the current
profile becomes less sharp and fewer finer structures appear in the system,
consistent with the predicted steep magnetic spectrum (Eq. (23)) and with the
fact that turbulence at these scales should be predominantly electrostatic.

#### Energy spectra.

The middle row shows the normalized perpendicular energy spectra of the
magnetic and ion-entropy fluctuations, $\hat{E}_{B}(k_{\perp})$ and
$\hat{E}_{\varphi}(k_{\perp})$. The main purpose of run K1 (left panel) is to
study the transition at around $k_{\perp}\rho_{i}\sim 1$, and so both the MHD
range $k_{\perp}\rho_{i}<1$ and the kinetic range $k_{\perp}\rho_{i}>1$ are
included in the simulation, with the compromise that neither of these two
ranges are sufficiently wide to form a clear inertial range. In spite of this
limitation, there is a clear steepening of the spectra at
$k_{\perp}\rho_{i}\sim 1$. We observe
$\hat{E}_{B}(k_{\perp})\approx\hat{E}_{\varphi}(k_{\perp})$ in the
sub-$\rho_{i}$ range, confirming the assumption of equipartition between
magnetic and density fluctuations (Eq. (11)).

The sub-$\rho_{i}$ range is explored in run K2a (middle panel) where the
transition at $k_{\perp}d_{e}\sim 1$ is well resolved. Both
$\hat{E}_{B}(k_{\perp})$ and $\hat{E}_{\varphi}(k_{\perp})$ closely follow a
$k_{\perp}^{-3}$ scaling until they deviate at scales below $d_{e}$. This is a
steeper spectrum than the predicted (and numerically confirmed)
$k_{\perp}^{-8/3}$ spectrum (Eq. (15)) in the isothermal limit. We attribute
this steepening of the spectra to electron Landau damping, which enables the
dissipation of a non-negligible energy fraction via its cascade to high order
Hermite moments through phase-mixing. This argument will be elaborated and
tested below. At scales below $d_{e}$ (run K3, right panel),
$\hat{E}_{B}(k_{\perp})\propto k_{\perp}^{-13/3}$ and
$\hat{E}_{\varphi}(k_{\perp})\propto k_{\perp}^{-7/3}$ are measured,
confirming the predictions for the electrostatic limit (Eq. (23)).

#### Energy dissipation.

The bottom row shows contour maps of the normalized (to its maximum value)
dissipation rate of the total free energy (defined in Eq. (6)) in the
$m$-$k_{\perp}$ phase space. Here the $m=0$ and $m=1$ moments correspond to
$\delta n_{e}/n_{0e}$ and $A_{z}$, respectively. These results directly
address one of the main questions in this paper — the competition between the
kinetic and fluid channels for energy dissipation. In run K1 (left panel), we
focus on the MHD range and the transition into sub-$\rho_{i}$ scales. The
strong energy dissipation at large ($m\gg 1$) Hermite moments turns on at
around the $\rho_{i}$ scale and extends towards smaller scales (larger
$k_{\perp}$), suggesting the domination of phase mixing over nonlinear
advection at those scales 333The large-$m$ dissipation at the system scale is
artificial, caused by the forcing.. Because $d_{e}$ is unresolved in this run,
the energy flux in position space reaches the (hyper-)viscous scale without
going through the full kinetic range. This truncation of electron kinetic
effects at large $k_{\perp}$ causes an artificially strong dissipation by the
hyper-viscosity (and resistivity), preventing us from making a direct
comparison of dissipation at large $m$ and large $k_{\perp}$ in this run. This
issue is resolved in run K2a (middle panel). With the understanding that the
large-$m$ dissipation turns on at $k_{\perp}\rho_{i}\sim 1$, we can study the
sub-$\rho_{i}$ range scale and focus on the range between $\rho_{i}$ and
$d_{e}$, with $d_{e}$ properly resolved. With electron kinetic physics fully
accounted for, we find that most of the energy is dissipated at large $m$
(rather than large $k_{\perp}$). This result is direct evidence of the absence
of plasma echo in the range between $\rho_{i}$ and $d_{e}$, and clearly
demonstrates that the kinetic channel is dominant for energy dissipation: the
free energy transfers to small scales in velocity space through phase mixing,
where it dissipates via collisions and heats the electrons.

The large-$m$ dissipation starting at the lowest $k_{\perp}$ turns off at
around the $d_{e}$ scale. To explore the sub-$d_{e}$ range, we analyze run K3
(right panel), and find that the energy dissipation at large $k_{\perp}$ is
clearly dominant. This is consistent with our prediction (in Sec. III.4) that
although the Landau damping rate of linear KAWs is large in the sub-$d_{e}$
range, the nonlinear advection is faster than the phase-mixing process,
causing the energy flux to mainly cascade to, and dissipate at, large
$k_{\perp}$ rather than large $m$.

Figure 4: Comparison between the magnetic spectrum from run K2c
($\hat{E}_{B}^{\rm K}$) and that from run Iso-1 ($\hat{E}_{B}^{\rm Iso}$)
multiplied by a scale-dependent factor $(1-2\gamma_{\rm KAW}/\omega_{\rm
KAW})$ that represents the subtraction of the Landau-damped energy within a
cascade time.

#### Landau damping of KAWs.

In order to confirm that the strong energy dissipation at high $m$ is indeed
caused by the Landau damping of KAWs, we test the following argument with our
numerical results. Let us assume that, at each scale $k_{\perp}^{-1}$, the
electromagnetic fluctuations are damped at the KAW linear damping rate
pertaining to that scale, $\gamma_{\rm KAW}(k_{\perp})$. Then, within a
nonlinear cascade time $1/\gamma_{nl}$ — which, by assumption of critical
balance, is comparable to the linear time scale, $1/\gamma_{nl}\sim
1/\omega_{\rm KAW}$ — a fraction of $\sim 2\gamma_{\rm KAW}/\omega_{\rm KAW}$
of the magnetic energy is damped and becomes the free energy in $g_{e}$. All
the mentioned quantities are scale-dependent. With this assumption, if we
consider a magnetic spectrum from a simulation with isothermal electrons,
$\hat{E}_{B}^{\rm iso}(k_{\perp})$, and subtract from it the energy that would
be Landau damped within a cascade time, it should match the spectrum
$\hat{E}_{B}^{\rm K}$ from the corresponding kinetic simulation (i.e., a
simulation that differs from the isothermal one only in allowing for kinetic
electrons, $g_{e}\neq 0$):

$\hat{E}_{B}^{\rm K}(k_{\perp})\approx\hat{E}_{B}^{\rm
iso}(k_{\perp})(1-2\gamma_{\rm KAW}/\omega_{\rm KAW}).$ (24)

This conjecture is confirmed in the comparison of spectra from run K2a and its
corresponding isothermal run (Iso-1), shown in Fig. 4, where the scale-
dependent $\gamma_{\rm KAW}$ and $\omega_{\rm KAW}$ are calculated from the
linear dispersion relation of KAWs Eq. (S19) (see Sec. .2 in SM for more
details). This agreement suggests that the magnetic spectra from kinetic runs,
shown in this figure and in Fig. 3, which we tentatively fitted with a
$k_{\perp}^{-3}$ scaling, is, in fact, not a simple power law, but rather a
more complicated function well approximated by Eq. (24). Incidentally, we note
that it should not be an artifact of numerical simulations that the
sub-$\rho_{i}$ spectra do not exhibit power-law scalings, since the scale
separation between $\rho_{i}$ and $d_{e}$ in our simulation is not far-off
from that in realistic systems such as the solar wind. This result is a direct
numerical confirmation of the weakened cascade model originally proposed by
Howes et al. [howes2008kinetic, howes2011weakened] and numerically tested by
TenBarge and co-workers [tenbarge2013collisionless, tenbarge2013current]. In a
broader context, this result provides a justification for a model
[podesta2010kinetic, howes2010prescription, passot2015model,
kunz2018astrophysical] that is widely used in studies of turbulent
astrophysical systems, which is based on the assumption that we just
confirmed: that the effect of phase mixing is merely to damp the turbulent
cascade of the electromagnetic fluctuations at a scale-dependent rate set by
linear electron Landau damping.

Figure 5: Normalized Hermite spectra from simulations K1-K3 (top) and K4a-e
(bottom).

#### An analytical solution for the Hermite-expansion coefficients of the
reduced electron distribution function $g_{e}$.

Strong energy dissipation at high Hermite mode numbers $m\gg 1$ requires a
significant energy cascade towards small scales in the velocity space and
should thus correspond to a shallow Hermite spectrum. This is consistent with
our numerical results shown in the top panel of Fig. 5. For runs K1 and K2a
with strong large-$m$ dissipation, we measure a shallow $m^{-1/2}$ spectrum.
In contrast, for run K3 where large-$m$ dissipation is weak and
large-$k_{\perp}$ dissipation dominates, a much steeper Hermite spectrum is
measured. The $m^{-1/2}$ spectrum in run K1 is consistent with our earlier
explanation that although the dissipation at large-$k_{\perp}$ is significant
in this case, it is due to the truncation of electron kinetic physics at the
hyper-viscous scale, and a significant phase-mixing-dominated energy cascade
in velocity space does exist at spatial scales larger than the hyper-viscous
scale.

In the bottom panel of Fig. 5 we show the Hermite spectra from runs with
different values of $\rho_{i}/d_{e}$. In these runs, $d_{e}$ is fixed at small
scales (but resolved). A large value of $\rho_{i}/d_{e}$ corresponds to a
wider electromagnetic sub-$\rho_{i}$ range where KAWs exist and are subject to
Landau damping, whereas runs with small $\rho_{i}/d_{e}$ are mostly in the MHD
range. As shown, decreasing values of $\rho_{i}/d_{e}$, i.e., a narrower
dynamical range for Landau-damped KAWs results in a progressively steeper
Hermite spectrum. This result further confirms the significant contribution of
the Landau damping of KAWs to velocity-space dissipation.

The $m^{-1/2}$ scaling is predicted by our analytical model (in Sec. III.3)
based on the timescale ordering between the KAW frequency (which is also the
cascade rate) and the phase mixing rate: $\omega_{\rm KAW}\ll
k_{\parallel}v_{{\rm th}e}/\sqrt{m}$; this leads to a lowest-order solution
relating different Hermite moments, Eq. (LABEL:eq:gm_solution). To further
test our theory, we measure the correlation between Hermite moments $g_{m-1}$
and $g_{m+1}$, a few examples of which are shown in Fig. 6, top and middle
panels. The measured correlations between $J_{z}$ and $g_{3}$, $g_{5}$ and
$g_{3}$, $g_{4}$ and $g_{2}$, and between $g_{6}$ and $g_{4}$ exhibit
remarkable agreement with our lowest-order solution
$g_{3}=2/\sqrt{3}(d_{e}/\rho_{s})J_{z}$ (in code units) and
$g_{m+1}=-\sqrt{m/(m+1)}g_{m-1}$ for $m\geq 3$ (Eq. (LABEL:eq:gm_solution)).

Figure 6: Scatter plots from simulation K1 showing the correlation between
$J_{z}$ and $g_{3}$ (top left), $g_{5}$ and $g_{3}$ (top right), $g_{4}$ and
$g_{2}$ (bottom left), and between $g_{6}$ and $g_{4}$ (bottom right).

#### Absence of plasma echoes.

The unimpeded Hermite cascade of energy that we observe implies that the
stochastic echo effect that we discussed earlier is negligible in our
simulations. It is revealing to understand why that should be so. The time-
scaling ordering leads to a critical Hermite mode number $m_{\rm
cr}\sim(\lambda/d_{e})^{2}/(2\tau^{2})$ (Eq. (21)), below which the phase
mixing is faster than the nonlinear advection, and above which the advection
can in principle dominate the phase mixing to establish the plasma echo. That
seems to indicate that, for a system with asymptotically large cutoff $M$, the
plasma echo should always occur, causing a returning energy flux to low
Hermite moments and thus strong energy dissipation at large $k_{\perp}$. Even
for a system with limited $M$ (as that in our simulations), at sufficiently
small scale $\lambda$, $m_{\rm cr}$ will be resolved ($m_{\rm cr}<M$), and
phase mixing should be suppressed at those scales. However, these arguments
are contradictory to our numerical results of the $m^{-1/2}$ Hermite spectrum
and efficient electron heating.

Figure 7: Comparison between normalized total energy dissipation spectra
(integrated over $m$) between run K2b (with $M=60$; blue solid curve) and run
K4b (with $M=30$; blue dashed line). Also in the figure is the dependence of
the critical Hermite moment order, $m_{\rm cr}$ (Eq. (21)), as a function of
$k_{\perp}\rho_{i}$. The horizontal lines indicate the cutoff of $m$ at 30 and
60. For both $M=30$ and $M=60$ runs, a wide range of wavenumbers exists that
in principle allows plasma echo to occur.

To further confirm the absence of echo, we perform two additional runs, K2b
and K4b, between which the only difference is the number of Hermite moments
($M=30$ for K4b and $M=60$ for K2b). The Hermite spectrum of K2b also exhibits
a $m^{-1/2}$ scaling (Fig. 5, top panel). In Fig. 7 we compare the dissipation
rate of total free energy between the K2b ($M=60$) and K4b ($M=30$) runs
(dashed and solid blue curves). The good overlap indicates that the features
of energy dissipation remain unchanged when including more Hermite moments and
thus extending the range of spatial scales over which a potential plasma echo
could arise. In the same figure we follow Eq. (21) and plot the dependence of
the critical Hermite mode number $m_{\rm cr}$ on $k_{\perp}\rho_{i}$ (red
curve). Based on the comparison between advection and phase-mixing rates (Eq.
(20)), the plasma echo is allowed in the region above the red curve and below
the cutoff $M$, which is a wide region for both runs. Therefore, the plasma
echo is absent not because it is excluded in our simulations due to limited
resolution in position or velocity space; instead, its absence must be due to
the turbulent dynamics.

Figure 8: The ratio of the nonlinear-advection rate to phase-mixing rate of
$g_{7}$ as a function of normalized current density $|J_{z}/J_{\rm rms}|$. The
color corresponds to the local energy density of $g_{7}$.

We proposed in Sec. III.2 that the absence of echo is due to the local
reduction of the advective nonlinearity around current sheets. The visuals of
$J_{z}$ (Fig. 3, top panels) show that indeed, at scales larger than $d_{e}$,
the regions with strong current density are elongated “sheet” structures. In
Fig. 8 we show the ratio of the nonlinear-advection rate to the phase-mixing
rate of $g_{7}$ versus the normalized current density (the results using other
Hermite moments are qualitatively similar); the color map of the data points
represents the local energy density of $g_{7}$. This plot clearly shows that
nonlinear advection is subdominant with respect to phase mixing at positions
with large current density. We conjecture that it is the specific
configuration of $A_{z}$, $\varphi$, and $g_{m}$ around the sheet-like
structure of current that leads to the local suppression of echo effects; this
interpretation is inspired by observations of strong phase-mixing in
supposedly similar geometries that arise in reconnection studies
[loureiro2013fast, numata2015, mccubbin2022characterizing]. (Note that simply
invoking the quasi-1D (in the perpendicular plane) morphology of current
sheets is not sufficient to explain the weakening of the advection-to-phase-
mixing rates, because such weakening would affect equally the advective
nonlinearity and the parallel streaming.) As a result of enhanced phase-
mixing, current sheets are also where the $g_{m}$ fluctuations are
concentrated (see Eq. (LABEL:eq:gm_solution)), as shown by the correlation
between $|g_{7}|^{2}$ and $|J_{z}/J_{\rm rms}|$ in Fig. 8. Therefore, current
sheets are energetically important for dissipation. The overlap in position
space between the suppression of plasma echo and the concentration of free
energy in $g_{e}$ results in the efficient energy transfer to small scales in
velocity space and, eventually, strong electron heating at the kinetic range.

To further visualize this correlation between the current density and electron
heating, in the color map of current density (Fig. 3 top panels), we show the
regions with strong electron heating using the contour (magenta lines) of the
electron heating rate (measured by the energy dissipation through hyper-
collisions at small scales in velocity space) that is twice its volume-
averaged value. Indeed we find that, for scales such that $\rho_{i}\simeq
k_{\perp}^{-1}>d_{e}$ (left and middle panels), significant electron heating
mostly occurs around current sheets, consistent with the above arguments. Such
correlation between heating and current sheets has been reported in numerical
simulations of kinetic turbulence [tenbarge2013current, navarro2016structure]
and observations [carbone2022high]. Here, we add to the previous body of work
by showing how such heating is possible: the reduction of the advective
nonlinearity around current sheets enables phase-mixing to proceed unimpeded
in those locations.

##  VI Conclusions

In this work, we conduct an analytical and numerical study of sub-ion-
gyroscale (sub-$\rho_{i}$) turbulence in the low-$\beta$ limit using an
analytical formalism known as the Kinetic Reduced Electron Heating Model
(KREHM) [zocco2011reduced]. Two main issues addressed by this work concern the
energy spectra of sub-$\rho_{i}$ turbulence, which is of a kinetic Alfvén wave
(KAW) nature in this low-$\beta$ regime, and the energy dissipation (and
electron heating) in such turbulence.

In the first part of the paper, we study the energy spectrum, spectral
anisotropy, and intermittency of the turbulent fluctuations in the limit of
isothermal electrons. In our numerical results, the unchanged
$k_{\perp}^{-8/3}$ magnetic (and density) spectrum obtained from simulations
with different flux-unfreezing mechanisms (electron inertia versus
(hyper-)resistivity) suggests that the spectrum is not set by tearing
mediation. By calculating the 3D, multi-point structure function of the
magnetic fields, we obtain a $\propto k_{\parallel}^{-7/2}$ parallel spectrum,
and the scalings $\ell\propto\lambda^{2/3}$ and $\xi\propto\lambda$, which are
in good agreement with recent solar wind measurements [wang2020observational,
zhang2022three]. The measured $\xi\propto\lambda$ scaling indicates isotropy
of the fluctuations in the plane perpendicular to the (local) magnetic field,
and is consistent with the absence of tearing mediation. We argue that,
instead, the steepening of the spectrum with respect to the straightforward
KAW cascade prediction of $k_{\perp}^{-7/3}$ is likely to be due to
intermittency, which we quantify using the higher $q$th-order structure
function $S_{q}(\delta r)$. The measured $S_{q}(\delta r)\sim\delta r^{q/3+1}$
scaling reveals that the energy-containing, intermittent structures are
quasi-2D, which confirms the main assumption made by boldyrev2012spectrum to
derive the $k_{\perp}^{-8/3}$ spectrum.

In the second part of the paper, we take into account electron kinetic effects
and study the phase-space dynamics and channels for energy dissipation. In the
$m$-$k_{\perp}$ phase space, where $m$ is the velocity-space mode number, the
energy dissipation is found to mostly occur at high $m$ through collisions,
implying that phase mixing dominates nonlinear advection of the free energy
and, thus, that the “kinetic channel” (as opposed to the “fluid channel”) is
the main route for energy dissipation and electron heating. We verify that
this argument is quantitatively supported, as follows: if we subtract from the
magnetic spectrum pertaining to an isothermal run the energy that would be
linearly Landau-damped within one eddy turn-over-time, scale-by-scale, we
obtain a spectrum that agrees remarkably well with that from a corresponding
run with kinetic electrons. A direct consequence of this result is that the
application of a Landau-fluid-type closure [Hammett1990] to KAW turbulence in
the low-$\beta$ limit is justified.

In Sec. III.3, we derive a lowest-order solution for the coefficients of the
Hermite expansion of the electron distribution function, $g_{m}$, under the
assumption that phase-mixing dominates over nonlinear advection:
$g_{m+1}=-\sqrt{m/(m+1)}g_{m-1}$, and $g_{3}\propto J_{z}$. This solution, as
well as the corresponding $m^{-1/2}$ Hermite spectrum, is confirmed by our
numerical results. Lastly, we find that the absence of stochastic plasma echo,
which would impede phase mixing, is due to the ubiquitous current sheets that
develop self-consistently in our system. The fluctuations of the Hermite
moments are concentrated in the vicinity of current sheets, while at these
specific positions, the advective nonlinearity is weakened because of the
specific configurations of the magnetic field and the flow around the sheets.
This observed anti-correlation between advective nonlinearity and current
density (and thus the density of free-energy) undermines plasma echo and
allows the free energy to transfer to higher $m$ by phase mixing.

It is important to note that for all our simulations, the long-wavelength MHD
range is either not included or too narrow for the progressive alignment of
the fluctuations [boldyrev2006spectrum] to establish itself. In reality,
fluctuations cascading from large scales down to the $\rho_{i}$ scale (i.e.,
in the MHD range) are expected to become progressively more anisotropic in the
field-perpendicular plane, and eventually transition to a tearing-mediated
range at a critical scale $\lambda_{\rm cr}>\rho_{i}$
[loureiro2017collisionless]. Below $\lambda_{\rm cr}$, the tearing mediation
is expected to decrease the alignment and isotropize the eddies (see Eq. (22)
in [boldyrev2017magnetohydrodynamic] and Sec. 6 in [mallet2017disruption]).
Therefore, the fluctuations entering the sub-$\rho_{i}$ range can be isotropic
only when there is a sufficiently wide range between $\lambda_{\rm cr}$ and
$\rho_{i}$. In weakly collisional turbulence, $\lambda_{\rm cr}$ is expected
to be close to $\rho_{i}$ (see Eqs. (16) and (18) in
[loureiro2017collisionless] and Sec. 4 in [mallet2017disruptionb]), and so the
fluctuations in the sub-$\rho_{i}$ range can be anisotropic. However, our
numerical results suggests that there is no dynamic alignment effect at
sub-$\rho_{i}$ scales, i.e., the residual anisotropy from the MHD range would
not increase in this range. Therefore, it is possible that the fluctuations
around the $\rho_{i}$ scale are sufficiently anisotropic for tearing
mediation; but the tearing mode should become stable as $k_{\perp}\rho_{i}$
increases. With greater computational resources than currently available,
simulations might be possible that cover both the MHD and the kinetic ranges
with sufficient resolution such as to enable a more self-consistent study of
turbulence, including the microscale feedback to the long-wavelength range
(via tearing), and a more accurate “large-scale” condition for the
sub-$\rho_{i}$ turbulence. Having said this, the agreement between our results
and several pieces of observational evidence is reassuring.

The efficient electron heating we found in the low-$\beta$ limit is
qualitatively consistent with previous numerical studies on electron versus
ion energization in kinetic turbulence, in which the electron to ion heating-
rate-ratio is a decreasing function of plasma $\beta$ (e.g.,
[howes2010prescription, kawazura2019thermal, zhdankin2019,
arzamasskiy2019hybrid]). We believe that the weakening of the advective
nonlinearity around current sheets — which enables phase mixing and heating —
remains valid beyond the low-$\beta$ limit that we explore here. In a
gyrokinetic study of reconnection [numata2015], it is found that ion heating
(via nonlinear phase-mixing) becomes progressively more important as $\beta$
increases. Extrapolating these ideas to kinetic turbulence suggests that, as
beta increases (but still below unity), current sheets remain critical as
energy dissipation sites, with the balance shifting from heating electrons at
low beta (via linear phase mixing) to ions (via nonlinear phase mixing).
Further studies are required to test this hypothesis.

In addition to critical plasma parameters such as the plasma $\beta$ and
temperature ratio, the effects of the “large-scale” conditions of the
turbulence on the energy dissipation are also to be investigated. For example,
the imbalance of turbulence relevant to the fast-wind streams has recently
been found to have profound effects on the heating mechanisms for ions and
electrons [meyrand2021violation, squire2022high], while the electron kinetic
effects are omitted. Kinetic studies taking into account both electron and ion
physics with relevant large-scale conditions are necessary to explain the
extensive measurements taken both at corona and the solar wind by currently
operating and future spacecrafts.

In summary, this paper presents a self-consistent explanation for the
electromagnetic spectra and electron heating in low-$\beta$ collisionless
turbulence by accounting for entwined physical processes in position and
velocity spaces. On the one hand, we show that intermittency and electron
Landau damping are the two physical mechanisms that underlie the steepening of
the spectrum with respect to the prediction from a Kolmogorov-like KAW
cascade. On the other, it is due to the presence of these spontaneously formed
sheet-like structures that the kinetic effects of phase mixing and Landau
damping are locally enhanced and thus play an important role in energy
dissipation. By explaining the above underling physics that connects the
position-space and velocity-space dynamics, this study clarifies the long-
standing discussion of the relation between current sheets, magnetic
reconnection, and electron heating in the low-$\beta$ kinetic turbulence.

#### Acknowledgments.

The authors thank L. Arzamasskiy, S. Boldyrev, G. Howes, N. Mandell, A. A.
Schekochihin, V. Zhdankin for insightful discussions. This work was supported
by the National Science Foundation (NSF) under CAREER award No. 1654168 (NFL
and MZ), by the National Aeronautics and Space Administration (NASA) under
award NNH19ZA001N-FINESST (MZ), and by the NSF-DOE Partnership in Basic Plasma
Science and Engineering Award No. PHY-2010136 (ZL). This research used
resources of the MIT-PSFC partition of the Engaging cluster at the MGHPCC
facility, funded by DOE award No. DE-FG02-91-ER54109 and the National Energy
Research Scientific Computing Center, a DOE Office of Science User Facility
supported by the Office of Science of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231 using NERSC award FES-ERCAP0020063.

## References

  * Borovsky and Funsten [2003] J. E. Borovsky and H. O. Funsten, Journal of Geophysical Research: Space Physics 108 (2003).
  * Goldreich and Sridhar [1995] P. Goldreich and S. Sridhar, The Astrophysical Journal 438, 763 (1995).
  * Tu and Marsch [1995] C.-Y. Tu and E. Marsch, Space Science Reviews 73, 1 (1995).
  * Parker [1983] E. N. Parker, The Astrophysical Journal 264, 635 (1983).
  * Balbus and Hawley [1998] S. A. Balbus and J. F. Hawley, Reviews of Modern Physics 70, 1 (1998).
  * Begelman and Fabian [1990] M. C. Begelman and A. C. Fabian, Monthly Notices of the Royal Astronomical Society 244, 26P (1990).
  * Sahraoui  _et al._ [2009] F. Sahraoui, M. L. Goldstein, P. Robert, and Y. V. Khotyaintsev, Physical Review Letters 102, 231102 (2009).
  * Alexandrova  _et al._ [2009] O. Alexandrova, J. Saur, C. Lacombe, A. Mangeney, J. Mitchell, S. J. Schwartz, and P. Robert, Physical Review Letters 103, 165003 (2009).
  * Alexandrova  _et al._ [2012] O. Alexandrova, C. Lacombe, A. Mangeney, R. Grappin, and M. Maksimovic, The Astrophysical Journal 760, 121 (2012).
  * Alexandrova  _et al._ [2013] O. Alexandrova, C. H. K. Chen, L. Sorriso-Valvo, T. S. Horbury, and S. D. Bale, Space Science Reviews 178, 101 (2013).
  * Chen  _et al._ [2019] C. H. K. Chen, K. G. Klein, and G. G. Howes, Nature Communications 10, 1 (2019).
  * Gary and Nishimura [2004] S. P. Gary and K. Nishimura, Journal of Geophysical Research: Space Physics 109 (2004).
  * Howes  _et al._ [2008a] G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin, Journal of Geophysical Research: Space Physics 113 (2008a).
  * Howes  _et al._ [2008b] G. G. Howes, W. Dorland, S. C. Cowley, G. W. Hammett, E. Quataert, A. A. Schekochihin, and T. Tatsuno, Physical Review Letters 100, 065004 (2008b).
  * Salem  _et al._ [2012] C. S. Salem, G. G. Howes, D. Sundkvist, S. D. Bale, C. C. Chaston, C. H. K. Chen, and F. S. Mozer, The Astrophysical Journal Letters 745, L9 (2012).
  * TenBarge  _et al._ [2012] J. M. TenBarge, J. J. Podesta, K. G. Klein, and G. G. Howes, The Astrophysical Journal 753, 107 (2012).
  * Franci  _et al._ [2018] L. Franci, S. Landi, A. Verdini, L. Matteini, and P. Hellinger, The Astrophysical Journal 853, 26 (2018).
  * Cerri  _et al._ [2019] S. S. Cerri, D. Grošelj, and L. Franci, Frontiers in Astronomy and Space Sciences , 64 (2019).
  * Grošelj  _et al._ [2019] D. Grošelj, C. H. K. Chen, A. Mallet, R. Samtaney, K. Schneider, and F. Jenko, Physical Review X 9, 031037 (2019).
  * Leamon  _et al._ [1998] R. J. Leamon, C. W. Smith, N. F. Ness, W. H. Matthaeus, and H. K. Wong, Journal of Geophysical Research: Space Physics 103, 4775 (1998).
  * Podesta  _et al._ [2010] J. J. Podesta, J. E. Borovsky, and S. P. Gary, The Astrophysical Journal 712, 685 (2010).
  * Chen  _et al._ [2013] C. H. K. Chen, S. Boldyrev, Q. Xia, and J. C. Perez, Physical Review Letters 110, 225002 (2013).
  * Kiyani  _et al._ [2012] K. H. Kiyani, S. C. Chapman, F. Sahraoui, B. Hnat, O. Fauvarque, and Y. V. Khotyaintsev, The Astrophysical Journal 763, 10 (2012).
  * Chen [2016] C. H. K. Chen, Journal of Plasma Physics 82 (2016).
  * Cho and Lazarian [2004] J. Cho and A. Lazarian, The Astrophysical Journal Letters 615, L41 (2004).
  * Cho and Lazarian [2009] J. Cho and A. Lazarian, The Astrophysical Journal 701, 236 (2009).
  * Schekochihin  _et al._ [2009] A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett, G. G. Howes, E. Quataert, and T. Tatsuno, The Astrophysical Journal Supplement Series 182, 310 (2009).
  * Bale  _et al._ [2005] S. D. Bale, P. J. Kellogg, F. S. Mozer, T. S. Horbury, and H. Reme, Physical Review Letters 94, 215002 (2005).
  * Servidio  _et al._ [2015] S. Servidio, F. Valentini, D. Perrone, A. Greco, F. Califano, W. H. Matthaeus, and P. Veltri, Journal of Plasma Physics 81 (2015).
  * Wan  _et al._ [2015a] M. Wan, W. H. Matthaeus, V. Roytershteyn, H. Karimabadi, T. Parashar, P. Wu, and M. Shay, Physical Review Letters 114, 175002 (2015a).
  * Told  _et al._ [2015] D. Told, F. Jenko, J. M. TenBarge, G. G. Howes, and G. Hammett, Physical Review Letters 115, 025003 (2015).
  * Cerri  _et al._ [2018] S. S. Cerri, M. W. Kunz, and F. Califano, The Astrophysical Journal Letters 856, L13 (2018).
  * Arzamasskiy  _et al._ [2019] L. Arzamasskiy, M. W. Kunz, B. D. G. Chandran, and E. Quataert, The Astrophysical Journal 879, 53 (2019).
  * Grošelj  _et al._ [2018] D. Grošelj, A. Mallet, N. F. Loureiro, and F. Jenko, Physical Review Letters 120, 105101 (2018).
  * Kiyani  _et al._ [2009] K. H. Kiyani, S. C. Chapman, Y. V. Khotyaintsev, M. W. Dunlop, and F. Sahraoui, Physical Review Letters 103, 075006 (2009).
  * Sahraoui  _et al._ [2013] F. Sahraoui, S. Y. Huang, G. Belmont, M. L. Goldstein, A. Retinò, P. Robert, and J. De Patoul, The Astrophysical Journal 777, 15 (2013).
  * Sundkvist  _et al._ [2007] D. Sundkvist, A. Retinò, A. Vaivads, and S. D. Bale, Physical Review Letters 99, 025004 (2007).
  * Karimabadi  _et al._ [2011] H. Karimabadi, J. Dorelli, V. Roytershteyn, W. Daughton, and L. Chacón, Physical Review Letters 107, 025002 (2011).
  * Osman  _et al._ [2014] K. T. Osman, W. H. Matthaeus, J. T. Gosling, A. Greco, S. Servidio, B. Hnat, S. C. Chapman, and T. D. Phan, Physical Review Letters 112, 215002 (2014).
  * Wan  _et al._ [2015b] M. Wan, W. H. Matthaeus, V. Roytershteyn, H. Karimabadi, T. Parashar, P. Wu, and M. Shay, Physical Review Letters 114, 175002 (2015b).
  * Chasapis  _et al._ [2015] A. Chasapis, A. Retinò, F. Sahraoui, A. Vaivads, Y. V. Khotyaintsev, D. Sundkvist, A. Greco, L. Sorriso-Valvo, and P. Canu, The Astrophysical Journal Letters 804, L1 (2015).
  * Camporeale  _et al._ [2018] E. Camporeale, L. Sorriso-Valvo, F. Califano, and A. Retinò, Physical Review Letters 120, 125101 (2018).
  * Boldyrev and Perez [2012] S. Boldyrev and J. C. Perez, The Astrophysical Journal Letters 758, L44 (2012).
  * Boldyrev [2006] S. Boldyrev, Physics Review Letters 96, 115002 (2006).
  * Loureiro and Boldyrev [2017a] N. F. Loureiro and S. Boldyrev, Physical Review Letters 118, 245101 (2017a).
  * Mallet and Schekochihin [2017] A. Mallet and A. A. Schekochihin, Monthly Notices of the Royal Astronomical Society 466, 3918 (2017).
  * Boldyrev and Loureiro [2017] S. Boldyrev and N. F. Loureiro, The Astrophysical Journal 844, 125 (2017).
  * Dong  _et al._ [2018] C. Dong, L. Wang, Y.-M. Huang, L. Comisso, and A. Bhattacharjee, Physical Review Letters 121, 165101 (2018).
  * Walker  _et al._ [2018] J. Walker, S. Boldyrev, and N. F. Loureiro, Physical Review E 98, 033209 (2018).
  * Loureiro and Boldyrev [2020] N. F. Loureiro and S. Boldyrev, The Astrophysical Journal 890, 55 (2020).
  * Loureiro and Boldyrev [2017b] N. F. Loureiro and S. Boldyrev, The Astrophysical Journal 850, 182 (2017b).
  * Parashar  _et al._ [2015] T. N. Parashar, C. Salem, R. T. Wicks, H. Karimabadi, S. P. Gary, and W. H. Matthaeus, Journal of Plasma Physics 81 (2015).
  * Dmitruk  _et al._ [2004] P. Dmitruk, W. H. Matthaeus, and N. Seenu, The Astrophysical Journal 617, 667 (2004).
  * Matthaeus and Velli [2011] W. H. Matthaeus and M. Velli, Space Science Reviews 160, 145 (2011).
  * Servidio  _et al._ [2011] S. Servidio, P. Dmitruk, A. Greco, M. Wan, S. Donato, P. A. Cassak, M. A. Shay, V. Carbone, and W. H. Matthaeus, Nonlinear Processes in Geophysics 18, 675 (2011).
  * Osman  _et al._ [2011] K. T. Osman, W. H. Matthaeus, A. Greco, and S. Servidio, The Astrophysical Journal Letters , L11 (2011).
  * Zhdankin  _et al._ [2013] V. Zhdankin, D. A. Uzdensky, J. C. Perez, and S. Boldyrev, The Astrophysical Journal 771, 124 (2013).
  * Chapman and Cowling [1990] 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_ (Cambridge university press, 1990).
  * Schekochihin  _et al._ [2008] A. A. Schekochihin, S. C. Cowley, W. Dorland, G. W. Hammett, G. G. Howes, G. G. Plunk, E. Quataert, and T. Tatsuno, Plasma Physics and Controlled Fusion 50, 124024 (2008).
  * Navarro  _et al._ [2016] A. B. Navarro, B. Teaca, D. Told, D. Groselj, P. Crandall, and F. Jenko, Physical Review Letters 117, 245101 (2016).
  * Grošelj  _et al._ [2017] D. Grošelj, S. S. Cerri, A. B. Navarro, C. Willmott, D. Told, N. F. Loureiro, F. Califano, and F. Jenko, The Astrophysical Journal 847, 28 (2017).
  * Servidio  _et al._ [2017] S. Servidio, A. Chasapis, W. Matthaeus, D. Perrone, F. Valentini, T. Parashar, P. Veltri, D. Gershman, C. Russell, B. Giles, _et al._ , Physical Review Letters 119, 205101 (2017).
  * Schekochihin  _et al._ [2016] A. A. Schekochihin, J. T. Parker, E. G. Highcock, P. J. Dellar, W. Dorland, and G. W. Hammett, Journal of Plasma Physics 82 (2016).
  * Adkins and Schekochihin [2018] T. Adkins and A. A. Schekochihin, Journal of Plasma Physics 84 (2018).
  * Eyink [2018] G. L. Eyink, Physical Review X 8, 041020 (2018).
  * Gould  _et al._ [1967] R. W. Gould, T. M. O’Neil, and J. H. Malmberg, Physical Review Letters 19, 219 (1967).
  * Malmberg  _et al._ [1968] J. H. Malmberg, C. B. Wharton, R. W. Gould, and T. M. O’Neil, Physical Review Letters 20, 95 (1968).
  * Meyrand  _et al._ [2019] R. Meyrand, A. Kanekar, W. Dorland, and A. A. Schekochihin, Proceedings of the National Academy of Sciences 116, 1185 (2019).
  * Parker  _et al._ [2016] J. T. Parker, E. Highcock, A. A. Schekochihin, and P. Dellar, Physics of Plasmas 23, 070703 (2016).
  * Zocco and Schekochihin [2011] A. Zocco and A. A. Schekochihin, Physics of Plasmas 18, 102309 (2011).
  * Kadomtsev and Pogutse [1973] B. B. Kadomtsev and O. P. Pogutse, Soviet physics, Journal of Experimental and Theoretical Physics 5, 575 (1973).
  * Strauss [1976] H. R. Strauss, The Physics of Fluids 19, 134 (1976).
  * Zank and Matthaeus [1992] G. P. Zank and W. H. Matthaeus, Journal of Plasma Physics 48, 85 (1992).
  * Note [1] In the conventional $\beta$-model [frisch1995turbulence] for intermittency widely used in hydrodynamics, the dimension of intermittent structures is a fractal dimension and is not necessarily the same as their spatial dimension. However, in boldyrev2012spectrum and other turbulence simulations, it is found that intermittent structures do appear to be 2D in position space.
  * Frisch [1995] U. Frisch,  _Turbulence: the legacy of A. N. Kolmogorov_ (Cambridge university press, 1995).
  * Mallet  _et al._ [2017a] A. Mallet, A. A. Schekochihin, and B. D. G. Chandran, Monthly Notices of the Royal Astronomical Society 468, 4862 (2017a).
  * Note [2] For the Harris sheet-type configuration, the growth rate of the fastest tearing mode becomes $\gamma_{t}\sim v_{A\lambda}d_{e}\rho_{s}/\lambda^{3}$ [loureiro2017collisionless], leading to the scale dependence of magnetic fluctuation $v_{A\lambda}\sim\varepsilon^{1/3}\lambda(d_{e}\rho_{s})^{-1/3}$, the magnetic spectrum $E_{B}(k_{\perp})\sim 4\pi\rho_{0}\varepsilon^{2/3}(d_{e}\rho_{s})^{-2/3}k_{\perp}^{-3}dk_{\perp}$, and the spectral anisotropy $k_{\parallel}\sim\varepsilon^{1/3}(d_{e}^{2}/\rho_{s})^{1/3}V_{A}^{-1}k_{\perp}$.
  * Howes  _et al._ [2018] G. G. Howes, A. J. McCubbin, and K. G. Klein, Journal of Plasma Physics 84 (2018).
  * Klein and Howes [2016] K. G. Klein and G. G. Howes, The Astrophysical Journal Letters 826, L30 (2016).
  * Howes  _et al._ [2017] G. G. Howes, K. G. Klein, and T. C. Li, Journal of Plasma Physics 83 (2017).
  * Loureiro  _et al._ [2013] N. F. Loureiro, A. A. Schekochihin, and A. Zocco, Physical Review Letters 111, 025002 (2013).
  * Numata and Loureiro [2015] R. Numata and N. F. Loureiro, Journal of Plasma Physics 81, 305810201 (2015).
  * Kanekar  _et al._ [2015] A. Kanekar, A. A. Schekochihin, W. Dorland, and N. F. Loureiro, Journal of Plasma Physics 81 (2015).
  * Loureiro  _et al._ [2016] N. F. Loureiro, W. Dorland, L. Fazendeiro, A. Kanekar, A. Mallet, M. S. Vilelas, and A. Zocco, Computer Physics Communications 206, 45 (2016).
  * Mallet  _et al._ [2016] A. Mallet, A. A. Schekochihin, B. D. G. Chandran, C. H. K. Chen, T. S. Horbury, R. T. Wicks, and C. C. Greenan, Monthly Notices of the Royal Astronomical Society 459, 2130 (2016).
  * Schekochihin [2020] A. A. Schekochihin, arXiv preprint arXiv:2010.00699 (2020).
  * Wang  _et al._ [2020] T. Wang, J. He, O. Alexandrova, M. Dunlop, and D. Perrone, The Astrophysical Journal 898, 91 (2020).
  * Zhang  _et al._ [2022] J. Zhang, S. Y. Huang, J. S. He, T. Y. Wang, Z. G. Yuan, X. H. Deng, Y. Y. Jiang, K .and Wei, S. B. Xu, Q. Y. Xiong, _et al._ , The Astrophysical Journal Letters 924, L21 (2022).
  * Loureiro  _et al._ [2007] N. F. Loureiro, A. A. Schekochihin, and S. C. Cowley, Physics of Plasmas 14, 100703 (2007).
  * Note [3] The large-$m$ dissipation at the system scale is artificial, caused by the forcing.
  * Howes  _et al._ [2011] G. G. Howes, J. M. Tenbarge, and W. Dorland, Physics of Plasmas 18, 102305 (2011).
  * TenBarge  _et al._ [2013] J. M. TenBarge, G. G. Howes, and W. Dorland, The Astrophysical Journal 774, 139 (2013).
  * TenBarge and Howes [2013] J. M. TenBarge and G. G. Howes, The Astrophysical Journal Letters 771, L27 (2013).
  * Howes [2010] G. G. Howes, Monthly Notices of the Royal Astronomical Society: Letters 409, L104 (2010).
  * Passot and Sulem [2015] T. Passot and P. L. Sulem, The Astrophysical Journal Letters 812, L37 (2015).
  * Kunz  _et al._ [2018] M. W. Kunz, I. G. Abel, K. G. Klein, and A. A. Schekochihin, Journal of Plasma Physics 84 (2018).
  * McCubbin  _et al._ [2022] A. J. McCubbin, G. G. Howes, and J. M. TenBarge, Physics of Plasmas 29, 052105 (2022).
  * Carbone  _et al._ [2022] V. Carbone, D. Telloni, F. Lepreti, and A. Vecchio, The Astrophysical Journal Letters 924, L26 (2022).
  * Hammett and Perkins [1990] G. W. Hammett and F. W. Perkins, Physical Review Letters 64, 3019 (1990).
  * Mallet  _et al._ [2017b] A. Mallet, A. A. Schekochihin, and B. D. G. Chandran, Journal of Plasma Physics 83 (2017b).
  * Kawazura  _et al._ [2019] Y. Kawazura, M. Barnes, and A. A. Schekochihin, Proceedings of the National Academy of Sciences 116, 771 (2019).
  * Zhdankin  _et al._ [2019] V. Zhdankin, D. A. Uzdensky, G. R. Werner, and M. C. Begelman, Phys. Rev. Lett. 122, 055101 (2019).
  * Meyrand  _et al._ [2021] R. Meyrand, J. Squire, A. A. Schekochihin, and W. Dorland, Journal of Plasma Physics 87 (2021).
  * Squire  _et al._ [2022] J. Squire, R. Meyrand, M. W. Kunz, L. Arzamasskiy, A. A. Schekochihin, and E. Quataert, Nature Astronomy , 1 (2022).
  * Coppi  _et al._ [1976] B. Coppi, R. Galvao, R. Pellat, M. N. Rosenbluth, and P. Rutherford, Soviet Journal of Plasma Physics 2, 533 (1976).
  * Furth  _et al._ [1963] H. P. Furth, J. Killeen, and M. N. Rosenbluth, The physics of Fluids 6, 459 (1963).
  * Cho and Vishniac [2000] J. Cho and E. T. Vishniac, The Astrophysical Journal 539, 273 (2000).

## Supplemental Materials

###  .1 Linear tearing mode with (hyper-)resistivity

In Sec. III.1 we invoked scalings pertaining to the (kinetic) tearing
instability driven by (hyper-)resistivity. We derive them here. We perform a
linear calculation of the tearing mode with an initial equilibrium described
by $\mbox{\boldmath{$B$}}_{\perp 0}(x)=\bar{B}_{\perp
0}f(x/a)\mbox{\boldmath{$\hat{y}$}}$ where $f(x/a)$ is an odd function
satisfying $f(0)=0$ and $a$ is the characteristic length scale of the gradient
of the equilibrium magnetic field. We consider the case where $\rho_{i}\gg a$
and so the GK Poisson’s law reduces to Eq. (1). We assume there is no
equilibrium flow. In this setup, the linearized KREHM equations in the
isothermal limit (i.e., $g_{e}=0$) become

$\gamma n_{e}=ik_{y}B_{\perp
0}(\partial_{x}^{2}-k_{y}^{2})A_{\parallel}-ik_{y}A_{\parallel}\bar{B}_{\perp
0}f^{\prime\prime}(x/a),$ (S1)
$\gamma(1+k_{\perp}^{2}d_{e}^{2})A_{\parallel}-ik_{y}\varphi B_{\perp
0}=-\rho_{s}^{2}ik_{y}n_{e}B_{\perp 0}+\eta_{H}k_{\perp}^{2}A_{\parallel}.$
(S2)

The tearing instability parameter
$\Delta^{\prime}\equiv\left[dA_{\parallel}/dx\right]_{0-}^{0+}$ is obtained by
solving the outer layer (MHD) equations: for $k_{y}a\ll 1$, it is
$\Delta^{\prime}a\propto(k_{y}a)^{-n}$, where $n=1$ corresponds to a Harris
sheet-like configuration for the equilibrium, and $n=2$ corresponds to a
sinusoidal profile.

In the inner layer, where $x\ll a$ and $\partial_{x}^{2}\gg k_{y}^{2}$, the
equilibrium magnetic field can be approximated as $B_{\perp
0}\approx\bar{B}_{\perp 0}x/a$. In the large $\Delta^{\prime}$ limit (the
“Coppi” modes [coppi1976resistive]), $\partial_{x}^{2}A_{\parallel}\approx
A_{\parallel}/\delta^{2}$, where $\delta$ is the thickness of the inner
boundary layer. In the opposite, small $\Delta^{\prime}$ limit (the “FKR”
modes [furth1963finite]),
$\partial_{x}^{2}A_{\parallel}\approx\Delta^{\prime}A_{\parallel}/\delta$. Let
us first consider the latter case. The balance of terms in the inner layer
equations becomes:

$\displaystyle\gamma
n_{e}\sim\bar{B_{0}}\frac{\delta}{a}k_{y}\frac{\Delta^{\prime}A_{\parallel}}{\delta}\sim\frac{\bar{B}_{\perp
0}}{a^{2}}k_{y}A_{\parallel}(k_{y}a)^{-n},$ (S3)
$\gamma(1+\frac{\Delta^{\prime}}{\delta}d_{e}^{2})A_{\parallel}\sim
k_{y}\rho_{\tau}^{2}\bar{B}_{\perp
0}\frac{\delta}{a}n_{e}-\eta_{H}\frac{\Delta^{\prime}}{\delta^{\alpha-1}}A_{\parallel},$
(S4)

where $\rho_{\tau}=\rho_{s}^{2}+\rho_{i}^{2}$. In the case of $\eta_{H}$=0,
frozen flux is broken by the electron inertia. Balancing the two terms on the
l.h.s. of Ohm’s law (Eq. (S4)) yields an expression for the width of the inner
layer:

$\displaystyle\delta\sim\Delta^{\prime}d_{e}^{2}\sim(k_{y}a)^{-n}d_{e}^{2}/a.$
(S5)

Using this expression as well as Eq. (S3) and Eq. (S4), we obtain the scaling
of the growth rate for FKR-type modes:

$\gamma\sim k_{y}\rho_{\tau}(k_{y}a)^{-n}\frac{d_{e}}{a}\frac{\bar{B}_{\perp
0}}{a}.$ (S6)

Let us now consider the large $\Delta^{\prime}$ limit, where
$\Delta^{\prime}a\sim(k_{y}a)^{-n}\rightarrow a/\delta$, and the growth rate
and inner layer thickness become

$\gamma\sim k_{y}\frac{\rho_{\tau}}{a}\frac{d_{e}}{a}\bar{B}_{\perp
0}\frac{d_{e}}{a}\sim k_{y}\frac{\rho_{\tau}}{a}\bar{B}_{\perp 0},$ (S7)
$\delta^{2}\sim a\frac{d_{e}^{2}}{a}\rightarrow\delta\sim d_{e}.$ (S8)

The growth rate of the most unstable mode is obtained by matching the growth
rate in the large and the small $\Delta^{\prime}$ limits:

$k_{y}\frac{\rho_{\tau}}{a}\bar{B}_{\perp 0}\frac{d_{e}}{a}(k_{y}a)^{-n}\sim
k_{y}\frac{\rho_{\tau}}{a}\bar{B}_{\perp 0},$ (S9)

from which follows that the most unstable wavenumber scales as

$(k_{\rm max}a)^{n}\sim d_{e}/a\rightarrow k_{\rm max}a\sim(d_{e}/a)^{1/n}.$
(S10)

We thus obtain

$\gamma_{\rm max}\sim\frac{\bar{B}_{\perp
0}}{a}\frac{\rho_{\tau}}{a}\left(\frac{d_{e}}{a}\right)^{1/n}.$ (S11)

In the case where the frozen flux is broken by the (hyper or Laplacian)
resistivity, we can formally replace the electron inertia term with the
resistivity term in the above derivation:

$\gamma\frac{\Delta^{\prime}}{\delta}d_{e}^{2}\rightarrow\eta_{H}\frac{\Delta^{\prime}}{\delta^{\alpha-1}},\quad\text{i.e.},\quad
d_{e}^{2}\rightarrow\eta_{H}\gamma^{-1}\delta^{2-\alpha},$ (S12)

where $\eta_{H}$ is the value of the resistivity and $\alpha$ is the order of
the spatial derivative in the resistive term. This then yields:

$\gamma_{\rm max}^{(2n+1)/2n}\sim\frac{\bar{B}_{\perp
0}}{a}\frac{\rho_{\tau}}{a}a^{-1/n}\eta_{H}^{1/2n}\delta^{(2-\alpha)/2n}.$
(S13)

The thickness of the inner layer should be

$\delta\sim d_{e}\sim(\eta_{H}/\gamma_{\rm
max})^{1/2}\delta^{(2-\alpha)/2}\rightarrow\delta\sim(\eta_{H}/\gamma_{\rm
max})^{1/\alpha}.$ (S14)

Combining Eqs. (S13) and (S14), we obtain the final expressions for
$\gamma_{\rm max}$ and $\delta$ as a function of the parameter for the
equilibrium profile $n$, and that for the order of the (hyper-)dissipation
$\alpha$:

$\gamma_{\rm max}\sim\left(\bar{B}_{\perp
0}\rho_{\tau}\right)^{n\alpha/(n\alpha+1)}a^{-\alpha(2n+1)/(n\alpha+1)}\eta_{H}^{1/(n\alpha+1)},$
(S15)

and

$\delta\sim\left(\bar{B}_{\perp
0}\rho_{\tau}\right)^{-n/(n\alpha+1)}a^{(2n+1)/(n\alpha+1)}\eta_{H}^{n/(n\alpha+1)}.$
(S16)

In the model of tearing-mediated turbulence, the nonlinear timescale is set by
the growth rate of tearing instability, expressed by Eq. (S15) (or Eq. (S11)
in the case of electron inertia). We use the expression for $\delta$
(Eq.(S16)) to estimate the resolution in the numerical simulation required to
resolve the inner layer thickness of the tearing mode set by the electron
inertia or the (hyper-) resistivity.

We now derive the magnetic (or density) energy spectrum at sub-$\rho_{i}$
scales following the tearing-mediated turbulence model (described in Sec.
III.1) for the case where the frozen flux is broken by the resistivity
(instead of the electron inertia). In that model, the eddy-turnover rate,
which satisfies the dimensional relation $\gamma_{nl}\sim\varepsilon/(\delta
B_{\perp\lambda}^{2}/8\pi)$, is set by the growth rate of tearing mode (Eq.
(S15)), where the characteristic length scale and amplitude of the equilibrium
magnetic fields correspond to those of eddies: $a\rightarrow\lambda$ and
$\bar{B}_{\perp 0}\rightarrow\delta B_{\perp\lambda}$. From
$\gamma_{nl}\sim\gamma_{\rm max}$ we obtain

$\delta B_{\perp\lambda}^{2}\propto\delta
B_{\perp\lambda}^{-n\alpha/(n\alpha+1)}\lambda^{\alpha(2n+1)/(n\alpha+1)},$
(S17)

and so the magnetic energy spectrum should scale as

$E_{B}(k_{\perp})dk_{\perp}\propto
k_{\perp}^{-(7n\alpha+2\alpha+2)/(3n\alpha+2)}dk_{\perp}.$ (S18)

In the case of Laplacian resistivity, i.e., $\alpha=2$, the $n=1$ profile
gives the $k_{\perp}^{-5/2}=k_{\perp}^{-2.5}$ spectrum, and the $n=2$ profile
gives $k_{\perp}^{-17/7}\approx k_{\perp}^{-2.4}$. In the case of hyper-
resistivity with $\alpha=6$ (the order used in our isothermal simulations),
the $n=1$ profile gives the $k_{\perp}^{-14/5}=k_{\perp}^{-2.8}$ spectrum, and
the $n=2$ profile gives $k_{\perp}^{-49/19}\approx k_{\perp}^{-2.6}$. In Sec.
III.1 we obtained the corresponding scalings for the case where it is instead
electron inertia breaking the frozen flux constraint: the $n=1$ profile then
yields a $k_{\perp}^{-3}$ spectrum whereas $k_{\perp}^{-8/3}$ is obtained for
the $n=2$ case. Comparing the above derived spectra, we reach the conclusion
that in numerical simulations where the flux unfreezing is caused by different
mechanisms, if the spectra were indeed set by the tearing modes, the
difference of the measured spectra should be sufficient to identify.

###  .2 Kinetic Alfvén waves

The linear dispersion relation for kinetic Alfvén waves (KAWs) in the KREHM
framework is [zocco2011reduced]:

$\left[\zeta^{2}-\frac{\tau}{z_{c}}\frac{k_{\perp}^{2}d_{e}^{2}/2}{1-\Gamma_{0}(k_{\perp}^{2}\rho_{i}^{2}/2)}\right][1+\zeta
Z(\zeta)]=\frac{1}{2}k_{\perp}^{2}d_{e}^{2},$ (S19)

where $\zeta\equiv\Omega_{\rm KAW}/|k_{\parallel}|v_{{\rm th}e}$ ($\Omega_{\rm
KAW}=\omega_{\rm KAW}+i\gamma_{\rm KAW}$), $\tau\equiv T_{0i}/T_{0e}$ is the
temperature ratio, $z_{c}$ is the charge ratio,
$\Gamma_{0}(\alpha)=I_{0}(\alpha)e^{-\alpha}$ with $I_{0}$ being the zeroth-
order modified Bessel function, and
$Z(\zeta)\equiv\pi^{-1/2}\int_{-\infty}^{\infty}e^{-t^{2}}/(t-\zeta)dt$ is the
plasma dispersion function. In what follows, we specify $z_{c}=1$ and
$\tau=1$.

Fig. S1 shows the scale ($k_{\perp}\rho_{i}$) dependence of the oscillation
frequency, $\omega_{\rm KAW}$, and the (Landau) damping rate, $\gamma_{\rm
KAW}$, of KAWs obtained from Eq. (S19) for $\rho_{i}/d_{e}=40$. At around
$k_{\perp}d_{e}\sim 1$, $\omega_{\rm KAW}/k_{\parallel}v_{{\rm th}e}\sim 1$,
while at $k_{\perp}d_{e}\ll 1$, $\omega_{\rm KAW}\ll k_{\parallel}v_{{\rm
th}e}$. In this case, at $k_{\perp}d_{e}\ll 1$, the nonlinear eddy turn over
rate of the KAW turbulence, which is comparable to the KAW frequency
$\omega_{\rm KAW}$, is much smaller than the rate of linear phase mixing,
$k_{\parallel}v_{{\rm th}e}$. This timescale ordering leads to the zeroth-
order solution we describe in Sec. III.3.

As argued in Sec. III.3, the validity of our argument of strong electron
heating is predicated on the existence of a large enough dynamical range
$\rho_{i}\gg 1/k_{\perp}\gg d_{e}$ such that both condition for the zeroth-
order solution,

$\displaystyle\frac{\omega_{\rm KAW}}{k_{\parallel}v_{{\rm
th}e}}\sim\frac{k_{\parallel}V_{A}k_{\perp}\rho_{s}}{k_{\parallel}v_{{\rm
th}e}}\sim\frac{V_{A}}{v_{{\rm
th}e}}k_{\perp}\rho_{s}\sim\frac{d_{e}}{\rho_{i}}k_{\perp}\rho_{s}\sim
k_{\perp}d_{e}\ll 1,$ (S20)

and the considerable damping rate,

$\displaystyle-\frac{\gamma_{\rm KAW}}{\omega_{\rm KAW}}>\frac{1}{2\pi},$
(S21)

are satisfied.

Figure S1: Oscillation frequency ($\omega_{\rm KAW}$) and damping rate
($\gamma_{\rm KAW}$) of KAWs (normalized to $k_{\parallel}v_{{\rm th}e}$),
obtained from Eq. (S19) for $\rho_{i}/d_{e}=40$ (corresponding to run H2a ).
Vertical dotted line indicates the $d_{e}$ scale. Figure S2: Ratio of damping
rate ($\gamma_{\rm KAW}$) to frequency ($\omega_{\rm KAW}$) of KAWs for
varying $\rho_{i}/d_{e}$, obtained by solving Eq. (S19).

In Fig. S2 we plot the ratio of damping rate to frequency $\gamma_{\rm
KAW}/\omega_{\rm KAW}$ at different values of $\rho_{i}/d_{e}$. The horizontal
dashed line indicates $\gamma_{\rm KAW}/\omega_{\rm KAW}=1/2\pi$, i.e., the
value at which Landau damping of KAWs becomes significant; and the vertical
dotted lines show the corresponding values of $k_{\perp}\rho_{i}$ at which
this threshold is reached for each value of $\rho_{i}/d_{e}$. As
$\rho_{i}/d_{e}$ increases, there is a progressively wider range of scales
where damping of KAWs cannot be ignored.

###  .3 Multi-point, three-dimensional structure function

In order to test the prediction of the parallel spectrum and spectral
anisotropy (given by boldyrev2012spectrum and loureiro2017collisionless), and,
furthermore, to quantify the intermittency of the system, we calculate the
five-point, three-dimensional structure function of the total magnetic fields
$\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}},t)$. The reason to use a multi-
point structure function (instead of a more common two-point version [cho2000,
mallet2016]) is to capture the correct scaling for fluctuations whose energy
spectra are close to or steeper than $\sim k^{-3}$ [cho2009simulations,
cerri2019kinetic] and to filter out more effectively the large-scale
fluctuations with wavenumbers below $k\sim\pi/\delta r$, where $\delta r$ is
the incremental scale. Following cho2009simulations, we define the five-point,
$q$th-order structure function for magnetic-field fluctuations as

$S_{q}(\delta\mbox{\boldmath{$r$}})\equiv\langle|\Delta\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}},\delta\mbox{\boldmath{$r$}})|^{q}\rangle_{\mbox{\boldmath{$r$}}}.$
(S22)

Here $\delta\mbox{\boldmath{$r$}}$ is the increment vector,
$\langle...\rangle_{\mbox{\boldmath{$r$}}}$ is the ensemble average over the
domain, and the (normalized) local field increment is calculated with the
mentioned five-point stencil:

$\displaystyle\Delta\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}},\mbox{\boldmath{$\delta
r$}})$
$\displaystyle=[\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}+2\delta\mbox{\boldmath{$r$}})-4\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}+\delta\mbox{\boldmath{$r$}})+6\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}})-$
(S23) $\displaystyle
4\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}-\delta\mbox{\boldmath{$r$}})+\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}-2\delta\mbox{\boldmath{$r$}})]/\sqrt{35}.$

The structure function $S_{q}(\delta\mbox{\boldmath{$r$}})$ with respect to
$\delta\mbox{\boldmath{$r$}}$ can be projected into local coordinates with
respect to a local, scale-dependent coordinate system
$(\hat{\ell},\hat{\xi},\hat{\lambda})$, where $\hat{\ell}$ is parallel to the
local mean magnetic fields defined as

$\displaystyle\mbox{\boldmath{$B$}}_{\rm loc}$
$\displaystyle=[\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}+2\delta\mbox{\boldmath{$r$}})+4\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}+\delta\mbox{\boldmath{$r$}})+6\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}})+$
(S24) $\displaystyle
4\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}+\delta\mbox{\boldmath{$r$}})+\mbox{\boldmath{$B$}}(\mbox{\boldmath{$r$}}-2\delta\mbox{\boldmath{$r$}})]/16,$

$\hat{\xi}$ is the direction of the local perpendicular fluctuation
$\Delta\mbox{\boldmath{$B$}}_{\perp}=\mbox{\boldmath{$B$}}_{\rm
loc}\times(\Delta\mbox{\boldmath{$B$}}\times\mbox{\boldmath{$B$}}_{\rm loc})$,
and $\hat{\lambda}=\hat{\ell}\times\hat{\xi}$. We denote $\theta_{B}$ as the
angle between $\delta\mbox{\boldmath{$r$}}$ and $\mbox{\boldmath{$B$}}_{\rm
loc}$, and $\theta_{\Delta B_{\perp}}$ as that between $\hat{\xi}$ and the
projection of $\delta\mbox{\boldmath{$r$}}$ on the plane perpendicular to
$\mbox{\boldmath{$B$}}_{\rm loc}$. The structure function in the three
orthogonal directions can then be defined as

$S_{q}(\ell)\equiv
S_{q}(\ell;0^{\circ}<\theta_{B}<\delta\theta,0^{\circ}<\theta_{\Delta
B_{\perp}}<90^{\circ}),$ (S25) $S_{q}(\xi)\equiv
S_{q}(\xi;90^{\circ}-\delta\theta<\theta_{B}<90^{\circ},0^{\circ}<\theta_{\Delta
B_{\perp}}<\delta\theta),$ (S26) $S_{q}(\lambda)\equiv
S_{q}(\lambda;90^{\circ}-\delta\theta<\theta_{B}<90^{\circ},90^{\circ}-\delta\theta<\theta_{\Delta
B_{\perp}}<90^{\circ}),$ (S27)

where $\delta\theta$ is a finite angular threshold. When calculating the
directional structure functions, we reduce $\delta\theta$ until their scalings
appeared converged. Using the second-order structure function ($q=2$) and
equating the value between pairs of $S_{2}(\ell)$, $S_{2}(\xi)$, and
$S_{2}(\lambda)$, we obtain the (spectral anisotropy) scalings between $\ell$,
$\xi$, and $\lambda$.