# High-frequency heating of the solar wind triggered by low-frequency
turbulence

Jonathan Squire,1∗ Romain Meyrand,1 Matthew W. Kunz,2,3
Lev Arzamasskiy,4 Alexander A. Schekochihin,5,6 Eliot Quataert2

1Physics Department, University of Otago, Dunedin 9010, New Zealand,
2Department of Astrophysical Sciences, Princeton University, Peyton Hall,
Princeton, NJ 08544, USA
3 Princeton Plasma Physics Laboratory, PO Box 451, Princeton, NJ 08543, USA
4 School of Natural Sciences, Institute for Advanced Study, Princeton, NJ
08544, USA
5 The Rudolf Peierls Centre for Theoretical Physics, University of Oxford,
Clarendon Laboratory, Parks Road,
Oxford, OX1 3PU, UK
6 Merton College, Oxford OX1 4JD, UK
∗To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>
> The fast solar wind’s high speeds and nonthermal features require that
> significant heating occurs well above the Sun’s surface. Two leading
> theories seem incompatible: low-frequency “Alfvénic” turbulence, which
> transports energy outwards and is observed ubiquitously by spacecraft but
> struggles to explain the observed dominance of ion over electron heating;
> and high-frequency ion-cyclotron waves (ICWs), which explain the nonthermal
> heating of ions but lack an obvious source. Here, we argue that the recently
> proposed “helicity barrier” effect, which limits electron heating by
> inhibiting the turbulent cascade of energy to the smallest scales, can unify
> these two paradigms. Our six-dimensional simulations show how the helicity
> barrier causes the large-scale energy to grow in time, generating small
> parallel scales and high-frequency ICW heating from low-frequency
> turbulence, while simultaneously explaining various other long-standing
> observational puzzles. The predicted causal link between plasma expansion
> and the ion-to-electron heating ratio suggests that the helicity barrier
> could contribute to key observed differences between fast- and slow-wind
> streams.

The basic mechanisms that heat the solar corona and accelerate the solar wind
remain mysterious despite intensive study over many decades1. A successful
theory must explain how energy contained in photospheric motions and magnetic
fields can be liberated to cause extreme and sudden heating of the coronal
plasma, along with its acceleration to velocities well in excess of the escape
velocity of the Sun. Adding to the complexity, the coronal plasma is
collisionless — the mean-free path of protons can be large compared to the
largest observed structures — meaning it can be far out of local thermal
equilibrium. This freedom opens up a wide array of channels for plasma heating
— ions might be heated more than electrons (or vice versa), or particles might
gain energy preferentially in a particular direction with respect to the local
magnetic field 2. Such differences can have significant macroscopic
consequences.

The dominant heating mechanism(s) must be consistent with an extensive array
of measurements taken both remotely, from the low corona itself, and _in situ_
, from spacecraft spread throughout the solar wind. In fast-wind streams,
these data indicate that the heating must be spatially extended out to several
solar radii in order to drive observed wind speeds 3. It must preferentially
heat protons over electrons 4, while heating heavier ions (e.g., alpha
particles) even more effectively 5. It must heat protons in the direction
perpendicular to the local magnetic field significantly more than in the
parallel direction, in order to explain temperature anisotropies 6. And, its
features and/or aftereffects should be observable in the measured field
fluctuations and particle distribution, particularly at the low altitudes now
being explored by Parker Solar Probe (PSP) 7.

A paradigm that can — at least in principle — satisfy the above requirements
is heating through “Alfvénic” turbulence. Low-frequency Alfvénic motions in
the low corona are observed to contain sufficient energy to power the wind 8,
9, and there are well-developed theories for how such motions become turbulent
following reflection from large-scale density gradients 10, 11, 12. This
turbulence transfers energy into successively smaller-scale motions
perpendicular to the magnetic field (larger $k_{\perp}$, where $k_{\perp}$ is
the inverse perpendicular scale), ultimately dissipating to heat the plasma.
The difficulty is that most theories predict that, in the strongly magnetized
limit relevant to the solar corona (the low-$\beta$ limit, where $\beta$ is
the ratio of thermal to magnetic pressure), such low-frequency,
high-$k_{\perp}$ structures dissipate to heat predominantly _electrons_ 13,
14. Other low-frequency plasma motions, such as compressive waves, generally
cause parallel heating of ions 15. Both possibilities are inconsistent with
observations. More promisingly, for turbulence of sufficient amplitudes,
“stochastic heating” 16 can heat ions through a random walk on ion-gyroscale
electric-field fluctuations. Although it can plausibly explain key
observations17, 18, questions remain, such as its possible quenching due to
flattening of the distribution function19, 20 and the influence of cross
helicity21. Another possibility — that ions are heated by kinetic-Alfvén-wave
(KAW) turbulence at sub-gyroradius scales 19, 22 — remains less well
understood and may be inefficient at low $\beta$ 20.

In the opposite limit of short field-parallel wavelengths (large $k_{\|}$),
high-frequency ion-cyclotron waves (ICWs) provide a simpler mechanism to cause
strong perpendicular ion heating 6, 23. At wavenumbers approaching $k_{\|}\sim
d_{i}^{-1}$, where $d_{i}$ is the ion inertial length, their frequency
approaches the ion gyrofrequency, where the cyclotron resonance causes highly
efficient energy transfer from electromagnetic fields to ion velocities 24,
25. ICWs are observed ubiquitously _in situ_ 7 and can suprathermally heat
minor ions in a way that is observationally compelling 26, 27. However, a
sufficiently energetic direct solar source of ICWs is highly unlikely28, and
the Alfvénic cascade does not efficiently transfer energy to small parallel
scales 29, seemingly ruling out their turbulent origin. Although their
occurrence in data can be explained by kinetic instabilities30, in most
theories this implies they would cool, rather than heat, the plasma24, except
perhaps in the presence of strong nonthermal particle beams31.

If combined, these two heating paradigms — via Alfvénic turbulence or ICWs —
can conceivably satisfy the fast-wind heating requirements described above,
maintaining an abundant source of perpendicular ion heating well above the
solar surface. Here we assess whether a newly discovered effect, termed the
“helicity barrier”32, can fulfil this role by obstructing the dissipation of
collisionless Alfvénic turbulence into electron heat. Using six-dimensional,
high-resolution, hybrid-kinetic simulations, we explore the effect of the
helicity barrier on collisionless turbulent heating, choosing parameters to
match as closely as possible the conditions observed in fast-wind streams. We
assess the relevance of our results to the solar wind by comparing detailed
features of the turbulent spectra and ion distribution function to
observations from PSP and other spacecraft.

##### The helicity barrier

Solar-wind turbulence is _imbalanced_ (possessing cross helicity), meaning
that it is energetically dominated by Alfvénic structures that propagate
outward from the Sun (designated $z^{+}$; the inward-propagating component is
designated $z^{-}$). The theory of highly perpendicular ($k_{\perp}\gg
k_{\|}$) perturbations in a collisionless plasma predicts that, at scales
larger than the ion gyroradius $\rho_{i}$, imbalanced turbulent Alfvénic
energy in $z^{\pm}$ can cascade towards smaller scales (larger $k_{\perp}$),
as required to heat the plasma 15. In constrast, at sub-gyroradius scales
$k_{\perp}\rho_{i}\gtrsim 1$, magnetic helicity conservation implies the
opposite — imbalanced energy must cascade _inversely_ , towards larger scales
15, 33, 32. The effect is directly analogous to 2D hydrodynamic turbulence,
where enstrophy conservation causes energy to cascade inversely, creating
large-scale vortex structures. However, unlike in hydrodynamics, in
low-$\beta$ plasmas, cross helicity at $k_{\perp}\rho_{i}\lesssim 1$
transforms conservatively into magnetic helicity at $k_{\perp}\rho_{i}\gtrsim
1$ (the system conserves a generalized helicity 32). This implies that an
imbalanced energy flux arriving at $k_{\perp}\rho_{i}\sim 1$ from large scales
cannot cascade to arbitrarily small scales.

Mathematically, it is helpful to separate the turbulent energy flux
$\varepsilon$ into components associated with the large-scale outward
($\varepsilon^{+}$) and inward ($\varepsilon^{-}$) propagating fluctuations:
$\varepsilon=\varepsilon^{+}+\varepsilon^{-}$. Generalized-helicity
conservation prohibits the conversion of $\varepsilon^{+}$ into
$\varepsilon^{-}$ at any scale (note, however, that at
$k_{\perp}\rho_{i}\gtrsim 1$, $\varepsilon^{+}$ is associated with a mixture
of outward- and inward- propagating KAWs). Because a forward cascade must be
balanced ($\varepsilon^{+}\approx\varepsilon^{-}$) at
$k_{\perp}\rho_{i}\gtrsim 1$, only a small portion ${\sim}2\varepsilon^{-}$ of
the energy flux can cascade to small scales where it will heat electrons 14.
The rest of the flux, ${\sim}(\varepsilon^{+}-\varepsilon^{-})$, is stuck — it
hits the “helicity barrier” and thus remains at scales
$k_{\perp}\rho_{i}\lesssim 1$. If the system is forced continuously, this
large-scale energy grows in time with a decreasing parallel correlation length
32, as expected from critical balance 15. We show that this growth eventually
funnels the turbulent energy into a spectrum of ICW fluctuations, heating the
ions, which absorb the majority of the energy flux.

### Numerical method

Our simulation uses the Pegasus++ code 34, which solves the hybrid-kinetic
equations with isothermal electrons using the particle-in-cell method. The
system is strongly magnetized with mean magnetic field
$\bm{B}_{0}=-B_{0}\hat{\bm{z}}$, Alfvén speed $v_{\rm A}=B_{0}/\sqrt{4\pi
nm_{i}}$, and initial ion $\beta$, $\beta_{i}\equiv 8\pi
nk_{B}T_{i}/B_{0}^{2}=0.3$ ($m_{i}$, $n$, and $k_{B}T_{i}=m_{i}v_{\rm
th}^{2}/2$ are the ion mass, number density, and temperature, respectively
with $v_{\rm th}$ the ion thermal velocity and $k_{B}$ the Boltzmann
constant). Perpendicular ($x$ and $y$ directed) ion-velocity fluctuations
$\bm{u}_{\perp}$ and magnetic fluctuations $\bm{B}_{\perp}$ are driven at
large scales and correlated to create imbalance, with
$\bm{z}^{+}\equiv\bm{u}_{\perp}+\bm{B}_{\perp}/\sqrt{4\pi
nm_{i}}\gg\bm{z}^{-}\equiv\bm{u}_{\perp}-\bm{B}_{\perp}/\sqrt{4\pi nm_{i}}$,
($\bm{z}^{+}$ perturbations propagate in the $+\hat{\bm{z}}$ direction). The
energy-injection rate $\varepsilon$ and cross-helicity injection rate
$\varepsilon_{H}=\varepsilon^{+}-\varepsilon^{-}=0.9\varepsilon$ are constant
in time. Plasma heating is strongly influenced by the amplitude $\delta
B_{\perp}/B_{0}$ and spectral anisotropy $k_{\|}/k_{\perp}$ of fluctuations
with $k_{\perp}\rho_{i}\sim 1$. In order to reach realistic values without
simulating the enormous scale separation of the real solar wind, we use a
highly elongated box with dimensions $L_{\|}=L_{z}=6L_{\perp}$; given the box
size $k_{\perp 0}\rho_{i}\equiv 2\pi\rho_{i}/L_{\perp}\approx 0.05$, this
gives conditions near $k_{\perp}\rho_{i}\sim 1$ that are comparable to those
observed 35. The elongated domain also implies that the timescales we probe
are rapid compared to the solar wind’s outer scales and its expansion rate,
justifying the external forcing to represent a turbulent flux from larger
scales36 and our neglect of expansion effects. The simulation’s resolution is
$N_{\perp}^{2}\times N_{\|}=392^{2}\times 2352$ cells, so that the smallest
resolved scales have $k_{\perp,{\rm max}}\rho_{i}\simeq\pi
N_{\perp}\rho_{i}/L_{\perp}=10$. Most other simulation parameters, including a
(hyper-)resistivity that dissipates small-scale magnetic energy, are chosen to
match a previous balanced turbulence simulation with $\beta_{i}=0.3$ 19. This
allows direct comparison of their spectra and heating. Further details are
provided in Methods.

### Results

Figure 1: The time evolution confirms the formation of the helicity barrier.
(A) Outward/inward-propagating fluctuation amplitudes ($z^{\pm}_{\rm rms}$)
and imbalance ($\sigma_{c}$); (B) Energy budget and heating, illustrated
through the ion heating rate ($Q_{i}$) and the small-scale resistive
dissipation ($\varepsilon_{\eta}$). The balanced part of the injected energy
(${\sim}2\varepsilon^{-}$) saturates early by $t\simeq 3\tau_{\rm
A}=3L_{\|}/v_{\rm A}$ (blue shaded region); this is demonstrated by the
saturation of $z^{-}$ (blue line in A) and the resistive dissipation (purple
line in B), which absorbs only a small fraction of the input energy
$\varepsilon_{\eta}\approx 2\varepsilon^{-}\ll\varepsilon$. Eventually, by
$t\approx 14\tau_{\rm A}$ (orange shaded region), ion heating absorbs the
remainder of the injected energy so that
$Q_{i}\approx\varepsilon-\varepsilon_{\eta}$ (brown line in B), halting the
growth of $z^{+}$. A numerical cooling effect has been removed to compute the
energy budget (see Methods and Extended Data Fig. 1).

The simulation’s time evolution, shown in Fig. 1, exhibits several features
that are expected from the helicity barrier 32 but not from other theories of
imbalanced turbulence 37, 32. Panel A shows the growth of the root-mean-square
ampliutudes $z^{\pm}_{\rm rms}/v_{\rm
A}\equiv\langle(z^{\pm})^{2}\rangle^{1/2}/v_{\rm A}$ and imbalance (normalized
cross helicity) $\sigma_{c}=2\langle\sqrt{4\pi
nm_{i}}\bm{u}\bm{\cdot}\bm{B}\rangle/\langle 4\pi
nm_{i}\bm{u}^{2}+\bm{B}^{2}\rangle$ (where $\langle\,\dots\,\rangle$ denotes a
box average). Because $z^{+}_{\rm rms}/v_{\rm A}\approx 2u_{\perp,{\rm
rms}}/v_{\rm A}\approx 2B_{\perp,{\rm rms}}/B_{0}$, the final $\delta
B_{\perp}/B_{0}\approx 0.26$ is nearly twice that of the balanced simulation
($\delta B_{\perp}/B_{0}\approx 0.14$). We see that $z^{-}$ saturates quickly,
by time $t\simeq 3\tau_{\rm A}$, while saturation of $z^{+}$ occurs only after
$t\simeq 14\tau_{\rm A}$ ($\tau_{\rm A}\equiv L_{\|}/v_{\rm A}$ is the Alfvén
time). This is expected: energy in $z^{-}$ can cascade to
$k_{\perp}\rho_{i}\gtrsim 1$ through standard KAW turbulence, while most of
the energy in $z^{+}$ cannot cascade past $k_{\perp}\rho_{i}\sim 1$ due to the
helicity barrier. The imbalance saturates at $\sigma_{c}\approx 0.98$, which,
though large, is regularly observed by PSP 38. Panel B shows the ion heating
rate $Q_{i}=V\partial_{t}\langle 3m_{i}v_{\rm th}^{2}/2\rangle$ (where $V$ is
the simulation volume) and the (hyper-)resistive dissipation rate
$\varepsilon_{\eta}$. The latter is a proxy for electron heating because it
absorbs the energy that cascades to the smallest scales below where ions can
respond to the fluctuations. We see that $\varepsilon_{\eta}$ saturates
rapidly, together with $z^{-}$ and well before the total energy. Its value,
$\varepsilon_{\eta}\approx 2\varepsilon^{-}$ approximates the balanced portion
of the injected flux, as expected. The remainder of the energy input,
$\varepsilon-\varepsilon_{\eta}$, must eventually go into ion heating $Q_{i}$.
However, in order for this to happen, $z^{+}$ must grow significantly, which,
notably, occurs without change to $\varepsilon_{\eta}$ or $z^{-}$. The
simulation saturates with $Q_{i}\approx\varepsilon-\varepsilon_{\eta}$ at
$t\approx 14\tau_{\rm A}$; we run it for another ${\approx}4.2\tau_{\rm A}$ in
this steady state.

Figure 2: Fluctuation spectra exhibit a steep transition region around
$k_{\perp}\rho_{i}=1$. Spectra of the magnetic field
$\mathcal{E}_{\bm{B}}=\mathcal{E}_{B_{x}}+\mathcal{E}_{B_{y}}+\mathcal{E}_{B_{z}}$,
electric field
$\mathcal{E}_{\bm{E}}=\mathcal{E}_{E_{x}}+\mathcal{E}_{E_{y}}+\mathcal{E}_{E_{z}}$,
and KAW-normalized density $n_{\rm KAW}=\sqrt{\beta_{i}(1+2\beta_{i})}n$,
binned in $k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$. The lower panel shows the
local scaling exponents, $\alpha={\rm d}\ln\mathcal{E}/{\rm d}\ln k_{\perp}$;
at small scales (approximately the shaded region) spectra are affected by
particle noise (see Extended Data Fig. 2). The inset shows the extremely steep
$k_{\|}$ spectrum of $B_{\perp}$
($\mathcal{E}_{B_{\perp}}=\mathcal{E}_{B_{x}}+\mathcal{E}_{B_{y}}$) near
$k_{\|}d_{i}=1$. Thin black lines show $\mathcal{E}_{\bm{B}}$ for an analogous
balanced-turbulence simulation 19.

Perpendicular magnetic ($\mathcal{E}_{\bm{B}}$), electric
($\mathcal{E}_{\bm{E}}$), and density ($\mathcal{E}_{n_{\rm KAW}}$)
fluctuation spectra are shown in Fig. 2. The signature of the helicity barrier
is the sharp spectral break at $k_{\perp}\rho_{i}<1$ in $\mathcal{E}_{\bm{B}}$
and $\mathcal{E}_{\bm{E}}$, which migrates to larger scales with time 32.
“Double-kinked” spectra — a ${\sim}k_{\perp}^{-5/3}$ range at large scales, a
steep ${\sim}k_{\perp}^{-4}$ transition range bracketing
$k_{\perp}\rho_{i}=1$, and a flatter range of KAW turbulence
(${\sim}k_{\perp}^{-2.8}$) at yet smaller scales — have been observed for
decades39 including by PSP close to the Sun40, but have lacked a clear
theoretical explanation. So far as we are aware, this is the first self-
consistent kinetic simulation to exhibit this feature. Although
$\mathcal{E}_{\bm{B}}$ in Fig. 2 lacks a ${\sim}k_{\perp}^{-2.8}$ range due to
box resolution (resistivity), it clearly re-flattens at
$k_{\perp}\rho_{i}\approx 2$, and the width of the transition range is
comparable to observations 40. The $\mathcal{E}_{\bm{E}}$ spectrum has a
similar shape, flattening after a steep drop around $k_{\perp}\rho_{i}\sim 1$.
The spectrum of $n_{\rm KAW}=\sqrt{\beta_{i}(1+2\beta_{i})}n$ is predicted 15
and observed 41 to satisfy $\mathcal{E}_{n_{\rm
KAW}}\approx\mathcal{E}_{\bm{B}}$ in KAW turbulence because linear KAWs
satisfy $\delta B_{\perp}/B_{0}\approx(\delta
n/n)\sqrt{\beta_{i}(1+\beta_{i})}\approx(\delta
B_{\|}/B_{0})\sqrt{1+1/\beta_{i}}$ and $\delta B^{2}=\delta
B_{\perp}^{2}+\delta B_{\|}^{2}$ — this occurs here for
$k_{\perp}\rho_{i}\gtrsim 1$, providing further evidence for the KAW-like
nature of the sub-$\rho_{i}$ turbulence. The parallel spectrum (see Fig. 2
inset and Methods) exhibits a ${\sim}k_{\|}^{-2}$ range at large scales,
followed by a very steep ${\sim}k_{\|}^{-6}$ range bracketing $k_{\|}d_{i}\sim
1$, which flattens at smaller scales; again these features match recent PSP
observations 42. Balanced turbulence19, unlike the imbalanced case, does not
exhibit a steep transition range in either $k_{\perp}$ or $k_{\|}$, which is
likely also the case in the solar wind43. Note that spectra are adversely
affected by particle noise for $k_{\perp}\rho_{i}\gtrsim 3$ (see Methods and
Extended Data figures 1-2).

Figure 3: Evidence for ICW fluctuations in the saturated state. (A) Structure
of $E_{\perp}$; superimposed white lines show magnetic-field lines projected
onto the plane. (B) Perpendicular magnetic-fluctuation spectrum
($\mathcal{E}_{B_{\perp}}=\mathcal{E}_{B_{x}}+\mathcal{E}_{B_{y}}$) in
($k_{\perp},\,k_{\|}$), illustrating the path of energy to $k_{\|}d_{i}\sim 1$
scales (labelled “Energy flow”) and the parallel ICW bump at
$k_{\|}d_{i}\simeq 0.8$. The angled white dotted lines show the canonical
critical balance result15 based on $z^{+}_{\rm rms}$ (top) or the box outer
scale (bottom). (C) Normalized magnetic helicity spectrum
$k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}$(see text) in
($k_{\perp},\,k_{z}$), showing the reversed helicity signatures of KAWs
($k_{\perp}\gg k_{\|}$) and ICWs ($k_{\|}\gg k_{\perp}$). For perpendicular
structures, $\mathcal{H}_{\bm{k}}\approx 0$ at both large scales
($k_{\perp}\rho_{i}\ll 1$), because large-scale Alfvén waves are nonhelical,
and in the KAW range ($k_{\perp}\rho_{i}\gg 1$), because the KAW turbulence
must be balanced. In between, in the transition range around
$k_{\perp}\rho_{i}\simeq 1$, $\mathcal{H}_{\bm{k}}$ has a local maximum
because the Alfvénic fluctuations are both dispersive and imbalanced.

Evidence for the presence of ICWs in the saturated state is shown in Fig. 3.
The projected magnetic-field lines and electric field (panel A) reveal the
coexistence of $k_{\|}d_{i}\sim 1$ parallel structure with the sub-$\rho_{i}$
striations of KAW turbulence in the perpendicular plane. More quantitatively,
panel B shows the two-dimensional $(k_{\perp},k_{\|})$ spectrum of $B_{\perp}$
(see Methods and Extended Data figure 3). Several features are manifest:
first, unlike in the analogous balanced-turbulence simulation 19, the outer-
scale parallel correlation length is significantly smaller than $L_{\|}$
because it decreases with increasing amplitude to maintain critical balance,
$k_{\|}v_{\rm A}\sim k_{\perp}z^{+}$ (the region of maximal spectral power
moves upwards in time as $z^{+}$ grows); second, the cone of maximal spectral
power appears to steepen as the turbulence moves to smaller scales, creating
small parallel scales faster than the canonical result15 $k_{\|}\propto
k_{\perp}^{2/3}$ (or $k_{\|}\propto k_{\perp}^{1/2}$ for aligned
turbulence37); third, there is a clear spectral bump at $k_{\|}d_{i}\simeq
0.8$ and $k_{\perp}<k_{\|}$, the signature of ICWs. By integrating the energy
spectrum over modes with $k_{\perp}\leq k_{\|}$, we estimate that these ICW
modes contain ${\simeq}1\%$ of the total energy; this fraction grows by a
factor ${\simeq}50$ from $t\simeq 5\tau_{\rm A}$ to saturation and exceeds
that of saturated balanced turbulence by a factor ${\simeq}30$. We confirm
that these modes are indeed ICWs in panel C, which shows the normalized
magnetic-helicity spectrum $k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}$, where
$\mathcal{H}_{\bm{k}}=i(B_{x}^{*}B_{y}-B_{y}^{*}B_{x})/k_{z}$, binned in
($k_{\perp},\,k_{z}$). As predicted 44 and observed 45, ICWs (with
$k_{\perp}\lesssim k_{\|}$) are characterised by
$k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}<0$, while
$k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}>0$ for $k_{\perp}\gtrsim k_{\|}$
near $k_{\perp}\rho_{i}\simeq 1$. This results from the intrinsic polarization
of ICWs and their propagation direction;
$k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}\approx-1$ indicates that ICWs
propagate almost exclusively in the $+\bm{\hat{z}}$ direction like the large-
scale $z^{+}$. At $k_{\perp}>k_{\|}$,
$k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}>0$ is a result of the dominantly
Alfvénic perturbations becoming dispersive at $k_{\perp}\rho_{i}\sim 1$;
helicity grows to $k\mathcal{H}_{\bm{k}}/\mathcal{E}_{\bm{B}}\approx 0.5$ at
$k_{\perp}\rho_{i}\simeq 1$ before decreasing again because the small-scale
KAW cascade is balanced. This feature, which is a theoretical corollary of the
helicity barrier, has been commonly observed by PSP and other spacecraft 45
and correlates with the transition-range spectral slope as predicted46.

Evidence that plasma heating occurs through ICWs is provided in Fig. 4. Quasi-
linear cyclotron-heating theory 24 is based on the idea that ions and ICWs
interact strongly if the wave frequency $\omega_{\bm{k}}$ is resonant with the
Doppler-shifted ion gyromotion. When ICWs exist across a range of $k_{\|}$,
the process flattens the ion distribution function $f_{i}(w_{\perp},w_{\|})$
along specific “scattering contours,” which can be computed25 from
$\omega_{\bm{k}}$ for waves of a particular $k_{\perp}$ ($w_{\perp}$ and
$w_{\|}$ are the field-perpendicular and -parallel velocities of ions in the
frame moving with the plasma). Theory suggests that, because the scattering
contours steepen with increasing $k_{\perp}/k_{\|}$, heating by oblique ICWs
generates an $f_{i}$ that increases along the scattering contours of parallel
ICWs 47. The consequences are twofold: first, oblique-ICW heating generates
parallel ICWs, explaining the dominance of $k_{\perp}\ll k_{\|}$ modes in Fig.
3B; second, in quasi-steady state, $f_{i}$ is nearly flat along the _parallel_
ICW scattering contours. We plot these, along with $f_{i}(w_{\perp},w_{\|})$,
in Fig. 4A. There is exceptionally good agreement at saturation for
$w_{\|}\lesssim-w_{\|,d_{i}}$ particles, which are those that can resonate
with $k_{\|}d_{i}\lesssim 1$ ICWs propagating in the $+\hat{\bm{z}}$ direction
(see Methods). The time evolution is also telling: the quasi-linear flattening
starts at large $|w_{\|}|$ and moves upwards as time advances, which is
expected because there is more power in low-$k_{\|}$ modes that resonate with
high-$|w_{\|}|$ ions. Panel B shows the perpendicular energy diffusion
coefficient $D^{E}_{\perp\perp}$, computed for $w_{\|}<-w_{\|,d_{i}}$ using
the time evolution of $f_{i}$48 and validated by computing
$\langle\bm{E}_{\perp}\bm{\cdot}\bm{w}_{\perp}\rangle$ directly from particle
trajectories 20. Quasi-linear ICW theory predicts24 $D^{E}_{\perp\perp}\propto
w_{\perp}^{2}$, as seen in Fig. 4, while a stochastically heated plasma has20
$D^{E}_{\perp\perp}\sim{\rm const}$ for $w_{\perp}\sim v_{\rm th}$. This
quantitatively confirms the dominance of quasi-linear ICW heating and we find
no other evidence for stochastic heating in this simulation, although it is
possible that it could govern saturation under different conditions (e.g.,
lower $\beta$; see Methods and Extended Data figure 4). Finally, we see a
clear flattening of $f_{i}(w_{\perp},w_{\|})$ at $w_{\|}\approx v_{\rm A}$ and
small $w_{\perp}$, which forms a modestly super-Alfvénic beam feature in the
direction of dominant wave propagation, similar to those observed in the fast
solar wind 49, 2, 50. By comparing $\partial f_{i}/\partial t$ with
$\langle{E}_{\|}{w}_{\|}\rangle$ (see Extended Data figure 5), we confirm that
this arises through Landau damping of Alfvén waves as their phase velocity
increases near $k_{\perp}\rho_{i}\sim 1$ 51.

Figure 4: The evolution of the ion distribution function
$f_{i}(w_{\perp},w_{\|})$ shows how ICWs heat ions. (A) Structure of $f_{i}$
at different times. At saturation, $f_{i}$ aligns nearly perfectly with the
scattering contours of parallel ICWs (dashed lines) for $w_{\|}\lesssim-
w_{\|,d_{i}}$ (red line), which is the region of velocity space where there
exists resonant ICWs. (B) Perpendicular energy diffusion coefficient
$D^{E}_{\perp\perp}$; solid lines include only $w_{\|}<-w_{\|,d_{i}}$; the
dashed line includes all $w_{\|}$, with the lower values of
$D^{E}_{\perp\perp}$ indicating less heating for $w_{\|}>0$.

### Discussion

Together, Figs. 3 and 4 provide strong evidence for heating through ICWs
generated by imbalanced Alfvénic turbulence regulated by the helicity barrier.
Given the excellent match to _in-situ_ observations of spectra, helicity, and
distribution functions, we suggest the same is true in the fast solar wind,
reconciling the paradigms of low-frequency Alfvénic turbulence and ICW
heating. If, as suggested above, ICWs are produced due to the decrease in
outer parallel scale with increasing amplitude, then the turbulence will
saturate once it has grown sufficiently for $k_{\|}d_{i}\sim 1$ scales to be
reached (where Alfvén waves become ICWs) before the $k_{\perp}\rho_{i}\sim 1$
scales (where Alfvén waves become KAWs). If large-scale fluctuations are
critically balanced with spectrum $\mathcal{E}_{\bm{B}}\propto
k_{\perp}^{-5/3}$, this suggests a critical outer-scale (${\sim}L_{\perp}$)
saturation amplitude for appreciable ICW heating given by $\delta z^{+}/v_{\rm
A}\sim A\beta_{i}^{1/2}(\rho_{i}/L_{\perp})^{-1/3}$, with our simulation
giving a proportionality factor $A\simeq 0.2$ (an observationally testable
prediction). A corollary is that the energy-injection rate sets only the
timescale to reach saturation, not the amplitude, unlike standard viscous
damping. Correspondingly — and in contrast to other possible ion-heating
mechanisms — if the fluctuations’ amplitude is too small at some time (or in
some region of a radially stratified wind) then energy is not deposited into
electrons. Rather, the helicity barrier halts the cascade, storing the energy
in fluctuations that keep growing, eventually to heat ions through ICWs. This
mechanism is expected to become more robust with decreasing $\beta$, at least
for $\beta\gg m_{e}/m_{i}$ (for $\beta\sim m_{e}/m_{i}$, electron heating
effects are important at $k_{\perp}\rho_{i}\sim 1$). Our finding that the
helicity barrier occurs even at modest $\beta$ ($\beta\approx 0.3$) indicates
that it should apply nearly everywhere in the corona. Thus, the ratio of
electron to ion heating — an important input to larger-scale models of the
solar wind and other astrophysical objects — is simply
$Q_{i}/Q_{e}\approx\varepsilon_{H}/(\varepsilon-\varepsilon_{H})$ for
saturated imbalanced Alfvénic turbulence, nearly independent of $\beta$ for
$m_{e}/m_{i}\ll\beta\lesssim 1$.

A heating ratio that is controlled by imbalance is appealing for understanding
the generation of fast- and slow-wind streams by reflection-driven Alfvénic
turbulence. Fast wind — which emerges from coronal holes with lower expansion
factors 52 and thus less wave reflection and larger
$\varepsilon_{H}/\varepsilon$ — is observed to have hotter protons, strong
minor-ion heating, cooler electrons, larger imbalance, and steep transition-
range spectra. Slow wind — which arises from less ordered fields involving
closed loops and/or larger expansion factors, and thus has smaller
$\varepsilon_{H}/\varepsilon$ (or even $\varepsilon_{H}/\varepsilon\ll 1$ in a
closed-field region) — is observed to have cooler protons, little suprathermal
minor-ion heating, hotter electrons, less imbalance, and does not usually
exhibit a steep kinetic transition range. While reflection-driven turbulence
models can already reproduce the observed correlation of expansion factor with
wind speed based on the radial location of energy deposition53, 54, the
additional physics afforded by the helicity barrier yields interesting
implications. First, it would naturally explain the other aforementioned
correlations of electron and ion thermodynamics with wind speed. Second,
because of electrons’ high thermal speeds, a given quantity of energy
deposited into ion heat generally drives a higher asymptotic wind velocity
than if electrons are heated4. We propose that the helicity barrier could act
as a “switch” effect, supplementing other turbulent heating physics: low-
expansion regions with a robust barrier would heat predominantly ions at large
radii, ideal conditions for generating fast wind speeds; high-expansion
factors would break the barrier, depositing energy into electron heat and
exacerbating the inefficient acceleration that results from heating at lower
radii. Thus, by linking plasma thermodynamics to magnetic-field expansion, it
is plausible that the helicity barrier plays an important role in generating
the bimodal speed distribution of the solar wind.

## Methods

### Hybrid-kinetic simulation method

The equations of the hybrid-kinetic model solved by Pegasus++ are 55, 34,

$\displaystyle\frac{\partial f_{i}}{\partial
t}+\bm{v}\cdot\mbox{\boldmath{$\nabla$}}f_{i}+\left[\frac{e}{m_{i}}\left(\bm{E}+\frac{\bm{v}}{c}\,\mbox{\boldmath{$\times$}}\,\bm{B}\right)+\frac{\bm{F}_{\perp}^{U}}{m_{i}}\right]\cdot\frac{\partial
f_{i}}{\partial\bm{v}}=0,$ (1a) $\displaystyle\frac{\partial\bm{B}}{\partial
t}=-c\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,(\bm{E}+\bm{F}_{\perp}^{B})+\eta_{4}\mbox{\boldmath{$\nabla$}}^{4}\bm{B},$
(1b)
$\displaystyle\bm{E}=-\frac{\bm{u}_{{i}}\,\mbox{\boldmath{$\times$}}\,\bm{B}}{c}+\frac{({\mbox{\boldmath{$\nabla$}}}\,\mbox{\boldmath{$\times$}}\,\bm{B})\,\mbox{\boldmath{$\times$}}\,\bm{B}}{4\pi
en_{{i}}}-\frac{T_{e}\mbox{\boldmath{$\nabla$}}n_{{i}}}{en_{{i}}}.$ (1c)

The kinetic equation (1a) is solved in a six-dimensional $(\bm{x},\bm{v})$
space using a particle-in-cell method, while Faraday’s law (1b) and the
kinetic Ohm’s law (1c) are solved on a three-dimensional grid. A single ion
species of charge $e$ and mass $m_{i}$ is assumed, $c$ is the speed of light,
and $\bm{F}_{\perp}^{U}$ and $\bm{F}_{\perp}^{B}$ are the forcing terms, which
are described belows. In equation (1c), $n_{i}=\int{\rm d}\bm{v}f_{i}$ and
$\bm{u}_{i}=n_{i}^{-1}\int{\rm d}\bm{v}\,\bm{v}f_{i}$ are computed from
$f_{i}$. Equations (1a)–(1c) are derived from the two-species Vlasov equation
by expanding the electron equation in $m_{e}/m_{i}\ll 1$, assuming quasi-
neutrality $n_{e}=n_{i}=n$ and isothermal electrons (temperature $T_{e}$).
Pegasus++ uses a second-order accurate predictor-predictor-corrector scheme to
enforce the kinetic Ohm’s law 34 and has been highly optimized for efficient
operation on large supercomputing systems. The hyper-resistivity in equation
(1b) is not intended to represent a true physical resistivity (it does not
contribute to the electric field in the particle push), but is included to
dissipate energy in the magnetic field at the smallest grid scales. This is a
proxy for electron heating in the model.

#### Simulation parameters

As discussed in the main text, the simulation domain is elongated by a factor
of $6$ ($L_{x}=L_{y}=67.5d_{i}$ and $L_{z}=404.7d_{i}$, but with cubic grid
cells, $N_{x}=N_{y}=392$, $N_{z}=2352$) in order to realize realistic solar-
wind conditions near $k_{\perp}\rho_{i}\sim 1$. These conditions can be
estimated roughly by taking the outer-scale fluctuations to be in approximate
critical balance $\delta B_{\perp}/B_{0}\simeq L_{\perp}/L_{\parallel}$, and
assuming a magnetic spectrum $\mathcal{E}_{B}(k_{\perp})\propto
k_{\perp}^{-5/3}$ and $k_{\|}\propto k_{\perp}^{2/3}$ in the MHD inertial
range 56, 15. Given the outer perpendicular scale $k_{\perp 0}\rho_{i}\equiv
2\pi\rho_{i}/L_{\perp}\approx 0.05$, this suggests that at
$k_{\perp}\rho_{i}\sim 1$ the spectral anisotropy is
$k_{\perp}/k_{\|}\simeq\tan 86.5^{\circ}$ with fluctuation amplitude $\delta
B_{\perp}/B_{0}\simeq 0.06$. This is comparable to observed values (see Ref.
35, figure 1) justifying the appropriateness of our study to solar wind
heating.111In fact, the observations of Ref. 35 show a somewhat _larger_
amplitude (less anisotropic) than this estimate; but our simulation also has a
larger amplitude, because, due to the helicity barrier, the outer-scale
amplitude grows noticeably beyond the estimated $\delta B_{\perp}/B_{0}\simeq
L_{\perp}/L_{\parallel}$.  We also note that the elongated simulation domain
implies that the inferred outer-scale — the scale where $\delta B_{\perp}\sim
B_{0}$ — has an extremely long turnover time
$\simeq(L_{\|}/L_{\perp})^{2}\tau_{A}$ (where $\tau_{\rm A}=L_{\|}/v_{\rm
A}$), which is also comparable to the heating time $\tau_{\rm
heat}\sim(3/2m_{i})T_{i}V/\varepsilon$ ($V=L_{\perp}^{2}L_{\|}$ is the
volume). Because these outer-scale timescales exceed the duration of the
simulation (${\simeq}18\tau_{\rm A}$), it is apt to consider the outer-scale
forcing in our simulation as representing a turbulent flux of energy arriving
from larger scales, even if the outer scales are decaying36 (as relevant the
solar wind). Similarly, we note that the direct effect of solar-wind expansion
is negligible at these scales. A simple estimate can be obtained by matching
the simulation’s ion-inertial scale to that in the solar wind, defining the
expansion time as $\tau_{\rm exp}\sim R/U$, where $U$ the bulk solar-wind
velocity and $R$ the heliocentric radius. Using parameters similar to PSP’s
first perihelion7, 57, which also had $\beta_{i}\approx 0.3$, yields
$\tau_{\rm A}/\tau_{\rm exp}\simeq 8\times 10^{-4}(U/350{\rm
kms}^{-1})(R/35R_{\odot})^{-1}(B/80{\rm nT})^{-1}$, illustrating that the
expansion of the box is negligible ($\simeq 1.5\%$) over the simulation’s
duration.

Other parameters of the imbalanced simulation are chosen to match the balanced
turbulence simulation of Ref. 19, which had a resolution of
$N_{\perp}^{2}\times N_{\|}=200^{2}\times 1200$, a smaller box with $k_{\perp
0}\rho_{i}\equiv 2\pi\rho_{i}/L_{\perp}\approx 0.1$, and saturated amplitude
$\delta B_{\perp}/B_{0}\approx 0.14\simeq L_{\perp}/L_{\|}=1/6$. The energy-
injection rate is computed from $\varepsilon/V=C_{A}\delta u_{\perp
0}^{2}/\tau_{\rm A}$, with $\delta u_{\perp 0}/v_{\rm A}=L_{\perp}/L_{\|}$ and
a Kolmogorov constant $C_{A}=0.29$ to match that measured from the balanced
simulation. This implies a lower energy injection per unit volume in the
imbalanced case, compensating for the slower outer-scale motions in its larger
box; stated differently, with this $C_{A}$, turbulence in the larger $k_{\perp
0}\rho_{i}\approx 0.05$ box would saturate with $\delta B_{\perp}/B_{0}\approx
0.14$ if the forcing were balanced, which also implies it would have a smaller
amplitude (smaller $k_{\|}/k_{\perp}$) at $k_{\perp}\rho_{i}\sim 1$. The
simulation is initialized by randomly drawing particle velocities from a
stationary Maxwellian distribution with temperature $T_{i}=T_{e}$. We use
$N_{\rm ppc}=216$ particles per cell and the full-f method 34, with particles
initially evenly distributed within the cell. In computing the gridded moments
of $f_{i}$, we use two filter passes in order to reduce the impact of the
particle noise 34.

The value of the hyper-resistivity, $\eta_{4}\approx 2.4\times
10^{-5}d_{i}^{4}\Omega_{i}$, is also chosen to match that of the balanced
simulation (here $\Omega_{i}$ is the ion gyrofrequency and $d_{i}$ is the ion
inertial length). Its value is not intended to represent reality, but just to
absorb magnetic energy that cascades to the grid scales of the simulation. We
have tested the impact of this choice by restarting the simulation in the
saturated regime (at $t\simeq 16.3\tau_{\rm A}$) with $\eta_{4}\approx
1.6\times 10^{-5}d_{i}^{4}\Omega_{i}$. This modification extends
$\mathcal{E}_{B_{\perp}}(k_{\perp})$ to smaller scales where it flattens
further, as expected, while changing $\varepsilon_{\eta}$ only slightly and
making no noticeable difference to the diagnostics presented in Figs. 3 and 4
of the main text. We are thus confident that the chosen $\eta_{4}$ is
appropriate and that $\varepsilon_{\eta}$ is not adversely affected by grid-
scale effects.

#### Forcing

The plasma is driven at the largest scales in the box with the forcing terms
$\bm{F}_{\perp}^{U}$ and $\bm{F}_{\perp}^{B}$ in equations (1a) and (1b),
respectively. Given that our elongated box is designed to represent a small
patch of a much larger system, these terms are supposed to mimic crudely the
effect of the larger-scale turbulence on our box’s outer scale. As in Ref. 32,
we opt to design the forcing to inject energy and cross helicity at a constant
rate in time, which necessarily requires adjusting $\bm{F}_{\perp}^{U}$ and
$\bm{F}_{\perp}^{B}$ to respond to the state of the plasma. We thus define
$\bm{F}_{\perp}^{U}=f^{U}\mbox{\boldmath{$F$}}_{0}$ and
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\bm{F}_{\perp}^{B}=f^{B}\mbox{\boldmath{$F$}}_{0}$,
where the forcing function $\mbox{\boldmath{$F$}}_{0}$ is divergence-free and
purely perpendicular to $\mbox{\boldmath{$B$}}_{0}$ (no $z$ component), which
implies the forcing is nearly purely Alfvénic (its compressive part is small).
$\bm{F}_{0}$ is evolved in time as an Ornstein-Uhlenbeck process with
correlation time $t_{\rm corr}=\tau_{\rm A}/2$ for each mode with
$2\pi/L_{j}\leq k_{j}\leq 4\pi/L_{j}$, where $j$ represents each direction
$x$, $y$, and $z$. At each time step, we compute
$n\bm{u}_{i}\,\mbox{\boldmath{$\cdot$}}\,\bm{F}_{\perp}^{U}$ and
$\bm{B}\,\mbox{\boldmath{$\cdot$}}\,(\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\times$}}\,\bm{F}_{\perp}^{B})$
and adjust the values of $f^{U}$ and $f^{B}$ so as to make the injected energy
and cross helicity equal to their desired values ($\varepsilon$ and
$\varepsilon_{H}$, respectively, with $\varepsilon_{H}=0.9\varepsilon$). This
process requires inverting the curl to compute $\bm{F}_{\perp}^{B}$, which is
carried out using a Fourier transform; by adding the magnetic force in this
way and evolving $\bm{B}$ using the standard constrained-transport algorithm
of Pegasus++, we ensure
$\mbox{\boldmath{$\nabla$}}\,\mbox{\boldmath{$\cdot$}}\,\bm{B}=0$ to machine
precision.

This forcing method is an extension (to allow for imbalance) of the default
routines implemented in the Athena code 58, and has thus been used in a number
of previous works (e.g., Ref. 59). We have tested it by measuring separately
the energy and cross-helicity injection, which agree almost perfectly with the
input values, and by testing the full energy budget of the simulation (see
below). A possible downside of the method is that the forcing normalization,
$f^{U}$ and $f^{B}$, can change more rapidly than the spatial form of the
force $\bm{F}_{0}$. In the simulation, we see occasional sudden changes in
$f^{U}$ and $f^{B}$ that seem to be caused by the plasma flow and magnetic-
field perturbations (which dominantly propagate in the
$-\hat{\mbox{\boldmath{$z$}}}$ direction) moving out of phase with the large-
scale spatial structure of the forcing (determined by the slow evolution of
$\bm{F}_{0}$). In order to assess the impact of this effect, we restarted the
simulation at $t=10\tau_{\rm A}$ with a modified version of the forcing that
limited the change in $f^{U}$ and $f^{B}$ across a timestep $\delta t$ to
$\exp(\pm\Omega_{i}\delta t/10)$, where the $\pm$ was decided based on whether
$f^{U}$ and $f^{B}$ were above or below the optimal value that gave the input
energy and cross-helicity injection. This reduced the sudden changes (high-
frequency power) in $f^{U}$ and $f^{B}$ at the cost of causing
$\varepsilon(t)$ and $\varepsilon_{H}(t)$ to vary significantly (by around
$\pm 40\%$) in time. This change, however, made no noticeable difference to
the heating, distribution function, and spectra, lending us confidence in the
robustness of our results. Finally, we have confirmed that the helicity
barrier also forms robustly in the reduced model of Ref. 32 when the outer
scales are forced with white-in-time noise (as opposed to with constant energy
and cross-helicity injection). Thus, we do not expect our results to be
particularly sensitive to the design of the forcing.

Extended Data Fig. 1: Electric field noise and numerical cooling.
Contributions to the energy budget per unit volume of the imbalanced
simulation from the energy injection ($\varepsilon/V$; dashed red line),
increase in thermal energy $Q_{i}=\partial_{t}\langle m_{i}v_{\rm
th}^{2}/2\rangle$ (blue line), growth rate of mechanical energy
$\partial_{t}\langle nm_{i}u_{i}^{2}/2+B^{2}/8\pi\rangle$ (green line), and
resistive dissipation $\varepsilon_{\eta}/V$ ($V$ is the volume and
$\langle\,\dots\,\rangle$ denotes a box average). The black line shows the
total energy budget ${\rm
Total}=\varepsilon/V-\varepsilon_{\eta}/V-\partial_{t}\langle m_{i}v_{\rm
th}^{2}/2\rangle-\partial_{t}\langle nm_{i}u_{i}^{2}/2+B^{2}/8\pi\rangle$,
which is constant and negative, indicating numerical cooling that is
effectively independent of the turbulence or the heating of ions.  Extended
Data Fig. 2: The effect of particle noise on turbulence spectra. Perpendicular
($k_{\perp}$) spectra of the magnetic field ($\mathcal{E}_{\bm{B}}$), electric
field $\mathcal{E}_{\bm{E}}$, and KAW-normalized density $\mathcal{E}_{n_{\rm
KAW}}=\beta_{i}(1+2\beta_{i})\mathcal{E}_{n}$ in the saturated state (solid
lines) and at very early times (averaged over $t\leq 0.2\tau_{\rm A}$). The
latter is from before the turbulence has developed and is thus a proxy for the
noise floor in a given quantity. At the smallest scales,
$k_{\perp}\rho_{i}\gtrsim 3$, spectra are only modestly above the noise floor
and therefore uncertain.

#### Electric-field noise and numerical cooling

A persistent problem for the particle-in-cell method is the influence of
electric-field noise, which arises due to random density noise from the finite
number of numerical particles. In the full-f hybrid-kinetic method with the
Pegasus++ algorithm, a key impact of this noise is a numerical cooling, which
slowly drains thermal energy from the system. This was removed from Fig. 1 of
the main text for clarity, but can be accurately assessed by computing the
thermal energy budget, which is shown in Extended Data Fig. 1. The sum of the
various contributions to the total rate of change of energy, shown by the
thick black line, would equal zero if energy were conserved, but is instead
negative, indicating numerical cooling. This feature — although clearly
undesirable — can be avoided only by increasing the number of particles per
cell (thus increasing computational cost), or by using more filter passes
(thus decreasing the effective dynamic range). Its nearly constant value
throughout the simulation — including at very early times $t\lesssim\tau_{\rm
A}$, when the forcing is fully absorbed by large-scale mechanical energy and
there is no resistive dissipation — suggests its properties are mostly
separate from the development of turbulent heating in the system, justifying
tolerating its presence and removing it from our energy estimates.

Another effect of the particle noise is its direct influence on the spectra.
Noisy fluctuations cause an artificial bump in all quantities near the grid
scale, which accounts for some of the small-scale flattening of the electric-
and magnetic-field spectra in Fig. 2 of the main text. This can be quantified
by computing spectra from the very early stages of the simulation, before the
turbulence has developed, when the small-scale fluctuations are purely a
result of particle noise. We compare these to the saturated-turbulence spectra
in Extended Data Fig. 2. The general conclusion is that for
$k_{\perp}\rho_{i}\gtrsim 3$, noise makes a reasonable contribution to the
spectra, which are thus uncertain in this range. This effect is more severe in
this imbalanced simulation than other previous turbulence simulations with
Pegasus++ 19, 20 precisely because of the transition range, which makes the
sub-$\rho_{i}$ fluctuations very small in magnitude. Again, it can only be
ameliorated by increasing the number of particles per cell, which yields small
gains for large computational expense (the noise scales ${\propto}N_{\rm
ppc}^{1/2}$). If one subtracts the noise spectra from the saturated spectra —
likely a reasonable procedure because the noise power is dominated by modes of
high $k_{\perp}$ and $k_{z}$, while the majority of high-$k_{\perp}$ power at
saturation resides at lower $k_{z}$ — our main conclusions still hold, with
$\mathcal{E}_{\bm{B}}$ and $\mathcal{E}_{\bm{E}}$ re-flattening at
$k_{\perp}\rho_{i}\simeq 2$ and a ${\sim}k_{\perp}^{-0.8}$ KAW range in
$\mathcal{E}_{\bm{E}}$ (not shown).

Finally, we note that a different preliminary simulation with three-times
larger $\varepsilon$, four filter passes, and $N_{\rm ppc}=128$ reproduced the
key early time features discussed above and in the main text, including a
double-kinked electric-field spectrum (this simulation was run only until
$t\simeq 5\tau_{\rm A}$, which is before magnetic-field spectrum develops a
double kink). In addition to causing the large-scale energy to grow faster in
time, the larger $\varepsilon$ in this case caused the noise in spectra and
the numerical cooling to be proportionally much smaller. We are thus confident
that the physical features reported in the main text are robust.

### Measurement of the parallel spectrum

In Figs. 2 and 3 of the main text, we measure parallel spectra using a field-
line following method, which we describe here. While structure-function
methods are more commonly used to study anisotropic MHD turbulence 37, the
extremely steep parallel spectra (up to $\sim k_{\|}^{-6}$) caused by helicity
barrier are not well captured by structure functions 60. Our field-line-
following method computes a ($k_{\perp}$, $k_{\|}$) spectrum by first
constructing $N_{\rm lines}$ magnetic-field lines by solving ${\rm
d}\bm{r}/{\rm d}s=\hat{\mbox{\boldmath{$b$}}}(\bm{r})=\bm{B}/|\bm{B}|$ from
$s=0$ to $s=L_{\rm lines}$, where the periodicity of the system allows $L_{\rm
lines}>L_{\|}$ if desired. The field of interest $X(\bm{x})$ (e.g., $X=B_{y}$
or $X=E_{z}$) is then Fourier-filtered to a given bin in
$k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$, giving $X_{k_{\perp}}(\bm{x})$, which
is then interpolated onto the coordinates $\bm{r}(s)$ of the previously
computed field lines. One then computes the spectrum of
$X_{k_{\perp}}[\bm{r}(s)]$ in the field-line ($s$) direction, which becomes a
single $k_{\perp}$ slice of $\mathcal{E}_{X}(k_{\perp},k_{\|})$. Repeating
this process across a grid of $k_{\perp}$ yields the full 2D spectrum.222If
this process is applied to compute a $(k_{\perp},k_{z})$ spectrum [i.e.,
$\bm{r}(s)$ in the $\hat{\mbox{\boldmath{$z$}}}$ direction], it yields the
same result as computing the spectrum in the standard way with a 3D Fourier
transform. Through experimentation, we found $L_{\rm lines}=10L_{\|}$ and
$N_{\rm lines}=N_{x}N_{y}L_{\|}/L_{\rm lines}$ gives high-quality results
(though results are almost independent of $L_{\rm lines}$ for $L_{\rm
lines}\gtrsim 2L_{\|}$). In addition, a Hamming filter is used to compute the
$k_{\|}$ spectra because the field-line ($s$) direction is non-periodic.
Because the field-perpendicular plane is assumed to be the $x,y$ plane, the
method is valid in the “reduced-MHD” limit of $L_{\perp}/L_{\|}\sim\delta
B_{\perp}/B_{0}\ll 1$. A pure parallel spectrum, as in the inset of Fig. 2 of
the main text, is computed from $\int{\rm
d}k_{\perp}\,\mathcal{E}_{X}(k_{\perp},k_{\|})$, which recovers the standard
field-line parallel spectrum 61.

Extended Data Fig. 3: Measurement of the parallel spectrum. Two-dimensional
perpendicular magnetic-field spectrum
$\mathcal{E}_{B_{\perp}}(k_{\perp},k_{\|})=\mathcal{E}_{B_{x}}(k_{\perp},k_{\|})+\mathcal{E}_{B_{y}}(k_{\perp},k_{\|})$
from the balanced simulation. The method recovers the large- and small-scale
scalings of $k_{\|}$ with $k_{\perp}$ measured using structure functions (Ref.
19 figure 2), as well as the predicted 2D spectrum in the
$k_{\perp}\rho_{i}<1$ range (see Ref. 37, appendix B).

We have compared our method to the more conventional method of computing 2D
spectra via structure functions, finding reasonable agreement where both are
appropriate. Extended Data Fig. 3 shows
$\mathcal{E}_{B_{\perp}}(k_{\perp},k_{\|})=\mathcal{E}_{B_{x}}+\mathcal{E}_{B_{y}}$
for the balanced simulation, illustrating how we recovered the previously
discussed $k_{\|}\sim k_{\perp}^{2/3}$ and $k_{\|}\sim k_{\perp}$ ranges above
and below the ion-Larmor scale $k_{\perp}\rho_{i}\sim 1$ 19. Our method also
recovers the expected extremely steep scaling in $k_{\|}$
($\mathcal{E}_{B_{\perp}}\propto k_{\|}^{-6}$) predicted by theory at low
$k_{\perp}$ (see appendix B of 37; not shown). A more quantitative comparison
reveals that the most notable difference compared to the structure-function
method is that the 2D spectrum is shifted upwards somewhat (to larger
$k_{\|}$); thus, compared to the structure-function method, the field-line-
following method generally estimates $k_{\|}(k_{\perp})$ to be a factor of
$\simeq 2$ larger across all $k_{\perp}$.

Note that the helicity spectrum (Fig. 3C in the main text) was computed by
binning in $(k_{\perp},k_{z})$ rather than $(k_{\perp},k_{\|})$, which gave a
cleaner spectrum. The difficulty with $(k_{\perp},k_{\|})$ seems to relate to
$\mathcal{H}_{\bm{k}}$ explicitly involving $\bm{k}$ in its definition 62, 44,
which may cause problems with the non-orthogonal coordinate system
$(k_{x},k_{y},k_{\|})$. The $(k_{\perp},k_{z})$ spectrum of $B_{\perp}$ looks
broadly similar to its $(k_{\perp},k_{\|})$ spectrum, particularly in the
$k_{z}>k_{\perp}$ region of ICWs, so Fig. 3B and 3C of the main text can be
usefully compared on a qualitative level.

### Heating through quasi-linear diffusion on ICWs

#### Computation of the resonance contours for parallel ICWs

Here we describe in more detail the computation of the quasi-linear (QL) ICW
“scattering contours” shown in Fig. 4 of the main text. QL theory 24 is based
on the idea that particles interact strongly with a spectrum of ICWs if they
satisfy the cylotron-resonance condition

$\omega_{\bm{k}}-k_{\|}w_{\|}=\pm\Omega_{i}.$ (2)

Here $\omega_{\bm{k}}$ is the wave’s frequency and $w_{\|}$ is the ion’s
parallel velocity (both are measured in the fluid frame in which the plasma is
stationary); equation (2) expresses the condition for the Doppler-shifted wave
frequency to be resonant with the ion’s gyromotion, which causes the wave and
the particle to interact strongly. In the frame moving with the wave’s phase
velocity ${v}_{\rm ph}=\omega_{\bm{k}}/k_{\|}$, the magnetic-field
perturbation is constant in time. This means that the electric field must be
purely potential in this frame ($\bm{E}=-\mbox{\boldmath{$\nabla$}}\Phi$), and
a particle should approximately conserve its kinetic energy as it scatters
from the wave. Thus, if Eq. (2) has only one solution for each $w_{\|}$ — in
other words, when there exist only waves of a single $k_{\perp}$ (e.g., if
only parallel waves exist) or when $\omega_{\bm{k}}$ is a function of $k_{\|}$
only — particles diffuse only along specific “scattering contours” in the
$w_{\perp}$–$w_{\|}$ plane. If the waves are non-dispersive (i.e., if
$\omega_{\bm{k}}\propto k_{\|}$) these are semi-circles of constant particle
energy in the wave’s frame $\epsilon_{K,{\rm
wave}}=(w_{\|}-\omega_{\bm{k}}/k_{\|})^{2}+w_{\perp}^{2}$. More generally, the
scattering contours are defined by the null solutions
$\eta(w_{\perp},w_{\|})={\rm const.}$ of the QL diffusion operator 24, 47, 25,

$\left[1-\frac{w_{\|}}{v_{\mathrm{ph}}(w_{\|})}\right]\frac{\partial\eta}{\partial
w_{\perp}}+\frac{w_{\perp}}{v_{\rm ph}(w_{\|})}\frac{\partial\eta}{\partial
w_{\parallel}}=0,$ (3)

where $v_{\rm ph}(w_{\|})=\omega_{\bm{k}}/k_{\|}$ for the $w_{\|}$ that
satisfies the resonant condition (2). If $f_{i}$ is a decreasing function
along the scattering contours, the QL diffusion heats the plasma in the range
of $w_{\|}$ where there is appreciable power in waves of the corresponding
resonant $k_{\|}$; if $f_{i}$ is an increasing function along the scattering
contours, the distribution function will become unstable, growing the wave
power in the resonant range of $\bm{k}$.

In order to compute the scattering contours, we assume the cold-plasma ICW
dispersion relation for parallel propagating waves ($k_{\perp}=0$), which is

$\omega_{\bm{k}}=k_{\|}v_{\rm
A}\sqrt{1-\omega_{\bm{k}}/\Omega_{i}}\quad\Rightarrow\quad v_{\rm
ph}(k_{\|})=\frac{v_{\rm A}}{2\Omega_{i}}\left(\sqrt{k_{\|}^{2}v_{\rm
A}^{2}+4\Omega_{i}^{2}}-k_{\|}v_{\rm A}\right).$ (4)

A parametric solution for $\eta$ can be found by solving Eqs. (2) and (4) for
$v_{\rm ph}(w_{\|})$, and then using this in the solution of Eq. (3),
$w_{\perp}^{2}+w_{\|}^{2}-2\int{\rm d}w_{\|}^{\prime}\,v_{\rm
ph}(w_{\|}^{\prime})={\rm const}$. This yields the contours 63

$\frac{w_{\perp}^{2}}{v_{\rm A}^{2}}+\frac{1}{y^{2}}+\ln
y-\sinh^{-1}\frac{y}{2}={\rm const},$ (5)

where $y=k_{\|}v_{\rm A}/\Omega_{i}$ relates to $w_{\|}$ implicitly through
Eq. (2) (the full explicit solution is uninformatively complex).

The QL diffusion process relies on maintaining wave power in the relevant
range of $k_{\|}$, but the simulation exhibits a steep drop in the spectrum
for $k_{\|}d_{i}\gtrsim 1$ (see Fig. 3 of the main text).333This is likely
because linear ICWs become strongly damped at $k_{\|}d_{i}\simeq 1$ at
$\beta_{i}=0.3$ (as can be shown by solving the hot-plasma dispersion
relation). This implies a minimum $|w_{\|}|=w_{\|,d_{i}}$ above which $f_{i}$
should tend to flatten along the scattering contours and below which it should
not. This is computed from Eqs. (4) and (2) by solving for $w_{\|}$ at
$k_{\|}=d_{i}^{-1}$, giving

$w_{\|,d_{i}}\approx-\frac{1}{2}\left(3-\sqrt{5}\right)v_{\rm A}.$ (6)

The contours (5) and the cutoff (6), which are plotted in Fig. 4A of the main
text, provide an exceptionally good match to the shape of $f_{i}$.

#### Oblique ICWs

The above calculation assumes that only parallel waves exist in the plasma. In
the presence of oblique waves, the dependence of $v_{\rm ph}(k_{\|})$ on
$k_{\perp}$ means that a range of scattering contours exists for a given
$w_{\|}$. Refs. 47 and 64 provide a compelling argument as to why it is the
parallel ICW scattering contours that should determine the form of $f_{i}$,
even if oblique modes provide the primary heating power (see also Ref. 25).
They note that the dependence of $v_{\rm ph}(k_{\|})$ on $k_{\perp}$ is such
that higher-$k_{\perp}$ modes produce scattering contours that are steeper
[i.e., $(\partial\eta/\partial w_{\perp})/(\partial\eta/\partial w_{\|})$ is
larger]. A QL diffusion process along these contours will thus produce an
$f_{i}$ that is an increasing function along the resonance contours of
parallel modes [Eq. (3)], which will be unstable and generate parallel ICWs. A
spectrum of driven oblique modes will therefore generate parallel ICWs in the
process of heating — effectively a kinetic mechanism for spectral transfer
from oblique to parallel modes — creating an $f_{i}$ that is almost flat along
the _parallel_ ICW resonance contours.

This phenomenology provides a reasonable explanation for the behavior that we
observe in our simulation: in the saturated stage, there appears a significant
population of parallel ($k_{\perp}\approx 0$) ICWs (see main text Fig. 3B),
even though the turbulent cascade of energy to such modes is likely slow
compared to the power input into oblique modes (indeed, the power in parallel
modes is quite small earlier in the simulation).444A careful examination shows
that the contours of $f_{i}$ are very slightly steeper than the parallel-ICW
scattering contours, which should also be expected from this phenomenology:
the shape of $f_{i}$ results from a balance between heating from oblique modes
(which occurs on contours that are considerably steeper than those of $f_{i}$)
and cooling/instability from parallel modes. However, we caution that this
difference could easily have other explanations: e.g., differences between the
$\beta_{i}=0.3$ and cold-plasma ICW dispersion relations or higher-order
resonances.  This also suggests that the ICW heating process can continue even
after $f_{i}$ becomes perfectly flat along the parallel ICW resonance
contours. This is an important feature of this heating process for application
to the solar wind and seems to be what we observe in the saturated state:
after $t\approx 14\tau_{\rm A}$ heating continues but $f_{i}$ expands outwards
across the scattering contours. Potential complications with this scenario
arise from higher-order resonances (which are possible with oblique but not
with parallel ICWs) and electron damping by ICWs at $k_{\|}d_{i}\ll 1$ 65, but
these are generally expected to be unimportant to the overall physics 47.

Finally, we emphasize that by the end of the simulation, the average thermal
anisotropy of the full distribution is only $T_{\perp}/T_{\|}\approx 1.03$.
This value is small for two main reasons: first, because $\tau_{\rm
heat}\simeq 54\tau_{\rm A}$ is longer than the simulation duration and second,
because of the development of the beam, which contributes to the parallel
temperature. Given that this $T_{\perp}/T_{\|}$ is far below the usual bi-
Maxwellian ICW instability threshold of $T_{\perp}/T_{\|}\approx 1.7$ at
$\beta\approx 0.3$ that is often used in observational studies66, it is clear
that the detailed shape of the distribution function must be considered in
order to understand its stability to parallel ICWs and other wave modes.

#### Computation of $D_{\perp\perp}^{E}$

The perpendicular energy diffusion coefficient provides a useful quantitative
diagnostic of the plasma heating mechanism. It is defined by assuming that

$\frac{\partial f_{i}}{\partial t}=\frac{\partial}{\partial
e_{\perp}}\left(D^{E}_{\perp\perp}\frac{\partial f_{i}}{\partial
e_{\perp}}\right),$ (7)

where $e_{\perp}=w_{\perp}^{2}/2$ is the perpendicular kinetic energy per unit
mass. The assumption that the evolution of $f_{i}$ is described by (7), which
is simply a diffusion equation in perpendicular energy, hinges on
perpendicular heating dominating over parallel heating, as is indeed the case
in our simulation (at least in the $w_{\|}\lesssim 0$ part of velocity space).
In standard QL theory (ignoring the oblique–parallel mode interaction
discussed above), the resonance condition depends only on the parallel
velocity and the fluctuation spectrum, implying that the velocity diffusion
coefficient is independent of $w_{\perp}$ and thus that the energy diffusion
coefficient scales as $D_{\perp\perp}^{E}\propto w_{\perp}^{2}$. Numerically,
we compute $D_{\perp\perp}^{E}$ directly from Eq. (7) using

$D_{\perp\perp}^{E}=\left(\frac{\partial f_{i}}{\partial
e_{\perp}}\right)^{-1}\int_{0}^{e_{\perp}}{\rm
d}e_{\perp}^{\prime}\,\frac{\partial f_{i}}{\partial t},$ (8)

with the full expression integrated over a range of $w_{\|}$ (this is method
II of Ref. 48). We have also compared this result to a direct measurement of
the perpendicular heating from $D_{\perp\perp}^{E}=-(\partial
Q_{\perp}/\partial e_{\perp})/({\partial f_{i}}/{\partial e_{\perp}})$, where
$\partial Q_{\perp}/\partial e_{\perp}$ is a direct measure of the heating of
particles by electric fields based on $\partial Q_{\perp}/\partial
w_{\perp}=e\langle\bm{w}_{\perp}\,\mbox{\boldmath{$\cdot$}}\,\bm{E}_{\perp}f_{i}\rangle$,
and
$\langle\bm{w}_{\perp}\,\mbox{\boldmath{$\cdot$}}\,\bm{E}_{\perp}f_{i}\rangle$
is computed from ion velocities and electric fields during the simulation.
This measurement, although noisier, recovers very similar results to Eq. (8).
Further discussion can be found in Refs. 48 and 20, in particular in appendix
A of Ref. 20.

### Stochastic heating

Extended Data Fig. 4: Assessment of the influence of stochastic heating. We
show perpendicular spectra of the electric potential $\Phi$, computed from the
curl free part of $E$. Colored lines show various times from the imbalanced
simulation. The black line shows the equivalent balanced simulation, which is
averaged over the early period of the simulation (between $t=3.5\tau_{\rm A}$
and $t=4.5\tau_{\rm A}$) when stochastic-ion heating absorbs the majority of
the turbulent energy flux 19. Despite the larger turbulence amplitude in the
imbalanced simulation, the electric-potential fluctuations around
$k_{\perp}\rho_{i}\sim 1$ — those important for stochastic heating — are
smaller.

Here we address whether, in addition to quasi-linear ICW heating, stochastic
ion heating 16 might play an important role in turbulence with a helicity
barrier. The mechanism is of particular interest, given its prominence in
previous theoretical and observational studies 67, 68, 20. Possible scenarios
could involve multiple heating mechanisms operating at a particular time, or a
transition from one heating mechanism to another as $f_{i}$ changes shape in
time. We find no evidence for such behavior in this simulation:
$D_{\perp\perp}^{E}$ appears to maintain its QL scaling
$D_{\perp\perp}^{E}\propto w_{\perp}^{2}$ throughout, and $f_{i}$ does not
deviate from the scattering contours as it evolves in the saturated state.

We speculate that in this simulation, the lack of stochastic heating is simply
a consequence of its small electric-potential ($\Phi$) fluctuations around
$k_{\perp}\rho_{i}\sim 1$, which are required in order to make ion gyro-orbits
sufficiently random to cause heating 16. In Extended Data Fig. 4, we compare
the spectrum of $\Phi$ ($\mathcal{E}_{\Phi}$) in the imbalanced simulation and
the balanced simulation 19. $\mathcal{E}_{\Phi}(k_{\perp}\rho_{i}=1)$ grows
only modestly during the imbalanced simulation, despite the growth of
$\mathcal{E}_{\Phi}$ at larger scales, because of the steep drop at
$k_{\perp}\rho_{i}<1$ due to the helicity barrier. Coupled with its larger
box, we see that even though the imbalanced simulation saturates with larger
amplitude turbulence, its ion-gyroscale $\Phi$ fluctuations are smaller than
in the balanced run. Given that stochastic heating plays only a modest role in
this balanced run — it becomes subdominant after several turnover times due to
flattening of the core of $f_{i}$ 19, 20 — this difference in $\Phi$ may be
sufficient to render stochastic heating unimportant in imbalanced turbulence
at these parameters.

It is unclear whether stochastic heating will play a role in other regimes or
over longer timescales. So long as gyroscale fluctuations have sufficient
amplitudes, stochastic heating is expected to be more robust at lower $\beta$
because more heating occurs before it is quenched by the flattening of the
core of $f_{i}$ 16. This conclusion supported by the $\beta\approx 1/9$ hybrid
simulation of Ref. 20. On the other hand, test-particle simulations show a
strong reduction in the efficiency of heating of $w\approx v_{\rm A}$
($\beta\approx 1$) particles in imbalanced, compared to balanced, turbulence
21. If a similar reduction occurs also in the $w\ll v_{\rm A}$ low-$\beta$
regime, the effectiveness of stochastic heating in turbulence with a helicity
barrier may also be limited. Further work is needed. However, it is worth
noting that even if stochastic heating, rather than ICW heating, eventually
absorbs the turbulent energy flux, the helicity barrier could remain a key
ingredient in solar-wind heating. Just like for ICW heating, the barrier would
allow turbulent fluctuations to grow their amplitude sufficiently to enable
ion heating, rather than fluctuations dissipating their energy into electron
heating if their amplitude is initially too small.

### Landau damping and the ion beam

An interesting feature of the ion distribution function shown in Fig. 4 of the
main text is the plateau $w_{\perp}\approx 0$, $w_{\|}\sim v_{\rm A}$. This
forms a modestly super-Alfvénic beam with similar properties 49, 2 and
directionality 69 to those observed in the solar wind. Here we present
evidence that this feature is a result of Landau damping of perpendicular
Alfvén waves as they become dispersive (speed up) near $k_{\perp}\rho_{i}\sim
1$. Test-particle calculations have shown this process to be highly effective
51.

Extended Data Fig. 5: Development of the ion beam. We compare the rate of
change parallel thermal energy (solid lines; see text) with the work done on
particles by the parallel electric field $e\langle w_{\|}E_{\|}f_{i}\rangle$
(dotted lines). The thick dark-blue lines show the saturated state and the
orange-pink lines show $t=7\tau_{\rm A}$. The similarity of the magnitude and
general shape of the two measures of heating suggests that Landau damping is
responsible for the formation of the ion beam.

Extended Data Fig. 5 compares the measured parallel heating of the
distribution function,

$\left\langle\frac{\partial e_{\|}f_{i}}{\partial
t}\right\rangle(w_{\|})=\int{\rm
d}e_{\perp}\,\frac{1}{2}w_{\|}^{2}\frac{\partial f_{i}}{\partial t},$ (9)

to the parallel heating inferred from the work done by the parallel electric
field $e\langle w_{\|}E_{\|}f_{i}\rangle$, which is computed during the
simulation from particle trajectories. We see a clear peak in both quantities
at $w_{\|}\approx v_{\rm A}$. As discussed in Methods (see equation 8
therein), both methods measure plasma heating but they can differ by a total
derivative. Their similarity in Extended Data Fig. 5 — in particular the
similarity of their magnitudes even at different times during the simulation
when $\langle\partial(e_{\|}f_{i})/\partial t\rangle$ differs — suggests that
parallel-electric-field work, _viz.,_ Landau damping, is responsible for the
formation of the beam.

## Acknowledgments

We thank Bill Dorland, Ben Chandran, and Alfred Mallet for illuminating
discussions. J.S. and R.M acknowledge support from the Royal Society Te
Apārangi, New Zealand through Marsden Fund grant UOO1727 and Rutherford
Discovery Fellowship RDF-U001804. M.W.K. and E.Q. were supported by the
Department of Energy through the NSF/DOE Partnership in Basic Plasma Science
and Engineering, awards DE-SC0019046 and DE-SC0019047, with additional support
for E.Q. from a Simons Investigator Award from the Simons Foundation. L.A.
acknowledges the support of the Institute for Advanced Study, and the work of
A.A.S. was supported in part by UK EPSRC grant EP/R034737/1. This research was
part of the Frontera computing project at the Texas Advanced Computing Center,
which is made possible by National Science Foundation award OAC-1818253.
Further computational support was provided by the New Zealand eScience
Infrastructure (NeSI) high performance computing facilities, funded jointly by
NeSI’s collaborator institutions and through the NZ MBIE, and through PICSciE-
OIT TIGRESS High Performance Computing Center and Visualization Laboratory at
Princeton University. The funders had no role in study design, data collection
and analysis, decision to publish or preparation of the manuscript.

## References

  * 1 S. R. Cranmer and A. R. Winebarger.  The Properties of the Solar Corona and Its Connection to the Solar Wind.  Ann. Rev. Astron. Astrophys., 57:157–187, (2019).
  * 2 E. Marsch.  Kinetic physics of the solar corona and solar wind.  Living Rev. Solar Phys., 3: 1(2006).
  * 3 E. N. Parker.  Dynamical Theory of the Solar Wind.  Space Sci. Rev., 4(5-6):666–708, (1965).
  * 4 V. H. Hansteen and E. Leer.  Coronal heating, densities, and temperatures and solar wind acceleration.  J. Geophys. Res., 100(A11):21577–21594, (1995).
  * 5 J. L. Kohl, G. Noci, E. Antonucci, et al.  First Results from the SOHO Ultraviolet Coronagraph Spectrometer.  Solar Phys., 175(2):613–644, (1997).
  * 6 S. R. Cranmer, G. B. Field, and J. L. Kohl.  Spectroscopic Constraints on Models of Ion Cyclotron Resonance Heating in the Polar Solar Corona and High-Speed Solar Wind.  Astrophys. J., 518(2):937–947, (1999).
  * 7 S. D. Bale, S. T. Badman, J. W. Bonnell, et al.  Highly structured slow solar wind emerging from an equatorial coronal hole.  Nature, 576:237–242, (2019).
  * 8 B. De Pontieu, S. W. McIntosh, M. Carlsson, et al.  Chromospheric Alfvénic Waves Strong Enough to Power the Solar Wind.  Science, 318(5856):1574, (2007).
  * 9 S. Tomczyk, S. W. McIntosh, S. L. Keil, et al.  Alfvén Waves in the Solar Corona.  Science, 317(5842):1192, (2007).
  * 10 M. Velli, R. Grappin, and A. Mangeney.  Turbulent cascade of incompressible unidirectional Alfvén waves in the interplanetary medium.  Phys. Rev. Lett., 63(17):1807–1810, (1989).
  * 11 A. A. van Ballegooijen, M. Asgari-Targhi, S. R. Cranmer, and E. E. DeLuca.  Heating of the Solar Chromosphere and Corona by Alfvén Wave Turbulence.  Astrophys. J., 736(1):3, (2011).
  * 12 M. Shoda, T. K. Suzuki, M. Asgari-Targhi, and T. Yokoyama.  Three-dimensional Simulation of the Fast Solar Wind Driven by Compressible Magnetohydrodynamic Turbulence.  Astrophys. J. Lett., 880(1):L2, (2019).
  * 13 E. Quataert and A. Gruzinov.  Turbulence and particle heating in advection-dominated accretion flows.  Astrophys. J., 520(1):248–255, (1999).
  * 14 A. A. Schekochihin, Y. Kawazura, and M. A. Barnes.  Constraints on ion versus electron heating by plasma turbulence at low beta.  J. Plasma Phys., 85(3):905850303, (2019).
  * 15 A. A. Schekochihin, S. C. Cowley, W. Dorland, et al.  Astrophysical gyrokinetics: Kinetic and fluid turbulent cascades in magnetized weakly collisional plasmas.  Astrophys. J. Supp., 182(1):310, (2009).
  * 16 B. D. G. Chandran, B. Li, B. N. Rogers, et al.  Perpendicular Ion Heating by Low-frequency Alfvén-wave Turbulence in the Solar Wind.  Astrophys. J., 720(1):503–515, (2010).
  * 17 B. D. G. Chandran, T. J. Dennis, E. Quataert, and S. D. Bale.  Incorporating Kinetic Physics into a Two-fluid Solar-wind Model with Temperature Anisotropy and Low-frequency Alfvén-wave Turbulence.  Astrophys. J., 743(2):197, (2011).
  * 18 D. Vech, K. G. Klein, and J. C. Kasper.  Nature of Stochastic Ion Heating in the Solar Wind: Testing the Dependence on Plasma Beta and Turbulence Amplitude.  Astrophys. J. Lett., 850(1):L11, (2017).
  * 19 L. Arzamasskiy, M. W. Kunz, B. D. G. Chandran, and E. Quataert.  Hybrid-kinetic Simulations of Ion Heating in Alfvénic Turbulence.  Astrophys. J., 879(1):53, (2019).
  * 20 S. S. Cerri, L. Arzamasskiy, and M. W. Kunz.  On Stochastic Heating and Its Phase-space Signatures in Low-beta Kinetic Turbulence.  Astrophys. J., 916(2):120, (2021).
  * 21 B. Teaca, M. S. Weidl, F. Jenko, and R. Schlickeiser.  Acceleration of particles in imbalanced magnetohydrodynamic turbulence.  Phys. Rev. E, 90(2):021101, (2014).
  * 22 P. A. Isenberg and B. J. Vasquez.  Perpendicular Ion Heating by Cyclotron Resonant Dissipation of Turbulently Generated Kinetic Alfvén Waves in the Solar Wind.  Astrophys. J., 887(1):63, (2019).
  * 23 J. V. Hollweg and P. A. Isenberg.  Generation of the fast solar wind: A review with emphasis on the resonant cyclotron interaction.  J. Geophys. Res.: Space Phys., 107(A7):1147, (2002).
  * 24 C. F. Kennel and F. Engelmann.  Velocity Space Diffusion from Weak Plasma Turbulence in a Magnetic Field.  Phys. Fluids, 9(12):2377–2388, (1966).
  * 25 P. A. Isenberg and B. J. Vasquez.  A Kinetic Model of Solar Wind Generation by Oblique Ion-cyclotron Waves.  Astrophys. J., 731(2):88, (2011).
  * 26 J. C. Kasper, B. A. Maruca, M. L. Stevens, and A. Zaslavsky.  Sensitive Test for Ion-Cyclotron Resonant Heating in the Solar Wind.  Phys. Rev. Lett., 110(9):091102, (2013).
  * 27 G. Q. Zhao, H. Q. Feng, D. J. Wu, et al.  Dependence of Ion Temperatures on Alpha-Proton Differential Flow Vector and Heating Mechanisms in the Solar Wind.  Astrophys. J. Lett., 889(1):L14, (2020).
  * 28 J. V. Hollweg.  Compressibility of ion cyclotron and whistler waves: Can radio measurements detect high-frequency waves of solar origin in the corona?  J. Geophys. Res., 105(A4):7573–7582, (2000).
  * 29 J. V. Shebalin, W. H. Matthaeus, and D. Montgomery.  Anisotropy in MHD turbulence due to a mean magnetic field.  J. Plasma Phys., 29(3):525–547, (1983).
  * 30 L. D. Woodham, R. T. Wicks, D. Verscharen, et al.  Parallel-propagating Fluctuations at Proton-kinetic Scales in the Solar Wind Are Dominated By Kinetic Instabilities.  Astrophys. J. Lett., 884(2):L53, (2019).
  * 31 Y. Voitenko and M. Goossens.  Excitation of high-frequency Alfvén waves by plasma outflows from coronal reconnection events.  Solar Phys., 206(2):285–313, (2002).
  * 32 R. Meyrand, J. Squire, A. A. Schekochihin, and W. Dorland.  On the violation of the zeroth law of turbulence in space plasmas.  J. Plasma Phys., 87(3):535870301, (2021).
  * 33 J. Cho.  Magnetic Helicity Conservation and Inverse Energy Cascade in Electron Magnetohydrodynamic Wave Packets.  Phys. Rev. Lett., 106(19):191104, (2011).
  * 34 M. W. Kunz, J. M. Stone, and X.-N. Bai.  Pegasus: A new hybrid-kinetic particle-in-cell code for astrophysical plasma dynamics.  J. Comp. Phys., 259:154–174, (2014).
  * 35 C. H. K. Chen.  Recent progress in astrophysical plasma turbulence from solar wind observations.  J. Plasma Phys., 82(6):535820602, (2016).
  * 36 P. A. Davidson.  Turbulence : An Introduction for Scientists and Engineers.  Oxford University Press, (2004).
  * 37 A. A. Schekochihin.  MHD Turbulence: A Biased Review.  arXiv:2010.00699, (2021).
  * 38 M. D. McManus, T. A. Bowen, A. Mallet, et al.  Cross Helicity Reversals in Magnetic Switchbacks.  Astrophys. J. Supp., 246(2):67, (2020).
  * 39 R. J. Leamon, C. W. Smith, N. F. Ness, et al.  Observational constraints on the dynamics of the interplanetary magnetic field dissipation range.  J. Geophys. Res., 103(A3):4775–4788, (1998).
  * 40 T. A. Bowen, A. Mallet, S. D. Bale, et al.  Constraining Ion-Scale Heating and Spectral Energy Transfer in Observations of Plasma Turbulence.  Phys. Rev. Lett., 125(2):025102, (2020).
  * 41 C. H. K. Chen, S. Boldyrev, Q. Xia, and J. C. Perez.  Nature of Subproton Scale Turbulence in the Solar Wind.  Phys. Rev. Lett., 110(22):225002, (2013).
  * 42 D. Duan, J. He, T. A. Bowen, et al.  Anisotropy of Solar Wind Turbulence in the Inner Heliosphere at Kinetic Scales: PSP Observations.  Astrophys. J. Lett., 915(1):L8, (2021).
  * 43 S. Y. Huang, F. Sahraoui, N. Andrés, et al.  The Ion Transition Range of Solar Wind Turbulence in the Inner Heliosphere: Parker Solar Probe Observations.  Astrophys. J. Lett., 909(1):L7, (2021).
  * 44 G. G. Howes and E. Quataert.  On the Interpretation of Magnetic Helicity Signatures in the Dissipation Range Of Solar Wind Turbulence.  Astrophys. J. Lett., 709(1):L49–L52, (2010).
  * 45 S. Y. Huang, J. Zhang, F. Sahraoui, et al.  Kinetic Scale Slow Solar Wind Turbulence in the Inner Heliosphere: Coexistence of Kinetic Alfvén Waves and Alfvén Ion Cyclotron Waves.  Astrophys. J. Lett., 897(1):L3, (2020).
  * 46 G. Q. Zhao, Y. Lin, X. Y. Wang, et al.  Magnetic Helicity Signature and Its Role in Regulating Magnetic Energy Spectra and Proton Temperatures in the Solar Wind.  Astrophys. J., 906(2):123, (2021).
  * 47 B. D. G. Chandran, P. Pongkitiwanichakul, P. A. Isenberg, et al.  Resonant Interactions Between Protons and Oblique Alfvén/Ion-cyclotron Waves in the Solar Corona and Solar Flares.  Astrophys. J., 722(1):710–720, (2010).
  * 48 B. J. Vasquez, P. A. Isenberg, and S. A. Markovskii.  Proton Perpendicular Heating in Turbulence Simulations: Determination of the Velocity Diffusion Coefficient.  Astrophys. J., 893(1):71, (2020).
  * 49 E. Marsch, R. Schwenn, H. Rosenbauer, et al.  Solar wind protons: Three-dimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU.  J. Geophys. Res., 87(A1):52–72, (1982).
  * 50 J. L. Verniero, D. E. Larson, R. Livi, et al.  Parker Solar Probe Observations of Proton Beams Simultaneous with Ion-scale Waves.  Astrophys. J. Supp., 248(1):5, (2020).
  * 51 X. Li, Q. Lu, Y. Chen, et al.  A Kinetic Alfvén Wave and the Proton Distribution Function in the Fast Solar Wind.  Astrophys. J. Lett., 719(2):L190–L193, (2010).
  * 52 Y. M. Wang and J. Sheeley, N. R.  Solar Wind Speed and Coronal Flux-Tube Expansion.  Astrophys. J., 355:726, (1990).
  * 53 S. R. Cranmer.  Why is the Fast Solar Wind Fast and the Slow Solar Wind Slow? (Invited) A Survey of Geometrical Models.  In B. Fleck, T. H. Zurbuchen, and H. Lacoste, editors, Solar Wind 11/SOHO 16, Connecting Sun and Heliosphere, volume 592 of ESA Special Publication, page 159, (2005).
  * 54 B. D. G. Chandran.  An approximate analytic solution to the coupled problems of coronal heating and solar-wind acceleration.  J. Plasma Phys., 87(3):905870304, (2021).
  * 55 J. A. Byers, B. I. Cohen, W. C. Condit, and J. D. Hanson.  Hybrid Simulations of Quasineutral Phenomena in Magnetized Plasma.  J. Comp. Phys., 27(3):363–396, (1978).
  * 56 P. Goldreich and S. Sridhar.  Toward a theory of interstellar turbulence. Strong Alfvénic turbulence.  Astrophys. J., 438:763–775, (1995).
  * 57 J. C. Kasper, S. D. Bale, J. W. Belcher, et al.  Alfvénic velocity spikes and rotational flows in the near-Sun solar wind.  Nature, 576:228–231, (2019).
  * 58 J. M. Stone, T. A. Gardiner, P. Teuben, et al.  Athena: A new code for astrophysical MHD.  Astrophys. J. Supp., 178(1):137–177, (2008).
  * 59 J. W. Lynn, I. J. Parrish, E. Quataert, and B. D. G. Chandran.  Resonance broadening and heating of charged particles in magnetohydrodynamic turbulence.  Astrophys. J., 758(2):78, (2012).
  * 60 J. Cho and A. Lazarian.  Simulations of Electron Magnetohydrodynamic Turbulence.  Astrophys. J., 701(1):236–252, (2009).
  * 61 R. Meyrand, A. Kanekar, W. Dorland, and A. A. Schekochihin.  Fluidization of collisionless plasma turbulence.  Proc. Nat. Acad. Sci., 116(4):1185–1194, (2019).
  * 62 W. H. Matthaeus and M. L. Goldstein.  Measurement of the rugged invariants of magnetohydrodynamic turbulence in the solar wind.  J. Geophys. Res., 87(A8):6011–6028, (1982).
  * 63 P. A. Isenberg and M. A. Lee.  A dispersive analysis of bispherical pickup ion distributions.  J. Geophys. Res., 101(A5):11055–11066, (1996).
  * 64 P. Pongkitiwanichakul, B. D. G. Chandran, P. A. Isenberg, and B. J. Vasquez.  Resonant Interactions Between Protons and Oblique Alfvén/Ion-Cyclotron Waves.  In M. Maksimovic, K. Issautier, N. Meyer-Vernet, et al., editors, Twelfth International Solar Wind Conference, volume 1216 of AIP Conf. Ser., pages 72–75, (2010).
  * 65 S. P. Gary and J. E. Borovsky.  Alfvén-cyclotron fluctuations: Linear Vlasov theory.  J. Geophys. Res.: Space Phys., 109(A6):A06105, (2004).
  * 66 P. Hellinger, P. Trávníček, J. C. Kasper, and A. J. Lazarus.  Solar wind proton temperature anisotropy: Linear theory and wind/swe observations.  Geophysical Research Letters, 33(9):L09101, (2006).
  * 67 K. G. Klein and B. D. G. Chandran.  Evolution of The Proton Velocity Distribution due to Stochastic Heating in the Near-Sun Solar Wind.  Astrophys. J., 820(1):47, (2016).
  * 68 M. M. Martinović, K. G. Klein, J. C. Kasper, et al.  The Enhancement of Proton Stochastic Heating in the Near-Sun Solar Wind.  Astrophys. J. Supp., 246(2):30, (2020).
  * 69 J. He, Z. Pei, L. Wang, et al.  Sunward Propagating Alfvén Waves in Association with Sunward Drifting Proton Beams in the Solar Wind.  Astrophys. J., 805(2):176, (2015).