# Local proton heating at magnetic discontinuities in Alfvénic and non-
Alfvénic solar wind

C.A. González  Department of Physics, The University of Texas at Austin,
Austin, TX USA  J. L. Verniero  Code 672, NASA, Goddard Space Flight Center,
Greenbelt, MD 20771, USA  R. Bandyopadhyay  Department of Astrophysical
Sciences, Princeton, New Jersey 08544, USA  A. Tenerani  Department of
Physics, The University of Texas at Austin, Austin, TX USA

###### Abstract

We investigate the local proton energization at magnetic
discontinuities/intermittent structures and the corresponding kinetic
signatures in velocity phase space in Alfvénic (high cross helicity) and non-
Alfvénic (low cross helicity) wind streams observed by Parker Solar Probe. By
means of the Partial Variance of Increments method, we find that the hottest
proton populations are localized around compressible, kinetic-scale magnetic
structures in both types of wind. Furthermore, the Alfvénic wind shows
preferential enhancements of $T_{\parallel}$ as smaller scale structures are
considered, whereas the non-Alfvenic wind shows preferential $T_{\bot}$
enhancements. Although proton beams are present in both types of wind, the
proton velocity distribution function displays distinct features. Hot beams,
i.e., beams with beam-to-core perpendicular temperature
$T_{\bot,b}/T_{\bot,c}$ up to three times larger than the total distribution
anisotropy, are found in the non-Alfvénic wind, whereas colder beams in the
Alfvénic wind. Our data analysis is complemented by 2.5D hybrid simulations in
different geometrical setups, which support the idea that proton beams in
Alfvénic and non-Alfvénic wind have different kinetic properties and different
origins. The development of a perpendicular nonlinear cascade, favored in
balanced turbulence, allows a preferential relative enhancement of the
perpendicular plasma temperature and the formation of hot beams. Cold field-
aligned beams are instead favored by Alfvén wave steepening. Non-Maxwellian
distribution functions are found near discontinuities and intermittent
structures, pointing to the fact that the nonlinear formation of small scale
structures is intrinsically related to the development of highly non-thermal
features in collisonless plasmas. Our results contribute to understanding the
role of different coherent structures in proton energization and their
implication in collisionless energy dissipation processes in space plasmas.

Solar wind, Alfven waves, Interplanetary turbulence, Interplanetary particle
acceleration

††journal: ApJ

##  1 Introduction

Magnetic discontinuities have been observed for a long time throughout the
heliosphere (Colburn & Sonett, 1966; Parker, 1994; Burlaga, 1991; Tsurutani &
Ho, 1999; Vasquez et al., 2007), and particularly in turbulent solar wind
streams where discontinuities have been related to intermittent structures
(Greco et al., 2016, 2008a). Magnetic discontinuities and intermittent
structures may play an important role in energy dissipation and transport, and
in particle acceleration (Osman et al., 2010, 2012; Tessein et al., 2013).
Several mechanisms have been proposed for solar wind heating, including
resonant wave-particle interactions (such as Landau damping and ion-cyclotron
resonance (Narita & Marsch, 2015; Gary & Borovsky, 2008; Sahraoui et al.,
2010; Bowen et al., 2022)), magnetic pumping (Lichko et al., 2017), stochastic
heating (Chen et al., 2001; Johnson & Cheng, 2001; Chandran et al., 2010; Vech
et al., 2017) and intermittent dissipation (Dmitruk et al., 2004; Karimabadi
et al., 2013; Osman et al., 2014). However, solar wind streams differ by the
type of fluctuations and turbulence, and which heating mechanisms are favored
under different solar wind conditions remains to be understood.

Tangential (TD) and rotational discontinuities (RD) are the most common types
of discontinuities observed in the solar wind (Neugebauer, 2006; Paschmann et
al., 2013). In theory, TDs and RDs can be identified depending on the normal
component of the magnetic field and the changes in plasma parameters across
the structure (Hudson, 1970). TDs are stationary structures in the plasma
frame and have zero normal magnetic field component. These structures
represent boundaries between different plasma parcels and are in pressure
balance, without restrictions on the variations of magnetic field magnitude
across the discontinuity. TDs are commonly found at reconnection sites or,
more generally, at the boundaries of flux tubes (Greco et al., 2009; Servidio
et al., 2011; Zhdankin et al., 2012). On the other hand, RDs have a non-zero
normal magnetic field component and a field-aligned flow. RDs are propagating
structures that have been associated with phase-steepened Alfvén waves
(Tsurutani et al., 1994; Medvedev et al., 1997; Vasquez & Hollweg, 2001).

The classification of TDs and RDs is not straightforward, not only due to
single-spacecraft limitations to estimate accurately the normal component, but
also because discontinuities are non-ideal and can display properties of both
TDs and RDs (Neugebauer et al., 1984; Horbury et al., 2001; Knetter et al.,
2004; Artemyev et al., 2019). Nevertheless, previous work has shown that
different types of magnetic structures are statistically detected in different
types of wind. Namely, compressible structures are most likely observed in the
slow solar wind, while structures with small compressibility are typically
observed in the fast solar wind (Perrone et al., 2016, 2017). Although non-
Alfvénic and Alfvénic solar wind is traditionally associated with slow and
fast streams, recent observations showed the existence of “slow Alfvenic” wind
(D’Amicis et al., 2021). This suggests that the solar wind cannot be simply
categorized based solely on its velocity. The normalized residual energy (the
relative energy in kinetic and magnetic fluctuations), the Alfvén ratio (the
ratio between kinetic to magnetic fluctuations), and the normalized cross
helicity (the correlation between magnetic and velocity field fluctuations)
are some of the measures that quantify Alfvénicity. In this work, we classify
the solar wind based on the average value of cross-helicity. We consider the
solar wind as Alfvénic when it contains fluctuations with cross helicity close
to unity. In contrast, the non-Alfvénic wind refers to fluctuations with cross
helicity approaching zero.

Several studies have shown that plasma temperature enhancements are associated
with intermittent structures/discontinuities identified via the partial
variance of increment (PVI) method (Osman et al., 2010; Qudsi et al., 2020;
Sioulas et al., 2022b; Phillips et al., 2023). The PVI method enables
identification of current sheet structures, but it does not distinguish
whether those structures are predominantly compressible or rotational.
However, TDs and RDs are fundamentally different types of structures that
affect particle energization in different ways, and how compressible (TD-like)
and incompressible (RD-like) structures in the solar wind interact with
particles ultimately leading to heating (and their velocity-space signatures)
remains to be investigated.

In this work, we address such a problem by investigating proton heating at
different types of intermittent structures and/or discontinuities in the inner
heliosphere, by using Parker Solar Probe (PSP) data and 2.5D hybrid-kinetic
simulations. Toward this goal, we conduct a comparison study between Alfvénic
and non-Alfvénic streams and investigate the dependence of temperature
anisotropies on Alfvénicity and on the type of magnetic structures.
Additionally, hybrid simulations, with simplified geometry settings, are used
to investigate what processes are likely to contribute to proton heating and
to the generation of the observed proton temperature anisotropy in each type
of wind.

This paper is organized as follows. In section 2, we describe the data and
methods, and we report our data analysis. An overview of the properties of
magnetic field fluctuations in Alfvénic and non-Alfvénic wind is reported in
section 2.1. In section 2.2 we discuss the correlation between PVI and proton
temperatures and temperature anisotropy, and in sec. 2.3 we consider two case
studies to show the typical signatures in proton velocity space at different
magnetic structures in the two types of wind. Results from numerical
simulations and a comparison between numerical outputs and PSP data are
reported in section 3. The summary and discussion are in section 4.

##  2 PSP observations

###  2.1 Overview of fluctuations’ properties for E6-E9

In this work we have used PSP data from Encounter 6 through Encounter 9
(E6-E9), which occurred from September 9, 2020, to August 15, 2021 covering a
range of radial distances from the sun $0.07<R<0.263$ AU. We use magnetic
field data from the flux-gate magnetometer at 4 samples/cycle resolution
(FIELDS; Bale et al. (2016)). The 3D proton velocity distribution functions
(L2) and its moments (L3) at 3s resolution are obtained by the electrostatic
analyzer (SPAN-I; Livi et al. (2022)), part of the Solar Wind Electrons Alphas
and Protons instrument suite (SWEAP; Kasper et al. (2016)); Since SPAN-I is
partially obstructed by PSP’s thermal protection shield, we use the quasi-
thermal noise (QTN) measured by the FIELDS Radio Frequency Spectrometer
(Moncuquet et al. (2020)) to have a more accurate determination of plasma
density. The total parallel and perpendicular proton temperatures, $T_{\bot}$
and $T_{\parallel}$, are obtained by projecting the temperature tensor in the
parallel and perpendicular directions with respect to the magnetic field.

Figure 1: PSP Encounter 6 data: spacecraft radial distance $R$ (top panel);
proton and electron number density $n_{i,e}$ (second panel); magnetic field
${\bf b}$, its magnitude $B$ and proton bulk velocity ${\bf V}$ in the
instrument frame (third and fourth panels); normalized cross-helicity
$\sigma_{c}$, residual energy $\sigma_{r}$ and the cosine of the angle between
the velocity and magnetic field fluctuations $\cos\theta_{vb}$ (bottom panel).

For reference, in Fig. 1 we provide an overview of E6. The top four panels
show the spacecraft radial distance ($R$), the proton and electron number
density ($n_{i,e}$), and the magnetic field (${\bf b}$ and $B=|{\bf b}|$) and
proton bulk velocity (${\bf V}$) in the instrument frame ($x$-component
pointing toward the sun). The fifth panel shows the Alfvénic properties of the
wind represented by the normalized cross-helicity $\sigma_{c}$ and residual
energy $\sigma_{r}$ of fluctuations, given respectively by

$\sigma_{c}=\frac{\langle|{\bf z}^{+}|^{2}\rangle-\langle|{\bf
z}^{-}|^{2}\rangle}{\langle|{\bf z}^{+}|^{2}\rangle+\langle|{\bf
z}^{-}|^{2}\rangle}\quad\text{and}\quad\sigma_{r}=\frac{\langle|\mathbf{\delta
V}|^{2}\rangle-\langle|\mathbf{\delta
b}|^{2}\rangle/(\mu_{0}\rho)}{\langle|\mathbf{\delta
V}|^{2}\rangle+\langle|\mathbf{\delta b}|^{2}\rangle/(\mu_{0}\rho)},$ (1)

where brackets $\langle...\rangle$ denote the moving average over a time
window, ${\bf z}^{\pm}=\mathbf{\delta V}\pm\mathbf{\delta
b}/\sqrt{\mu_{0}\rho}$ are the Elsässer variables defined by the velocity and
magnetic field fluctuations with respect to the mean field and by the QTN mass
density $\rho$ assuming quasi-neutrality. In this work we fix the averaging
window to 2 hours, that corresponds to several times the correlation time of
the magnetic field in the inner heliosphere (Parashar et al., 2020; Sioulas et
al., 2022b). The cross-helicity measures the relative dominance of inward and
outward propagating Alfvénic fluctuations (${\bf z}^{\pm}$), and the residual
energy quantifies the partition between the fluctuations’ kinetic and magnetic
energy. The bottom panel of Fig. 1 shows the cosine of the angle $\theta_{vb}$
between the velocity and magnetic field, defined as

$\cos{\theta_{vb}}=\frac{\mathbf{\delta V}\cdot\mathbf{\delta
b}}{|\mathbf{\delta V}||\mathbf{\delta b}|}.$ (2)

As will be discussed below, fluctuations during E6 are remarkably Alfvénic,
with high cross-helicity intervals (defined by $|\sigma_{c}|>0.75$) that last
several days (particularly during and after perihelion). In addition, a
shorter interval of low-Alfvénicity is observed in the vicinity of the
heliospheric current sheet (September 25th and 26th, 2020), a complex
structure with dominant magnetic energy and nearly zero cross-helicity. This
agrees with the observed non-Alfvénic wind reported near the Heliospheric
Current Sheet (HCS) in E4 (Chen et al., 2021).

Figure 2: Fluctuations’ properties for E6-E9. Top panels: PDF of normalized
cross-helicity $\sigma_{c}$(top-left panel), residual energy $\sigma_{r}$
(middle panel) and the cosine angle of velocity and magnetic field
fluctuations $\cos\theta_{vb}$ (top-right panel). Different Encounters
correspond to a different line style as indicated in the legend at the top.
The shaded blue and yellow areas indicate values $|\sigma_{c}|>0.75$ and
$|\sigma_{c}|<0.25$, respectively. Bottom panels: PDF of the normalized
fluctuations amplitude $\delta b/|\langle{\bf b}\rangle|$ (bottom left) and of
the ratio of magnetic field and magnetic field strength variances $C_{b}$
(bottom right) for the Alfvénic wind (blue lines, corresponding to the blue
shaded area in the top left panel) and for the non-Alfvénic wind (yellow
lines, corresponding to the yellow shaded area). The values in the legend of
the bottom panels shows the mean value of each distribution.

In Fig. 2 we show the statistical analysis of Alfvénic properties and
compressibility of solar wind fluctuations from E6 through E9. The top panels
show the Probability Density Function (PDF) of $\sigma_{c},\sigma_{r}$ and
$\cos{\theta_{ub}}$ for each Encounter. Although data display a range of
values $-1\leq\sigma_{c}\leq 1$, intervals of imbalanced fluctuations with
$|\sigma_{c}|>0.5$ can be identified, with peaks of the distribution of
$\sigma_{c}$ at $|\sigma_{c}|\geq 0.75$ (see also Fig. 1). In general,
fluctuations are magnetically dominated, a trend which is consistent
throughout all Encounters (Chen et al., 2020; Shi et al., 2021), with mean
value $\sigma_{r}\approx-0.5$ and with only a few intervals being
characterized by an excess of kinetic energy. Lastly, the PDF of
$\cos{\theta_{ub}}$ peaks at $\pm 1$, corresponding to fluctuations with
aligned ${\bf V}$ and ${\bf B}$ in all the Encounters. The Alfvénic wind
(those intervals with $|\sigma_{c}|\geq 0.75$) present strong alignment, i.e.,
$|\cos\theta_{vb}|\sim 0.89$, while the non-Alfvénic wind (those intervals
with $|\sigma_{c}|\leq 0.25$) presents a moderate alignment with a broader
distribution of angles having mean value $|\cos\theta_{vb}|\sim 0.57$. In the
latter case, such alignment corresponds with the relaxation and suppression of
nonlinearities in standard balanced turbulence that has been previously
observed in both simulations and solar wind measurement (Matthaeus et al.,
2008, 2012).

The bottom panels of Fig. 2 show the PDF of the amplitude of magnetic field
fluctuations,

$\frac{\delta b}{|\langle{\bf
b}\rangle|}\coloneqq\sqrt{\frac{\langle(b_{x}-\langle
b_{x}\rangle)^{2}+(b_{y}-\langle b_{y}\rangle)^{2}+(b_{z}-\langle
b_{z}\rangle)^{2}\rangle}{\langle b_{x}\rangle^{2}+\langle
b_{y}\rangle^{2}+\langle b_{z}\rangle^{2}}},$ (3)

and of the ratio of the variances of the magnetic field and magnetic field
strength (Villante & Vellante, 1982),

$C_{b}\coloneqq\sqrt{\frac{\langle(B-\langle
B\rangle)^{2}\rangle}{\langle(b_{x}-\langle b_{x}\rangle)^{2}+(b_{y}-\langle
b_{y}\rangle)^{2}+(b_{z}-\langle b_{z}\rangle)^{2}\rangle}},$ (4)

the latter providing a proxy for magnetic compressibility. In this analysis,
we have sorted the solar wind into Alfvénic (intervals with
$|\sigma_{c}|>0.75$, corresponding to the blue-shaded areas in Fig. 2, top
left) and non-Alfvénic wind (intervals with $|\sigma_{c}|<0.25$, corresponding
to the yellow-shaded area). For the non-Alfvénic wind we have also removed
data corresponding to the HCS crossing, where $|\langle{\bf b}\rangle|\approx
0$ (for the time window considered of 2-hours), to avoid fictitious long tails
of the PDF.

From Fig. 2, bottom left panel, it can be observed that the distributions of
the normalized fluctuations’ amplitude are remarkably different between
Alfvénic (blue) and non Alfvénic (yellow) wind. In the former case, relative
fluctuations’ amplitudes are bounded roughly by $\delta b/|\langle{\bf
b}\rangle|\lesssim 2$. This is in line with previous observations of Alfvénic
fluctuations in Helios and Ulysses, where saturation of amplitudes (although
defined differently) was found, a feature that is consistent with spherically
polarized fluctuations (Matteini et al., 2018). In the non-Alfvénic case,
instead, relative fluctuations’ amplitudes are large and there is no
constraint on their value (E9 being the only ambiguous case). We interpret the
different PDF of the normalized amplitude as a signature of a higher level of
compressibility in the non-Alfvénic wind (at the timescale of 2 hours). The
bottom right panel of Fig. 2 shows the distribution of $C_{b}$ for both types
of wind. As expected, $C_{b}$ in the non-Alfvénic wind is systematically
larger than in the Alfvénic wind, with the PDF mean in the non-Alfvénic wind
reaching up to 5 times that of the Alfvénic wind.

To summarize, we have considered the statistical properties of fluctuations
over a timescale of 2 hours. From this part of our data analysis, we conclude
that fluctuations in the Alfvénic ($|\sigma_{c}|>0.75$) and non-Alfvénic
($|\sigma_{c}|<0.25$) wind at distances $R<0.25$ AU display statistical
properties similar to those traditionally observed further away from the sun,
namely, almost incompressible and saturated fluctuations in the Alfvénic wind
and highly compressible fluctuations in the non-Alfvénic wind.

###  2.2 Correlations between PVI and proton thermal properties in different
types of wind

To study the connection between proton heating and discontinuities or, more
generally, intermittent structures in different plasma conditions, we utilized
the Partial Variance of Increments (PVI) method by calculating the following
quantity (Greco et al., 2008b):

$PVI=\frac{|\Delta\mathbf{B}(t,\tau)|}{\sqrt{\langle|\Delta\mathbf{B}(t,\tau)|^{2}\rangle}},$
(5)

where $\Delta\mathbf{B}(t,\tau)=\ \mathbf{b}(t+\tau)-\mathbf{b}(t)$ is the
magnetic field increment vector. The PVI method identifies non-Gaussian
features in the magnetic field and events with values $PVI>3$ are typically
associated with coherent structures or discontinuities. Although the PVI
defined in Eq. 5 captures variations in the magnetic field, it does not
distinguish RDs, TDs or switchbacks, that sometimes display properties of both
RDs and TDs (Larosa et al., 2021). For this reason, we also considered the PVI
of the magnetic field strength, hereafter called $MAG$-$PVI$, by taking the
variation of $B$, i.e., by substituting $\Delta{\bf
B}(t,\tau)\Rightarrow\Delta B(t,\tau)=B(t+\tau)-B(t)$ in Eq. (5). A comparison
of events identified with $PVI$ and $MAG$-$PVI$ allows us to gain insights on
whether compressible or incompressible magnetic structures are statistically
associated to local particle energization. Finally, because of the resolution
of SPAN-i data, we used a time lag of $\tau=12\ s$, that lies within the range
of scales where the transition from the inertial to the kinetic regime takes
place, to correlate discontinuities at kinetic scales and proton heating (or
variations of plasma moments).

After constructing the time series of $PVI$ (and $MAG$-$PVI$), we computed the
normalized probability distribution PDF of various functions of temperature
$f(T)$ conditioned on a given range of $PVI$ (and $MAG$-$PVI$) values

$\mathcal{PDF}(f(T)|\theta_{1}\leq PVI\leq\theta_{2}).$ (6)

We have calculated the conditional PDF of the $(i)$ total proton temperature
($f(T)=T$), $(ii)$ variation of temperature over the time lag $\tau=12\,s$
($f(T)=T(t+\tau)-T(t)\equiv\delta T$), $(iii)$ parallel and perpendicular
temperature ($f(T)=T_{\parallel},\,T_{\perp}$ where $\parallel$ and $\bot$ are
defined with respect to the local magnetic field direction
$\mathbf{\hat{b}}=\mathbf{b}/B$), and $(iv)$ temperature anisotropy
($f(T)=T_{\perp}/T_{\parallel}$). We performed the same PVI analysis with
$\tau=1\ s$ and then resampled the PVI data to match the SPAN-i resolution,
but we did not find significant quantitative differences.

Figure 3: Conditional PDF with respect to $PVI$ (top panels) and $MAG$-$PVI$
(bottom panels) of proton temperatures using data from E6 through E9: PDF of
parallel proton temperature $T_{\parallel}$ (left panels); PDF of
perpendicular proton temperature $T_{\bot}$ (middle panels); PDF of proton
temperature anisotropy $T_{\bot}/T_{\parallel}$ (right panels). Solid lines
correspond to Alfvénic wind and dashed lines to non-Alvénic wind. Different
color lines show PDF that corresponds to different ranges of $PVI$
($MAG$-$PVI$): black lines ($0<PVI\leq 1.5$), blue lines ($1.5<PVI\leq 3.0$),
and red lines ($PVI>3.0$). The numbers on the legends indicate the mean value
of each distribution.

In agreement with other studies, we find that the highest values of $T$ and
$\delta T$ are associated with the highest PVI values in both winds (not shown
here; see also, e.g., Osman et al. (2010); Qudsi et al. (2020); Sioulas et al.
(2022b)), supporting the idea that heating occurs at localized structures.
Results of our PVI analysis are reported in Fig. 3 for parallel and
perpendicular temperatures.

Fig. 3 shows the conditional PDF with respect to $PVI$ (top panels) and
$MAG$-$PVI$ (bottom panels) of $T_{\parallel}$ (left), $T_{\bot}$ (middle) and
$T_{\bot}/T_{\parallel}$ (right). The PDFs have been obtained by combining
data from all Encounters (E6-E9) and by using the cross helicity to sort
Alfvénic ($|\sigma_{c}|>0.75$, solid lines) and non-Alfvénic
($|\sigma_{c}|<0.25$, dashed lines) intervals. Red, blue and black colors
correspond to different ranges of $PVI$ ($MAG$-$PVI$), that we split into
$(0<PVI\leq 1.5)$, $(1.5<PVI\leq 3.0)$, and $(PVI>3.0)$ respectively. The
legend in each panel reports the mean value of each distribution. Finally, E6
was also analyzed individually, but we found the same trends as those shown in
Fig. 3 and discussed below.

As can be seen from the PDFs of $T_{\bot}/T_{\parallel}$ (Fig. 3, right
panels), the non-Alfvénic wind is more anisotropic than Alfvénic wind and, on
average,
$(T_{\bot}/T_{\parallel})_{nonAlfv}<(T_{\bot}/T_{\parallel})_{Alfv}<1$. As we
will see, such an anisotropy is due to the ubiquitous presence of proton
beams. While proton beams have been commonly observed at all radial distances
in the Alfvénic wind before PSP (e.g., Marsch (2012)), the frequent
observation of field-aligned beams in the non-Alfvénic wind (including the so-
called “hammerhead” distributions; Verniero et al. (2020, 2022)) is one of the
unexpected results of PSP. Further inspection of the temperature anisotropy
PDFs shows that Alfvénic and non-Alfvénic wind display an opposite correlation
between $PVI$ (and $MAG$-$PVI$) values, and $T_{\bot}/T_{\parallel}$. For the
$PVI$ case (Fig. 3, top right panel), the mean values of the non-Alfvénic wind
anisotropy lie in the range $0.71\lesssim T_{\perp}/T_{\parallel}\lesssim
0.73$ with increasing $PVI$ (dashed lines). In the Alfvénic wind (solid
lines), we find instead $0.88\lesssim T_{\perp}/T_{\parallel}\lesssim 0.97$
with decreasing $PVI$. The $MAG$-$PVI$ case (Fig. 3, bottom right panel) shows
similar values of anisotropy and trends.

Insights on parallel and perpendicular temperature can be found by cross-
checking PDFs of $T_{\bot}/T_{\parallel}$ with the PDFs of $T_{\parallel}$ and
$T_{\bot}$ (Fig. 3, left and middle panels). Interestingly, we find that the
ensembles of events selected through each $MAG$-$PVI$ value have consistently
higher average temperatures than the ensembles of events selected with the
traditional $PVI$ values. With increasing $MAG$-$PVI$ values, the average
temperatures are in the range $(4.792<T_{\parallel}<5.768)\times 10^{5}$ K and
$(4.601<T_{\bot}<5.405)\times 10^{5}$ K in Alfvénic wind. In the non-Alfvénic
wind we find $(4.097<T_{\parallel}<4.596)\times 10^{5}$ K and
$(2.582<T_{\bot}<3.124)\times 10^{5}$. By contrast, the largest temperatures
found with $PVI>3$ are $T_{\parallel}=5.358\times 10^{5}$ K and
$T_{\bot}=4.569\times 10^{5}$ K (Alfvénic), and $T_{\parallel}=4.125\times
10^{5}$ K and $T_{\bot}=2.785\times 10^{5}$ K (non Alfvénic).

To quantify the trends of perpendicular and parallel temperature with
$MAG$-$PVI$, we calculated the relative temperature enhancement between the
smallest and the largest range of $MAG$-$PVI$ considered. In the Alfvénic wind
we find $\Delta T_{\parallel}/T_{\parallel}\simeq 20\%$ and $\Delta
T_{\bot}/T_{\bot}\simeq 17\%$, indicating that the relative enhancements of
parallel temperature are larger than those in the perpendicular temperature as
smaller scales are considered. The opposite trend is found in the non-Alfvénic
wind, where $\Delta T_{\parallel}/T_{\parallel}\simeq 12\%$ and $\Delta
T_{\bot}/T_{\bot}\simeq 21\%$. Such an opposite trend of temperature
enhancements explains why the anisotropy $T_{\bot}/T_{\parallel}$ shifts away
from (towards) unity as PVI values increase for the Alfvénic (non-Alfvénic)
wind.

From this analysis we conclude, first, that structures with the largest
variations in $B$ (magnetic compressible structures) contribute to the highest
$T_{\parallel}$ and $T_{\bot}$ in both types of wind, suggesting that the
largest temperature enhancements occur in compressible structures also in the
Alfvénic wind. Second, the Alfvénic wind shows a preferential enhancement of
$T_{\parallel}$ as smaller scale structures are considered, whereas the non-
Alfvénic wind shows a preferential enhancement of $T_{\bot}$.

###  2.3 Proton kinetic features at different magnetic structures

Figure 4: Data from two intervals during E6. The non-Alfvénic interval is
shown on the left panels and the Alfvénic one on the right. From top to
bottom: proton mass density ($\rho\ \ [kg/m^{3}]$); magnetic field (${\bf b}$
and $B\ [nT]$); proton bulk velocity (${\bf V}\ \ [km/s^{-1}]$); parallel and
perpendicular temperature ($T_{\parallel}$ and $T_{\bot}\ \ [eV]$);
temperature anisotropy ($T_{\bot}/T_{\parallel}$) and plasma beta ($\beta$).
The last two panels show fluctuations properties ($\sigma_{r}$, $\sigma_{c}$
and $\cos\theta_{vb}$) and the PVI values.

To investigate the local energization of protons at
discontinuities/intermittent structures and kinetic features of the
distribution function in velocity space, we manually inspected 1-hour
intervals during quiet solar wind and without switchbacks. We selected two
sub-intervals that characterize the Alfvénic and the non-Alfvénic wind
(according to high and low cross-helicity), respectively. The two selected
sub-intervals occurred during E6 at approximately the same radial distance
from the sun ($30$ solar radii). Different fields and particle quantities are
reported in Fig. 4, where panels on the left correspond to the non-Alfvénic
wind and those on the right to the Alfvénic wind.

The non-Alfvénic sub-interval occurred near the heliospheric current sheet
crossing, on September 25th, 2020 from 03:10:00 to 04:10:00 UT, with an
average cross-helicity $\sigma_{c}\simeq-0.25$. This interval presents many
discontinuities with large PVI values. The structures are substantially
compressible with variations of $B$ and $\beta$ enhancements.

Figure 5: The reduced proton VDF as a function of parallel
($f(t,v_{\parallel})$) and perpendicular ($f(t,v_{\bot_{1,2}})$) velocity
components for the non-Alfvénic (left column) and Alfvénic (right column) sub-
intervals. The bottom panels show two examples of the reconstructed 3D VDF
$f(v_{\parallel},v_{\bot_{1}},v_{\bot,2})$ and its 2D projections,
$f(v_{\parallel},v_{\bot_{1,2}})$, at the time of an event with $PVI>3$ in
each type of wind.

The Alfvénic interval occurred on September 30th, 2020 between 01:00:00 to
02:00:00 UT with an overall
$\mathinner{\\!\left\lvert\sigma_{c}\right\rvert}\simeq 0.88$. The
fluctuations show the typical highly-correlated velocity-magnetic field that
characterizes wave-like fluctuations with rotation of $\cos\theta_{ub}$ near
small-scale structures. In general, the variations of the magnitude of the
magnetic field and the PVI values are smaller compared to the discontinuities
in the non-Alfvénic wind. In both cases, however, the most substantial
variations in proton temperature occur within large $PVI$ and $MAG$-$PVI$
structures and correspond to a net enhancement of $T_{\parallel}/T_{\bot}$, as
can be seen by comparing the fourth or fifth panel and the bottom panel.

The changes in temperature at these small timescales are connected to the
local particle heating processes at magnetic structures. The signatures in
proton VDF (velocity distribution function) of local heating for these two
representative sub-intervals are reported in Fig. 5, where we show the reduced
VDF ($f(t,v_{\parallel}),f(t,v_{\bot_{1,2}})$ for the non-Alfvénic (left
panel) and for the Alfvénic (right panels) sub-intervals described above (and
reported in Fig. 4). We also show the reconstructed 3D VDFs in the plane
$\\{(v_{\parallel},v_{\bot_{1}},v_{\bot_{2}})\\}$, and its projection into the
$(v_{\parallel},v_{\bot_{1,2}})$ planes, at the bottom panels for each type of
wind. The reduced VDF has been obtained by first transforming the initial 3D
energy distribution function from the $\\{E,\theta,\varphi\\}$ space (energy,
elevation, and azimuth) to velocity space in the field-aligned coordinate
system $\\{v_{\parallel},v_{\bot_{1}},v_{\bot_{2}}\\}$. Details can be found
in Appendix A.

As can be seen from Fig. 5, the proton VDF shows strong deviations from local
thermal equilibrium and gyrotropy in the vicinity of large PVI values. We
observe that protons have enhanced perpendicular velocities in those regions,
and parallel propagating proton beams around the local Alfvén speed are
observed near discontinuities in both Alfvénic and non-Alfvénic intervals. We
note some differences in those two distributions. In the non-Alfvénic case the
beam is hot in the perpendicular direction. Such proton beams with large
$T_{\bot,b}/T_{\bot,c}$ are consistent with the “hammerhead” distribution
(Verniero et al., 2022). In the Alfvénic case the beam is instead cooler and
more focused in the parallel direction. This can be seen by inspecting Fig. 6,
which shows a scatter plot of the beam-to-core drift velocity normalized by
Alfvén speed versus $(T_{\bot,b}/T_{\bot,c})^{*}$, which is the ratio of the
perpendicular temperature of the beam and the core
($T_{\bot,{b}}/T_{\bot,{c}}$) normalized by the combined beam+core
$T_{\bot,{b+c}}/T_{\parallel,{b+c}}$ values. Details on the fitting method can
be found in the appendix B. For reference, the crosses show data point at $t=$
03:35:39 UT (non-Alfvénic wind) where $(T_{\bot,b}/T_{\bot,c})^{*}=2.81$ and
at $t=$ 01:22:33 UT, where $(T_{\bot,b}/T_{\bot,c})^{*}=1.16$. To help
intuition, we have split the scatter plot into four quadrants that we use as a
reference to characterize the shape of the VDF. On the right side of the
diagram, the values on the lower quadrant represent VDFs with a focused beam,
and values on the upper quadrant represent VDFs with hotter beams or
“hammerhead” distributions. The values on the left side, the beam and core are
not well separated, and one can arbitrarily attribute them to distributions
with no beams. A sketch of the shape of the distribution as one moves above
and below the line $({T_{\bot_{b}}/T_{\bot_{c}}})^{*}=1$ is also presented on
the bottom right corner of Fig. 6. Fig. 5 suggests that beams in the Alfvénic
and non Alfvénic wind have different origins and display distinct properties.

Figure 6: Scatter plot of the beam-to-core drift velocity normalized by Alfvén
speed versus the ratio of the perpendicular temperature of the beam and the
core for the same Alfvenic and non-Alfvenic intervals presented before. The
black crosses in the blue and red contours correspond to the local values at
$t=$ 03:35:39 UT and $t=$ 01:22:33 UT, respectively.

##  3 Hybrid simulations

###  3.1 Initial conditions

We performed hybrid-kinetic simulations to investigate proton heating at
different types of magnetic structures under different plasma conditions, and
compared our results with observations. We considered a quasi-neutral hybrid
plasma model with massless and isothermal electrons, while protons are modeled
as kinetic particles using the particle-in-cell method (Matthews, 1994; Franci
et al., 2018). We adopted the following normalization: lengths are normalized
to the proton inertial length $d_{i}=c/\omega_{p}$ with $\omega_{p}=(4\pi
ne^{2}/m_{i})^{1/2}$ the proton plasma frequency and the time is normalized to
the inverse of the proton gyrofrequency
$\omega_{ci}^{-1}=(eB_{0}/m_{i}c)^{-1}$ and velocities to the Alfvén speed
$v_{A}=B_{0}/(4\pi n_{0}m_{i})^{1/2}$ with $B_{0}$ the mean magnetic field.
The plasma beta for both ions and electrons is defined as $\beta_{p,e}=8\pi
n_{0}T_{p,e}/B_{0}^{2}$. To avoid energy accumulation at the grid scale, we
have included explicit resistivity with a corresponding dissipation length
($l_{d}$) related to the Reynolds number and the box size ($L$) through
$R_{e}\sim(L/l_{d})^{4/3}$, which is chosen to be greater than the grid size.

We have considered two different setups in 2.5D geometry characterized by two
different types of turbulent fluctuations resembling non-Alfvénic and Alfvénic
wind, namely, $(i)$ a decaying 2D balanced turbulence simulation, and $(ii)$ a
wave simulation with an initial parallel propagating, Alfvénic broadband
spectrum. For case $(i)$, hereafter referred to as “turbulence simulation”, we
considered an initial large-amplitude fluctuation with $\sigma_{c}\approx 0$,
in energy equipartition and with an out-of-plane guide field
$\mathbf{B_{0}}=B_{0}\mathbf{\hat{z}}$. The initial condition consists of
large-amplitude perturbations of perpendicular wave modes in the range
$0.196<k_{\bot}d_{i}<0.49$ with random phases, similar to previous work
(Servidio et al., 2012; Franci et al., 2015; Cerri & Califano, 2017). For case
$(ii)$, referred to as “wave-like simulation”, we considered an in-plane mean-
field $\mathbf{B_{0}}=B_{0}\mathbf{\hat{x}}$ and the initial perturbation
satisfies the Walèn relation in the dispersionless limit $\delta{\bf
u}=-(\omega_{0}/k_{0})\delta{\bf b}$. The wave frequency $\omega_{0}$ is
determined from the normalized dispersion relation
$k_{0}^{2}=\omega_{0}^{2}/(1-\omega_{0})$ for left-handed polarized parallel
propagation waves. This initial condition corresponds to a broadband Alfvénic
fluctuation comprised of parallel modes in the range
$0.196<k_{\parallel}d_{i}<0.49$ (see also Malara et al. (2000); Matteini et
al. (2010); González et al. (2021)). In both cases, an initially isotropic and
Maxwellian plasma with proton beta $\beta_{i}=0.5$ and equal ion and electron
temperature $T_{i}/T_{e}=1$ are considered. The guide field is $B_{0}=1$ and
the same root-mean-square (rms) of the magnetic fluctuations $\delta
b_{rms}/B_{0}=0.63$ is chosen for both simulations. We adopt periodic boundary
conditions and fix a box of side $L=128\ d_{i}$ by using $1024^{2}$ grid
points with mesh size $0.125\ d_{i}$ and $8000$ particles per cell.

Because of the different geometry adopted, interactions between fields and
particles are fundamentally different in the two setups (Gary et al., 2020).
Naturally, our simulations are not meant to represent realistically the two
types of wind (both requiring 3D fields), but rather to isolate different
processes that might be dominant in each type of wind and that bear specific
signatures in phase space. In section 3.2 we describe the evolution of fields
and particles in the two setups. In section 3.3 we discuss the correlation
between anisotropy and PVI, and the signatures in phase space of proton
heating/acceleration in correspondence with at large PVI values.

###  3.2 Overview of numerical results

Figure 7: Contour plot of the out-of-plane current density $J_{z}$ for the
turbulence simulation (left) and the wave-like simulation (right). The black
dot shows the location at which the time series shown in Fig. 8 have been
taken.

In the turbulence simulation, once the system reaches the fully developed
turbulent state, the energy stored in the fields is progressively converted
into thermal energy, resulting in an average (over the spatial domain)
preferential perpendicular heating with $\langle
T_{\bot}/T_{\parallel}\rangle_{box}=1.37$ at $t=120\omega_{ci}^{-1}$. The
wave-like simulation, by contrast, displays a strong enhancement of the
spatially averaged $T_{\parallel}$, and the temperature anisotropy reaches
$\langle T_{\bot}/T_{\parallel}\rangle_{box}=0.6$ at
$t=180\omega_{ci}^{-1}$111Additional plots and movies showing the global
dynamics and overall proton heating in the simulations can be found in the
supplemental material..

These different heating processes in the two setups are expected, since the
turbulence simulation leads to a well-developed magnetic energy spectrum
perpendicular (in $k$-space) to the guide field, allowing for more channels of
particle heating, such as stochastic heating and particle scattering at
current sheets (Cerri et al., 2021; Sioulas et al., 2022a). Those processes
are suppressed in the wave-like simulation due to the different geometry
adopted. On the other hand, in the wave-like simulation, the available energy
carried by the wave is converted into kinetic and thermal energy via the rapid
disruption of the wave packet mediated by wave steepening along the guide
field (i.e., through the formation of gradients parallel to ${\bf B}_{0}$);
this leads to both the formation of a field-aligned proton beam at about the
local Alfvén speed, and local perpendicular heating at the steepened edges
(González et al., 2021). This mechanism is in turn suppressed in the
turbulence simulation. Figure 7 shows the contour plot of $J_{z}$ for the two
setups to show the different types of structures that form nonlinearly.

Figure 8: Single-point time evolution of fields from the turbulence simulation
(left panels) and wave-like simulation (right panels). From top to bottom:
proton density, magnetic field, proton bulk velocity, electric field, parallel
and perpendicular temperature, temperature anisotropy, and plasma beta.
Finally, the last two panels show $\sigma_{c}$, $\sigma_{r}$ and
$\cos\theta_{vb}$, and the PVI values.

To make a comparison with PSP data, in Figure 8 we show the time series of
fields and plasma quantities in the two setups. The time series have been
taken at the location marked by the black dot in Figure 7. From top to bottom,
Figure 8 shows the time series of $n$, ${\bf B}$ and $|B|$, ${\bf u}$, the
electric field ${\bf e}$, $T_{\bot}$, $T_{\parallel}$ and $\beta$,
$\cos\theta_{vb}$, $\sigma_{r}$ and $\sigma_{c}$, and the PVI values (see
section 3.3 for details on the PVI analysis). The turbulence simulation (left
panels), is characterized, on average by $|\sigma_{c}|<0.5$ and $\sigma_{r}<0$
(predominantly). The wave-like simulation (right panel) has $\sigma_{r}\approx
0$ and $|\sigma_{c}|\approx 1$. As can be seen, temperature enhancements occur
in correspondence with large PVI values in both cases. In the turbulence case,
those large PVI values are associated with structures such as flux tubes,
current sheets and small-scale plasmoids. These structures have a strong
perpendicular electric field, and some also display large out-of-plane
electric and magnetic fields (along the z-axis) with strong compressibility,
mediating particle acceleration and heating (Dmitruk et al., 2004; Wan et al.,
2015; Comisso & Sironi, 2022). As discussed above, in the wave-like simulation
large PVI values are associated with the steepened edges of Alfvénic
fluctuations. Also, rapid and large variations of $\cos\theta_{ub}$ near these
structures can be identified (sixth panel), similar to what is observed during
the Alfvénic sub-interval presented in section 2.3. In contrast to the
magnetic structures in the turbulent simulation, the steepened edges are
mostly rotational discontinuities with an embedded compressive component
(Matteini et al., 2010; González et al., 2021).

###  3.3 Correlations between magnetic structures and proton kinetic features

The PDFs conditioned on $PVI$ and $PVI$-$MAG$ are shown in Fig. 9 in the same
format as Fig. 3. To obtain the PDFs, we used a similar methodology to that
employed for PSP data, but we have taken $\tau=1\omega_{ci}^{-1}$ since we
have enough resolution to inspect proton scales. We have then taken the time
series of plasma quantities at a resolution of $0.1\omega_{ci}^{-1}$ at
$10^{4}$ fixed equidistant points in the simulation domain.

Figure 9: Conditional PDF with respect to $PVI$ (top panels) and $MAG$-$PVI$
(bottom panels) of proton temperatures using simulation data: PDF of parallel
proton temperature $T_{\parallel}$ (in normalized units; left panels); PDF of
perpendicular proton temperature $T_{\bot}$ (middle panels); PDF of proton
temperature anisotropy $T_{\bot}/T_{\parallel}$ (right panels). The solid
lines correspond to the wave-like simulation and the dashed lines to the
turbulence simulation. Different color lines show PDF that corresponds to
different ranges of $PVI$ and $MAG$-$PVI$: Black lines ($0<PVI\leq 1.5$), blue
lines ($1.5<PVI\leq 3.0$), and red lines ($PVI>3.0$). The numbers on the
legends are the mean value of each distribution.

In general, there is a positive correlation between PVI values with $T$ and
$\delta T$ (not shown), indicating that the hottest protons are found locally
near small-scale structures (in this case represented by $PVI>3$ values).
Although the turbulence simulation is hotter than the wave-like simulation
(this is true both locally at discontinuities and globally over the simulation
domain), protons remain close to Maxwellian in the turbulence simulation. On
the contrary, in the wave-like simulation, particles develop large
anisotropies localized at small-scale structures. In that case such a local
anisotropy, reaching an average value of $T_{\bot}/T_{\parallel}\approx 2$ for
$PVI>3$, is caused by strong particle pitch-angle scattering at the steepened
fronts, resulting in effective local proton perpendicular energization
(González et al., 2021; Malara et al., 2021; Gonzalez et al., 2023). Since
this effect is localized at the steepened edges, the background temperature
(represented by the lower PVI range) and the spatially averaged temperatures
over the entire domain display an opposite anisotropy,
$T_{\bot}/T_{\parallel}<1$, as discussed previously, due to the formation of a
field-aligned beam that fills the entire simulation domain.

Figure 10: Single-point proton VDF from simulations. The reduced proton VDF of
perpendicular ($f(t,v_{\bot_{1},\bot_{2}},t)$) and parallel
$f(t,v_{\parallel})$ of proton velocity for the turbulence simulation (left
panels) and wave-like simulation (right panels). The bottom panels show the
reduced VDF $f(v_{\parallel},v_{\bot})$ at the time indicated by the red
dashed line.

Some differences between $PVI$ (Fig. 9, top panels) and $MAG$-$PVI$ (Fig. 9,
bottom panels) results can be noticed. In the wave-like simulation, the
magnetic structures correspond essentially to rotational discontinuities with
relatively small variations in the magnitude of $B$ and, thus, the correlation
with the various PVI ranges is reduced in the $MAG$-$PVI$ case. This can be
seen from the relative temperature variations between the smallest and the
largest range of PVIs, showing the aforementioned trend ($PVI$: $\Delta
T_{\parallel}/T_{\parallel}\simeq 26\%$ and $\Delta T_{\bot}/T_{\bot}\simeq
150\%$; $MAG$-$PVI$: $\Delta T_{\parallel}/T_{\parallel}\simeq 3\%$ and
$\Delta T_{\bot}/T_{\bot}\simeq 59\%$). For the turbulence simulation, on the
contrary, higher relative temperature increase is found in the $MAG-PVI$ case,
similarly to observations (PVI: $\Delta T_{\parallel}/T_{\parallel}\simeq
4.44\%$ and $\Delta T_{\bot}/T_{\bot}\simeq 4.47\%$; $MAG$-$PVI$: $\Delta
T_{\parallel}/T_{\parallel}\simeq 5.79\%$ and $\Delta T_{\bot}/T_{\bot}\simeq
8.33\%$).

In Fig. 10, top two panels, we show the gyroaveraged reduced perpendicular and
parallel VDF ($f(t,v_{\bot})$ and $f(t,v_{\parallel})$, respectively and, in
the bottom panels, we show the reduced VDF $f(v_{\parallel},v_{\bot})$ at the
time indicated by the dashed line, marking the occurrence of an event with
$PVI>3$ in each case. The left column corresponds to the turbulence simulation
and the right to the wave-like. As can be seen, magnetic
discontinuities/intermittent structures are characterized by non-Maxwellian
features, namely, temperature anisotropies and beams at the Alfvén speed.
However, the turbulence simulation undergoes perpendicular heating in both the
core and the beams, so that $T_{\perp}\gtrsim T_{\parallel}$ locally and on
average, and beams are “hot”. The wave-like simulation instead displays a
strong perpendicular heating of the core, localized at the discontinuities,
and a cold beam at the Alfvén speed (more focused in the perpendicular
direction than in the turbulence setup).

To summarize, our simulations show that, in analogy with PSP data, higher
temperatures and temperature variations are found at higher PVI values in both
setups (both winds), and that colder beams (as found in Alfvénic wind) are
associated with the steepening of Alfvénic fluctuations while hot beams (as
found in non-Alfvénic wind) are generated in low-cross helicity turbulence.
The development of a standard turbulent cascade also favors preferential
perpendicular heating (as in the non-Alfvénic wind). Even though our
simulations highlight the different origins and properties of proton beams in
the two types of setups/winds, and qualitatively reproduce observed properties
of non-Alfvénic wind, the do not reflect all the local trends of temperatures
at magnetic structures observed in the Alfvénic wind. In particular, the PVI-
temperature analysis in the wave-like simulation shows a much larger relative
increase of perpendicular temperature than observed (larger than the increase
in parallel temperature), and compressibility of the Alfvénic wind at small
sacels is also not well reproduced numerically. These discrepancies are likely
due to the geometry contraints. However, we also mention the possibility that
instrument resolution may also play a role in underestimating strong
temperature variations at small scales.

Before concluding we also remark that our simulations do not reproduce
differences between the two types of winds in their average properties, such
as the fact that the Alfvénic wind is hotter, and that the global (average)
anisotropy is $T_{\perp}/T_{\parallel}<1$. This is expected since our setups
cannot reproduce fully 3D dynamics. Furthermore, such differences might well
be related to different coronal origins and/or a different radial evolution of
Alfvénic and non-Alfvénic wind, an aspect that is not addressed in this work.

##  4 Summary and conclusions

We have studied the correlation between proton temperatures and magnetic
discontinuities/ intermittent structures in different solar wind turbulence
conditions (high and low cross helicity, i.e., Alfvénic and non-Alfvénic wind)
by using PSP observations from E6-E9 and hybrid-kinetic simulations. Our main
results and conclusions are summarized as follows:

  * •

At large (MHD) scales, the Alfvénic wind is much less compressible than the
non-Alfvénic wind (see Fig. 2). However, our $PVI$ and $MAG$-$PVI$ analyses
show that the hottest protons are localized at kinetic-scale compressible
structures in both types of wind (see Fig. 3).

  * •

There is a statistical correlation between the highest proton total
temperature and coherent structures (quantified by PVI values), consistent
with previous studies (Qudsi et al., 2020; Sioulas et al., 2022b).
Furthermore, the Alfvénic wind shows a preferential enhancement of
$T_{\parallel}$ as smaller scale structures are considered, whereas the non-
Alfvénic wind shows a preferential enhancement of $T_{\bot}$ (see Fig. 3).

  * •

Proton beams are ubiquitous in both types of wind (leading to an average
anisotropy $T_{\parallel}>T_{\bot}$.) However, local kinetic features of
proton VDFs differ in the two winds (see Fig. 5 and 6). The non-Alfvénic wind
is characterized by “hot” beams ($T_{\bot,b}/T_{\bot,c}^{*}\gtrsim 2$)
resembling the “hammerhead” distributions. The Alfvénic wind is characterized
by “colder” beams ($T_{\bot,b}/T_{\bot,c}^{*}\lesssim 1$).

  * •

We find some similarities between the hybrid-kinetic simulations and in-situ
measurements despite the limitations of the reduced geometry adopted. The
field aligned proton beams that develop in our simulations display distinct
features, supporting the idea that proton beams in Alfvénic and non-Alfvénic
wind have different properties and different origins. Simulations suggest that
the development of a perpendicular cascade, favored in balanced turbulence,
allows a preferential relative enhancement of $T_{\bot}$, and the formation of
hot beams via nonlinear dynamics and reconnection. On the contrary, cold
field-aligned beams are favored by Alfvén wave steepening (see Fig. 9 and Fig.
10).

Additionally, we have shown for the first time 3D proton VDFs from PSP
displaying non-Maxwellian and non-gyrotropic features near discontinuities,
warranting a more general approach to fit VDFs than the widely adopted fit to
bi-Maxwellians. Furthermore, non-Maxwellian and non-gyrotropic proton VDFs
around discontinuities/intermittent structures are found in both winds,
confirming that nonlinearities and strong deviations from non-thermal
distributions are intrinsically related in collisonless plasmas (Valentini et
al., 2014), resulting in a universal heating channel regardless of the
Alfvénic properties of the solar wind. In conclusion, this work contributes to
understanding the distinctive role of coherent structures in heating
collisionless plasmas. To gain a comprehensive understanding, our simulation
results should be extended to include three-dimensional effects and encompass
a broader range of initial Alfvénic properties rather than the extreme cases
considered here ($\sigma_{c}\simeq 0$ and $\sigma_{c}\simeq-1$) This will be
the subject of future investigations.

## Acknowledgements

This research was supported by NASA grant #80NSSC18K1211 and NSF CAREER award
2141564. We acknowledge the Parker Solar Probe (PSP) mission for the use of
the data publicly available at the NASA Space Physics Data Facility and the
Texas Advanced Computing Center (TACC) at The University of Texas at Austin
for providing HPC resources that have contributed to the research results
reported within this paper. Simulations have been run on the Frontera
supercomputer http://www.tacc.utexas.edu. JLV acknowledges support from NASA
PSP-GI 80NSSC23K0208 and NASA LWS 80NSSC22K1014.

## References

  * Artemyev et al. (2019) Artemyev, A., Angelopoulos, V., Vasko, I., et al. 2019, Geophysical Research Letters, 46, 1185
  * Bale et al. (2016) Bale, S., Goetz, K., Harvey, P., et al. 2016, Space science reviews, 204, 49
  * Bowen et al. (2022) Bowen, T. A., Chandran, B. D., Squire, J., et al. 2022, Physical Review Letters, 129, 165101
  * Burlaga (1991) Burlaga, L. 1991, Journal of Geophysical Research: Space Physics, 96, 5847
  * Cerri et al. (2021) Cerri, S. S., Arzamasskiy, L., & Kunz, M. W. 2021, The Astrophysical Journal, 916, 120
  * Cerri & Califano (2017) Cerri, S. S., & Califano, F. 2017, New Journal of Physics, 19, 025007
  * Chandran et al. (2010) Chandran, B. D. G., Li, B., Rogers, B. N., Quataert, E., & Germaschewski, K. 2010, ApJ, 720, 503, doi: 10.1088/0004-637X/720/1/503
  * Chen et al. (2020) Chen, C., Bale, S., Bonnell, J., et al. 2020, The Astrophysical Journal Supplement Series, 246, 53
  * Chen et al. (2021) Chen, C., Chandran, B., Woodham, L., et al. 2021, Astronomy & Astrophysics, 650, L3
  * Chen et al. (2001) Chen, L., Lin, Z., & White, R. 2001, Physics of Plasmas, 8, 4713
  * Colburn & Sonett (1966) Colburn, D., & Sonett, C. 1966, Space Science Reviews, 5, 439
  * Comisso & Sironi (2022) Comisso, L., & Sironi, L. 2022, The Astrophysical Journal Letters, 936, L27
  * Dmitruk et al. (2004) Dmitruk, P., Matthaeus, W., & Seenu, N. 2004, The Astrophysical Journal, 617, 667
  * D’Amicis et al. (2021) D’Amicis, R., Perrone, D., Bruno, R., & Velli, M. 2021, Journal of Geophysical Research: Space Physics, 126, e2020JA028996
  * Franci et al. (2018) Franci, L., Hellinger, P., Guarrasi, M., et al. 2018, Journal of Physics: Conference Series, 1031, 012002
  * Franci et al. (2015) Franci, L., Landi, S., Matteini, L., Verdini, A., & Hellinger, P. 2015, The Astrophysical Journal, 812, 21
  * Gary et al. (2020) Gary, S. P., Bandyopadhyay, R., Qudsi, R. A., et al. 2020, The Astrophysical Journal, 901, 160
  * Gary & Borovsky (2008) Gary, S. P., & Borovsky, J. E. 2008, Journal of Geophysical Research: Space Physics, 113
  * Gonzalez et al. (2023) Gonzalez, C., Innocenti, M. E., & Tenerani, A. 2023, arXiv preprint arXiv:2301.07646
  * González et al. (2021) González, C. A., Tenerani, A., Matteini, L., Hellinger, P., & Velli, M. 2021, The Astrophysical Journal Letters, 914, L36, doi: 10.3847/2041-8213/ac097b
  * Greco et al. (2008a) Greco, A., Chuychai, P., Matthaeus, W., Servidio, S., & Dmitruk, P. 2008a, Geophysical Research Letters, 35
  * Greco et al. (2008b) —. 2008b, Geophysical Research Letters, 35
  * Greco et al. (2009) Greco, A., Matthaeus, W. H., Servidio, S., Chuychai, P., & Dmitruk, P. 2009, The Astrophysical Journal, 691, L111, doi: 10.1088/0004-637x/691/2/l111
  * Greco et al. (2016) Greco, A., Perri, S., Servidio, S., Yordanova, E., & Veltri, P. 2016, The Astrophysical Journal Letters, 823, L39
  * Horbury et al. (2001) Horbury, T. S., Burgess, D., Fränz, M., & Owen, C. J. 2001, Geophysical research letters, 28, 677
  * Hudson (1970) Hudson, P. 1970, Planetary and Space Science, 18, 1611
  * Johnson & Cheng (2001) Johnson, J. R., & Cheng, C. 2001, Geophysical research letters, 28, 4421
  * Karimabadi et al. (2013) Karimabadi, H., Roytershteyn, V., Wan, M., et al. 2013, Physics of Plasmas, 20, 012303
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Science Reviews, 204, 131
  * Knetter et al. (2004) Knetter, T., Neubauer, F., Horbury, T., & Balogh, A. 2004, Journal of Geophysical Research: Space Physics, 109
  * Larosa et al. (2021) Larosa, A., Krasnoselskikh, V., de Wit, T. D., et al. 2021, Astronomy & Astrophysics, 650, A3
  * Lichko et al. (2017) Lichko, E., Egedal, J., Daughton, W., & Kasper, J. 2017, The Astrophysical Journal Letters, 850, L28
  * Livi et al. (2022) Livi, R., Larson, D. E., Kasper, J. C., et al. 2022, The Astrophysical Journal, 938, 138, doi: 10.3847/1538-4357/ac93f5
  * Malara et al. (2021) Malara, F., Perri, S., & Zimbardo, G. 2021, Physical Review E, 104, 025208
  * Malara et al. (2000) Malara, F., Primavera, L., & Veltri, P. 2000, Physics of Plasmas, 7, 2866, doi: 10.1063/1.874136
  * Marsch (2012) Marsch, E. 2012, Space science reviews, 172, 23
  * Matteini et al. (2010) Matteini, L., Landi, S., Velli, M., & Hellinger, P. 2010, Journal of Geophysical Research: Space Physics, 115
  * Matteini et al. (2018) Matteini, L., Stansby, D., Horbury, T., & Chen, C. H. 2018, The Astrophysical Journal Letters, 869, L32
  * Matthaeus et al. (2008) Matthaeus, W., Pouquet, A., Mininni, P. D., Dmitruk, P., & Breech, B. 2008, Physical review letters, 100, 085003
  * Matthaeus et al. (2012) Matthaeus, W. H., Montgomery, D. C., Wan, M., & Servidio, S. 2012, Journal of Turbulence, 13, N37
  * Matthews (1994) Matthews, A. P. 1994, Journal of Computational Physics, 112, 102
  * Medvedev et al. (1997) Medvedev, M. V., Diamond, P. H., Shevchenko, V. I., & Galinsky, V. L. 1997, Phys. Rev. Lett., 78, 4934, doi: 10.1103/PhysRevLett.78.4934
  * Moncuquet et al. (2020) Moncuquet, M., Meyer-Vernet, N., Issautier, K., et al. 2020, The Astrophysical Journal Supplement Series, 246, 44
  * Narita & Marsch (2015) Narita, Y., & Marsch, E. 2015, The Astrophysical Journal, 805, 24
  * Neugebauer (2006) Neugebauer, M. 2006, Journal of Geophysical Research: Space Physics, 111, doi: https://doi.org/10.1029/2005JA011497
  * Neugebauer et al. (1984) Neugebauer, M., Clay, D., Goldstein, B., Tsurutani, B., & Zwickl, R. 1984, Journal of Geophysical Research: Space Physics, 89, 5395
  * Osman et al. (2014) Osman, K., Matthaeus, W., Gosling, J., et al. 2014, Physical Review Letters, 112, 215002
  * Osman et al. (2010) Osman, K. T., Matthaeus, W. H., Greco, A., & Servidio, S. 2010, The Astrophysical Journal, 727, L11, doi: 10.1088/2041-8205/727/1/l11
  * Osman et al. (2012) Osman, K. T., Matthaeus, W. H., Wan, M., & Rappazzo, A. F. 2012, Phys. Rev. Lett., 108, 261102, doi: 10.1103/PhysRevLett.108.261102
  * Parashar et al. (2020) Parashar, T. N., Goldstein, M. L., Maruca, B. A., et al. 2020, The Astrophysical Journal Supplement Series, 246, 58, doi: 10.3847/1538-4365/ab64e6
  * Parker (1994) Parker, E. N. 1994, Spontaneous current sheets in magnetic fields: with applications to stellar x-rays, Vol. 1 (Oxford University Press)
  * Paschmann et al. (2013) Paschmann, G., Haaland, S., Sonnerup, B., & Knetter, T. 2013, Annales Geophysicae, 31, 871, doi: 10.5194/angeo-31-871-2013
  * Perrone et al. (2016) Perrone, D., Alexandrova, O., Mangeney, A., et al. 2016, The Astrophysical Journal, 826, 196
  * Perrone et al. (2017) Perrone, D., Alexandrova, O., Roberts, O., et al. 2017, The Astrophysical Journal, 849, 49
  * Phillips et al. (2023) Phillips, C., Bandyopadhyay, R., McComas, D. J., & Bale, S. D. 2023, Monthly Notices of the Royal Astronomical Society: Letters, 519, L1
  * Qudsi et al. (2020) Qudsi, R. A., Maruca, B. A., Matthaeus, W. H., et al. 2020, The Astrophysical Journal Supplement Series, 246, 46, doi: 10.3847/1538-4365/ab5c19
  * Sahraoui et al. (2010) Sahraoui, F., Goldstein, M. L., Belmont, G., Canu, P., & Rezeau, L. 2010, Physical review letters, 105, 131101
  * Servidio et al. (2011) Servidio, S., Greco, A., Matthaeus, W. H., Osman, K. T., & Dmitruk, P. 2011, Journal of Geophysical Research: Space Physics, 116, doi: https://doi.org/10.1029/2011JA016569
  * Servidio et al. (2012) Servidio, S., Valentini, F., Califano, F., & Veltri, P. 2012, Physical review letters, 108, 045001
  * Shi et al. (2021) Shi, C., Velli, M., Panasenco, O., et al. 2021, Astronomy & Astrophysics, 650, A21
  * Sioulas et al. (2022a) Sioulas, N., Isliker, H., & Vlahos, L. 2022a, Astronomy & Astrophysics, 657, A8
  * Sioulas et al. (2022b) Sioulas, N., Velli, M., Chhiber, R., et al. 2022b, The Astrophysical Journal, 927, 140
  * Tessein et al. (2013) Tessein, J. A., Matthaeus, W. H., Wan, M., et al. 2013, The Astrophysical Journal, 776, L8, doi: 10.1088/2041-8205/776/1/l8
  * Tsurutani et al. (1994) Tsurutani, B., Ho, C., Smith, E., et al. 1994, Geophysical Research Letters, 21, 2267
  * Tsurutani & Ho (1999) Tsurutani, B. T., & Ho, C. M. 1999, Reviews of Geophysics, 37, 517, doi: https://doi.org/10.1029/1999RG900010
  * Valentini et al. (2014) Valentini, F., Servidio, S., Perrone, D., et al. 2014, Physics of Plasmas, 21
  * Vasquez et al. (2007) Vasquez, B. J., Abramenko, V. I., Haggerty, D. K., & Smith, C. W. 2007, Journal of Geophysical Research: Space Physics, 112, doi: https://doi.org/10.1029/2007JA012504
  * Vasquez & Hollweg (2001) Vasquez, B. J., & Hollweg, J. V. 2001, Journal of Geophysical Research: Space Physics, 106, 5661, doi: https://doi.org/10.1029/2000JA000268
  * Vech et al. (2017) Vech, D., Klein, K. G., & Kasper, J. C. 2017, The Astrophysical Journal Letters, 850, L11
  * Verniero et al. (2022) Verniero, J., Chandran, B., Larson, D., et al. 2022, The Astrophysical Journal, 924, 112
  * Verniero et al. (2020) Verniero, J. L., Larson, D. E., Livi, R., et al. 2020, The Astrophysical Journal Supplement Series, 248, 5, doi: 10.3847/1538-4365/ab86af
  * Villante & Vellante (1982) Villante, U., & Vellante, M. 1982, Solar Physics, 81, 367
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wan et al. (2015) Wan, M., Matthaeus, W. H., Roytershteyn, V., et al. 2015, Phys. Rev. Lett., 114, 175002, doi: 10.1103/PhysRevLett.114.175002
  * Woodham et al. (2021) Woodham, L., Horbury, T., Matteini, L., et al. 2021, Astronomy & Astrophysics, 650, L1
  * Zhdankin et al. (2012) Zhdankin, V., Boldyrev, S., & Mason, J. 2012, The Astrophysical Journal, 760, L22, doi: 10.1088/2041-8205/760/2/l22

##  Appendix A VDF interpolation method

To visualize the proton velocity distribution in velocity space we mapped the
initial 3D energy distribution function from energy, elevation, and azimuthal
angle coordinates $(E=1/2mv^{2},\theta,\phi)$ into velocity-space coordinates
in the instrument frame ($v_{x},v_{y},v_{z}$) using transformation from
spherical to cartesian coordinates:

$v_{x}=v\cos\theta\cos\phi,\ \ \ \ \ v_{y}=v\cos\theta\sin\phi,\ \ \ \ \
v_{z}=v\sin\theta,$

with $v$ the amplitude of the velocity field. The $v_{z}$ component is defined
as the rotational symmetry axis of the instrument, $v_{x}$ direction points
toward the sun, and the $v_{y}$ component is then defined orthogonal to them.
However, the results presented in section 2 are shown in the field-aligned
frame, which was obtained by applying a rotation matrix to the proton velocity
vectors in the instrument frame:

$\begin{bmatrix}v_{\parallel}\\\ v_{\bot_{1}}\\\
v_{\bot_{2}}\end{bmatrix}=\begin{bmatrix}\cos\psi&-k_{z}\sin\psi&k_{y}\sin\psi\\\
k_{z}\sin\psi&k_{y}^{2}+k_{z}^{2}\cos\psi&k_{y}k_{z}(1-\cos\psi)\\\
-k_{y}\sin\psi&k_{y}k_{z}(1-\cos\psi)&k_{z}^{2}+k_{y}^{2}\cos\psi\end{bmatrix}\begin{bmatrix}v_{x}\\\
v_{y}\\\ v_{z}\end{bmatrix},$

The above expression is obtained by using the Euler-Rodriguez formula that
corresponds to a rotation by an angle $\psi$ about an axis defined by the unit
vector $\mathbf{\hat{k}}$. Here we chose this axis to be defined as
$\mathbf{\hat{k}}=[0,-\hat{b_{z}},\hat{b_{y}}]$ and $\psi$ the angle between
the instrument frame pointing to the sun and the magnetic field vector
$\mathbf{\hat{b}}$. Further details can be found in Woodham et al. (2021). To
obtain a reduced representation (1D/2D) of the integrated velocity
distribution function that accurately accounts for the different velocity
ranges at different planes, we used 3D linear interpolation on a mesh with a
fixed velocity range. This method ensures the proper integration along the
preferred component that appropriately adds up the number of counts in each
velocity bin. The interpolation method is commonly used to approximate a
function from a set of discrete data and it returns a callable function that
can be used to evaluate the interpolated function at any point within a
defined interval. We employed the cubic interpolation method from the griddata
function in the scipy.interpolate module (Virtanen et al., 2020). For
illustration, Figure A1 presents a 3D render of the original and the
interpolated proton VDF. The Python code can be found in the following
repository.

Figure A1: Illustration of the interpolation method. Several VDFs that
correspond to different measurements in azimuthal angles (represented by the
range of colors from yellow to purple) and for different elevation angles
(arranged from left to right). For each VDF, it is shown the original signal
(solid lines) and the interpolated signal (dashed lines).

##  Appendix B Fitting method

To obtain additional information on the proton VDFs, the SPAN-I L2 sf00 data
product was fit to a sum of two 3D bi-Maxwellians, one for the “core” and one
for “beam,” following a similar procedure as described in (Verniero et al.,
2020). Note that the “beam” was constrained to lie parallel to the mean
magnetic field and the ”core” was labelled as the peak in phase-space density.
Since SPAN-I is partially obstructed by PSP’s Thermal Protection Shield, only
partial moments of the VDF are obtained. We therefore only include fits that
are at least 80% in the field-of-view (FOV). We quantify this by first
computing moments of the distribution from a VDF reconstructed from the
fitting parameters. Next, we take the sum of the fitted beam and core
densities and divide by the reconstructed moment density. To exceed a FOV
threshold of 80%, this number must exceed 0.8.