# On the role of solar wind expansion as a source of whistler waves:
scattering of suprathermal electrons and heat flux regulation in the inner
heliosphere.

A. Micera Solar-Terrestrial Centre of Excellence - SIDC, Royal Observatory of
Belgium, Brussels, Belgium.  Centre for Mathematical Plasma Astrophysics, KU
Leuven, Leuven, Belgium.  A. N. Zhukov Solar-Terrestrial Centre of Excellence
- SIDC, Royal Observatory of Belgium, Brussels, Belgium.  Skobeltsyn Institute
of Nuclear Physics, Moscow State University, Moscow, Russia.  R. A. López
Departamento de Física, Universidad de Santiago de Chile, Santiago, Chile  E.
Boella Physics Department, Lancaster University, Lancaster, UK.  Cockcroft
Institute, Daresbury Laboratory, Warrington, UK.  A. Tenerani Department of
Physics, The University of Texas at Austin, Austin, USA.  M. Velli University
of California Los Angeles, Department of Earth, Planetary, and Space Sciences,
Los Angeles, USA.  G. Lapenta Centre for Mathematical Plasma Astrophysics, KU
Leuven, Leuven, Belgium.  M. E. Innocenti Institut für Theoretische Physik,
Ruhr-Universität Bochum, Bochum, Germany.

###### Abstract

The role of solar wind expansion in generating whistler waves is investigated
using the EB-iPic3D code, which models solar wind expansion self-consistently
within a fully kinetic semi-implicit approach. The simulation is initialized
with an electron velocity distribution function modeled after Parker Solar
Probe observations during its first perihelion at 0.166 au, consisting of a
dense core and an anti-sunward strahl. This distribution function is initially
stable with respect to kinetic instabilities. Expansion drives the solar wind
into successive regimes where whistler heat flux instabilities are triggered.
These instabilities produce sunward whistler waves initially characterized by
predominantly oblique propagation with respect to the interplanetary magnetic
field. The excited waves interact with the electrons via resonant scattering
processes. As a consequence, the strahl pitch angle distribution broadens and
its drift velocity reduces. Strahl electrons are scattered in the direction
perpendicular to the magnetic field, and an electron halo is formed. At a
later stage, resonant electron firehose instability is triggered and further
affects the electron temperature anisotropy as the solar wind expands. Wave-
particle interaction processes are accompanied by a substantial reduction of
the solar wind heat flux. The simulated whistler waves are in qualitative
agreement with observations in terms of wave frequencies, amplitudes and
propagation angles. Our work proposes an explanation for the observations of
oblique and parallel whistler waves in the solar wind. We conclude that solar
wind expansion has to be factored in when trying to explain kinetic processes
at different heliocentric distances.

Plasma Astrophysics — Solar wind — Space plasmas

##  1 Introduction

Collisionless effects play a major role in the energy balance of the solar
wind plasma (e.g. Marsch, 2006). This is clear from the complex non-thermal
shapes of the solar wind particle Velocity Distribution Functions (VDFs). The
measurement of VDFs close to the Sun is crucial for two reasons. First, VDFs
may carry signatures of the mechanisms responsible for the solar wind heating
and acceleration at the early stages of its evolution (Ko et al., 1996; Berčič
et al., 2020). Second, information on the VDFs may be key to understanding the
processes affecting the further evolution of the solar wind plasma during its
propagation to 1 au and beyond. Observations by the Parker Solar Probe (PSP)
mission (Fox et al., 2016), that are for the first time carried out below 0.3
au, contribute to our understanding of these processes. In particular, the
SWEAP instrument (Solar Wind Electrons Alphas and Protons, see Kasper et al.
(2016)), provides a unique view of electrons in the near-Sun environment.
Electrons are the main contributors to the heat conduction that is one of the
processes regulating the flow of energy in the solar wind (e.g. Scime et al.,
1994).

During its first perihelion PSP was immersed in a slow but highly Alfvénic
solar wind stream emerging from a small equatorial coronal hole (Bale et al.,
2019). There, it detected electron VDFs composed of a cold dense core, and a
tenuous and hotter population, called “strahl” (German for “beam”), which
streams outward from the Sun along the interplanetary magnetic field (Halekas
et al., 2020; Berčič et al., 2020). In the solar wind another suprathermal
population of electrons is usually observed: the halo, which is strongly non-
Maxwellian, has higher temperatures and lower densities than the core, and it
is distributed at all pitch-angles (Verscharen et al., 2019b).

Many processes may contribute to shape this peculiar, three component VDF.
Since electrons are lighter than ions, solar gravity constrains them less
effectively. High energy electrons are then able to escape from the Sun into
interplanetary space. The electric field that arises extracts ions from the
solar atmosphere as well (Pierrard & Lemaire, 1996; Maksimovic et al., 1997).
In addition, electrons are strongly shaped by the magnetic forces acting in
the heliosphere. The divergence of the solar magnetic field and the
simultaneous conservation of the adiabatic invariants have a focusing effect
on the electrons and produce the strahl populations. On top of that, the
ubiquitous influence of waves and collisions on solar wind particles shapes
the electron VDF further.

Due to these marked asymmetries and non-thermal features, electrons carry the
heat flux in the solar wind (Feldman et al., 1975; Pilipp et al., 1987; Scime
et al., 1994). The electron heat flux bears information about the origin of
the solar wind and the mechanisms that generate it in the solar corona. By
observing the heat flux it is also possible to understand the topology of the
interplanetary magnetic field. During Encounter 1, PSP detected a
predominantly uni-polar anti-sunward heat flux (Halekas et al., 2021), as one
would expect to observe in the presence of a strahl going outward from the Sun
along open magnetic field lines. The heat flux also brings with it traces of
various processes underlying particle-particle and wave-particle iterations
that can significantly alter the structure of the electron VDFs. Indeed, by
analyzing the heat flux in the solar wind at different heliocentric distances
and in different solar wind regimes, it is possible to have an insight into
the mechanisms that regulate it by reducing the skewness of the electron VDFs
(Scime et al., 1994).

At 0.166 au, during PSP’s first perihelion, the heat flux, almost entirely
carried by the strahl, shows a clear anti-correlation with respect to the
electron plasma $\beta$, the ratio of the electron thermal pressure to the
magnetic pressure (Halekas et al., 2021). No significant correlation is
observed when one compares the heat flux to the collisonal age. This suggests
that non-collisional mechanisms are playing a predominant role in regulating
the electron heat flux during Encounter 1, rather than Coulomb collisions.

At small heliocentric distances of PSP’s first orbit, the halo fractional
density is substantially lower than that of the strahl (Halekas et al., 2020),
and lower than the already limited halo fractional density observed at 0.3 au
by the Helios spacecraft (Berčič et al., 2019). This is consistent with
observations made at different heliocentric distances, which show that while
the halo fractional density increases, the strahl fractional density decreases
when one moves further away from the Sun (Maksimovic et al., 2005; Stverák et
al., 2009; Gurgiolo et al., 2012). This also agrees with recent numerical
models which reproduce the halo formation from the scattering of the strahl
(Vocks et al., 2005; Roberg-Clark et al., 2019; Tang et al., 2020; Micera et
al., 2020b; Jeong et al., 2020).

The data taken by the FIELDS experiment (Bale et al., 2016) on-board PSP show
that the inner heliosphere is permeated with a wide range of whistler
fluctuations. Agapitov et al. (2020) and Cattell et al. (2021a) found an
appreciable presence of sunward whistler waves inside $0.3$ au, with
intermixed propagation angles with respect to the magnetic field. Both the
parallel and obliquely propagating whistler waves were characterized by large
amplitudes. In addition, according to Cattell et al. (2021b), in the time
intervals when whistler waves were observed, enhanced heat flux suppression
rates were found, suggesting that the scattering of electrons by these waves
may be a plausible mechanism to explain the regulation of the heat flux in the
near-Sun solar wind. Jagarlamudi et al. (2021) and Cattell et al. (2021b)
reported a significant broadening of the strahl pitch angle distribution in
coincidence with the detection of short-duration whistler wave trains by PSP,
confirming the fundamental role of the interaction between electrons and these
waves in regulating the near-Sun solar wind heat flux. It is clear that
understanding the origin of these waves and the consequences of their
interactions with electrons may shed some light on a key mechanism responsible
for the energy redistribution in the solar wind. At 1 au, in addition to the
oblique large-amplitudes whistler waves (Breneman et al., 2010; Cattell et
al., 2020), parallel small-amplitude whistler waves are also often observed
(Lacombe et al., 2014; Stansby et al., 2016; Tong et al., 2019).

There may be various potential sources of propagating whistler waves.
Whistlers can be driven by wave-wave interactions which lead to turbulent
cascades and to power law spectra of the magnetic field in the whistler regime
(Stawicki et al., 2001; Smith et al., 2006; Tang et al., 2020). However,
observations of whistler fluctuations are found to be regularly concentrated
in regions of the inner heliosphere where electrons are characterized by high
electron $\beta$ values, or at least higher than those occurring in the close
proximity (Jagarlamudi et al., 2020; Cattell et al., 2021b). A positive
correlation is also observed when one compares the values of electron $\beta$
with the magnitude of the pitch-angle scattering of the strahl (Pilipp et al.,
1987; Crooker et al., 2003; Berčič et al., 2019), suggesting that kinetic
instabilities self-generated by electrons are the possible source of the
observed whistler waves and possible candidates for shaping the non-thermal
features of the electron VDF. Whistler waves can be also generated by the so-
called whistler temperature anisotropy instability (Gary, 1993), which in turn
is produced by a temperature anisotropy of the electrons with
$T_{e,\perp}>T_{e,\parallel}$, where $\perp$ and $\parallel$ denote directions
perpendicular and parallel to the magnetic field. However, this
electromagnetic instability has the maximum growth rate parallel to the
magnetic field, which would not explain the generation of the observed highly
oblique whistler waves. Moreover, it turned out to have a negligible impact on
the scattering of the high-energy components of the electron VDF(Saito & Gary,
2007). Whistler heat flux instabilities are the most promising mechanism to
generate the whistler waves observed in the inner heliosphere. Indeed, Cattell
et al. (2020) found a clear dependence of the occurrence of the whistler waves
on the threshold of whistler heat flux instabilities. Halekas et al. (2021)
showed that during the first two PSP orbits the thresholds of whistler heat
flux instabilities (Vasko et al., 2019; Verscharen et al., 2019) were fully
bounding the observed electron heat flux values, and the strahl fractional
densities were constrained below the limit values predicted by these
instabilities.

In a weakly collisional magnetized plasma, heat flux instabilities can be
self-generated by counter-streaming populations of electrons. These conditions
occur frequently in the solar wind, where the strahl presents often an anti-
sunward drift with respect to the protons, and to maintain a balance in
current, the core also drifts, with a much lower speed, in the direction
opposite to that of the strahl (Feldman et al., 1975; Scime et al., 1994).
Core and strahl, with their relative drift velocities, can self-induce
microinstabilities which can reduce the heat flux below the limit values
provided by collisional models (Spitzer & Härm, 1953). The nature of the heat
flux instabilities can differ according to the plasma conditions. Depending on
the drift velocity of the strahl, electromagnetic parallel (Gary et al., 1975;
Gary & Feldman, 1977; Gary et al., 1994, 1999; Kuzichev et al., 2019; López et
al., 2019) or oblique (Komarov et al., 2018; Horaites et al., 2018; Vasko et
al., 2019; Verscharen et al., 2019; López et al., 2020; Micera et al., 2020b)
whistler heat flux instabilities (WHFIs) can be excited. The free energy of
the counter-streaming electron populations is converted into magnetic energy
in the form of parallel right-hand circularly polarized or oblique right-hand
elliptically polarized whistler-mode waves. A large variety of electrostatic
or hybrid instabilities (e.g. Gary, 1978; Shevchenko & Galinsky, 2010; Pavan
et al., 2013; Shaaban et al., 2018a; López et al., 2020) can also be self-
generated by the skewness of the electron VDF. However, these instabilities
were found to have significantly lower growth rates than WHFIs and do not
satisfy the resonant conditions to scatter the strahl into the halo
(Verscharen et al., 2019). We therefore focus here on WHFIs, especially on the
oblique WHFI, as a crucial source of whistler waves in the inner heliosphere.

The occurrence of microinstabilities in the solar wind is strictly related to
the evolution of its bulk parameters that vary as it propagates in the
heliosphere. If the solar wind expands spherically, then, in the absence of
wave-particle interactions, the temperatures should follow the double
adiabatic Chew–Goldberger–Low (CGL) model (Chew et al., 1956). According to
this model, due to the conservation of the two adiabatic invariants
$T_{\perp}/B$ and $T_{\parallel}B^{2}/n^{2}$, in a spherically expanding solar
wind, where both the radial magnetic field $B$ and the density $n$ decrease as
$R^{-2}$, with $R$ being the heliocentric distance, a marked temperature
anisotropy would result even starting from an initially isotropic distribution
($T_{\perp}\sim R^{-2}$ and $T_{\parallel}=const$).

However, due to kinetic instabilities, i.e. electron and proton firehose
instabilities, and to the resulting wave-particle interactions that produce a
simultaneous parallel cooling and perpendicular heating, these temperature
trends are not observed (Hellinger & Trávníček, 2008; Matteini et al., 2012;
Innocenti et al., 2019a; Micera et al., 2020a). In this respect, Innocenti et
al. (2020) highlighted the indirect role of solar wind expansion in heat flux
regulation: solar wind expansion can trigger or modify the evolution of
firehose instabilities, which can in turn contribute to heat flux regulation.

In the presence of multiple-component electron VDFs, firehose instabilities
(triggered as function of temperature anisotropy and plasma $\beta$ values)
are not the only instabilities that can alter the solar wind expansion trend:
heat flux instabilities, triggered by the presence of counter-streaming
electron populations, become a recurring possibility. Which specific heat flux
instability is triggered depends on the properties of counter-streaming
electrons.

For the whistler heat flux instabilities, the onset strongly depends on the
ratio between the strahl drift velocity and the Alfvén velocity,
$u_{s}/V_{A}$. When this ratio increases, the instability threshold becomes
lower and its growth rate is boosted (Verscharen et al., 2019; López et al.,
2020). In the inner heliosphere, where the magnetic field is predominately
radial ($B\sim R^{-2}$), the Alfvén velocity decreases as $R^{-1}$. If the
strahl drift velocity decreases with heliocentric distance slower than that,
we may expect that the threshold of the instability also decreases during the
solar wind passage through the inner heliosphere. This would mean that, as the
solar wind propagates, it would encounter more favorable conditions for the
quasi-continuous generation of unstable whistler-mode waves, as conjectured by
Verscharen et al. (2019).

In this work, we study the role of the solar wind expansion in triggering
electron kinetic instabilities and, hence, its role in regulating the heat
flux evolution in the solar wind. We investigate if the solar wind expansion
can drive the onset of whistler heat flux instabilities and hence contribute
to the generation of whistler waves in the inner heliosphere. We perform a
fully kinetic, Expanding Box Model Particle-in-Cell (PiC) simulation
(Innocenti et al., 2019b) and we infer how the impact of kinetic instabilities
varies during the solar wind propagation and how efficient they are in
limiting the heat flux at different heliocentric distances.

This paper is organised as follows. Section 2 describes the simulation
framework and setup, presents an overview of the code, the parameters used,
and initial conditions. In Section 3 we present and explain the main results.
In Section 4 we discuss the implications of our results for the understanding
of kinetic instabilities in the solar wind and for the interpretation of in
situ observations. Section 5 summarizes our results and reports the
conclusions.

##  2 Numerical Model

We study the onset and evolution of microinstabilities in the expanding solar
wind plasma self-consistently with EB-iPic3D simulations. The EB-iPic3D code
(Innocenti et al., 2019b) introduces for the first time the Expanding Box
Model (EBM) (Velli et al., 1992; Grappin & Velli, 1996; Tenerani & Velli,
2017) into a semi-implicit fully kinetic PiC framework (Brackbill & Forslund,
1982; Lapenta et al., 2006; Markidis et al., 2010; Innocenti et al., 2017).
This allows studies of the interplay between the solar wind plasma expansion
and kinetic processes down to the electron scales. In the EBM, a coordinate
transformation is used to incorporate expansion effects in the evolution
equations for fields and particles. This gives the possibility of reducing the
size of the computational domain and hence model the solar wind expansion at
an affordable computational cost (Innocenti et al., 2019b, a). In addition, in
semi-implicit PiC codes, the strict constraints on the temporal and spatial
scales needed to keep explicit PiC codes stable are almost completely removed
(Lapenta et al., 2017; Gonzalez-Herrero et al., 2018; Micera et al., 2020a).
These two innovative aspects together allow studies of the solar wind dynamics
from a kinetic point of view and to explore the long-term evolution and large-
scale effects of the expansion.

In the EBM, one follows a radially expanding solar wind plasma parcel as it
moves away from the Sun in the comoving frame. The average heliocentric
distance of the box center is therefore given by $R(t)=R_{0}+u_{0}~{}t$, where
$R_{0}$ is the initial distance from the Sun and $u_{0}$ is the constant
radial velocity at which the solar wind moves away from the Sun. The value
$\tau_{exp}=R_{0}/u_{0}$ defines the time scale of the expansion and thus
governs the rate at which the simulation parameters vary with the heliocentric
distance. All the expansion related quantities can be expressed via a
parametric dependence on $R(t)$.

We investigate the interaction between the adiabatic expansion and the
development of kinetic instabilities. Our objective is to follow how expansion
can lead to the generation of electromagnetic fluctuations that affect
particle VDFs. To do so, we initialise our simulation with realistic plasma
conditions, observed by PSP at its first perihelion, and we follow the non-
linear evolution of the electron VDF during its interaction with the waves
that are generated by the electrons themselves as the solar wind expands.

We carry out a two spatial dimensions, three velocity components (2D3V) fully
kinetic EBM simulation. The electron distribution function initially consists
of two populations. The strahl (subscript “s”) has a density $n_{s}$ that is
only $5\%$ of the total electron density $n_{e}=350$ cm-3 and a drift velocity
$u_{s}$ that is directed anti-sunward and, in agreement with observations,
equals $6900$ km/s (Halekas et al., 2020). The core (subscript “c”) contains
the remaining $95\%$ of the electron density and has a sunward drift $u_{c}$
relative to the proton rest frame, in order to satisfy zero-current condition
($n_{c}~{}u_{c}+n_{s}~{}u_{s}~{}=0$) and hence to ensure the quasi-neutrality
of the plasma. In agreement with recent PSP observations (i.e Halekas et al.,
2020; Berčič et al., 2020), we adopted an initially isotropic core with
$k_{B}T_{c}=30$ eV, where $k_{B}$ is the Boltzmann constant. The strahl is
characterised by a parallel temperature, $k_{B}T_{s\parallel}=173$ eV, twice
larger than the perpendicular temperature, in order to take into account its
limited angular extent. The ions, assumed to be only protons, have the real
proton-to-electron mass ratio $\mu=~{}m_{i}/m_{e}=~{}1836$, a zero drift
velocity in their reference frame, density $n_{i}=n_{e}$, and are initially
isotropic with $T_{i}=T_{c}$. The core and strahl initial VDFs are described
via drifting-Maxwellian and drifting-bi-Maxwellian distributions,
respectively:

$\displaystyle
f_{j}(v_{\parallel},v_{\perp},t=0)=\frac{(2\pi)^{-3/2}}{w_{j\perp}^{2}w_{j\parallel}}\exp\left(-\frac{v_{\perp}^{2}}{2w_{j\perp}^{2}}-\frac{(v_{\parallel}-u_{j})^{2}}{2w_{j\parallel}^{2}}\right),$
(1)

with $w_{j}=\sqrt{k_{B}T_{j}/m_{j}}$ and $u_{j}$ being respectively the
thermal and drift velocities of the species $j$ (core electrons, strahl
electrons and protons).

The initial background magnetic field is uniform and directed anti-sunward,
$\bm{B}_{0}=B_{0}\hat{e}_{x}$, with $B_{0}=60$ nT, so that the initial proton
Alfvén speed $v_{A}=B_{0}/\sqrt{4\pi n_{e}m_{i}}=0.00023\,c$, with $c$ being
the speed of the light in vacuum ($x$ is the parallel spatial coordinate, and
$y$ the perpendicular one). The boundary conditions adopted are periodic for
both particles and fields.

The simulation domain sizes are $L_{x}=L_{y}=8\,d_{i}$, resolved with
$800\times 800$ cells, where $d_{i}=c/\omega_{pi}$ is the ion skin depth, and
$\omega_{pi}=\sqrt{4\pi e^{2}n_{i}/m_{i}}=24615$ s-1 the ion plasma frequency
(at the beginning of the simulation), with $e$ being the elementary charge. We
use $1024$ particles per cell per species (the core and strahl electron
populations are treated as separated species in the code). The temporal step
is $\Delta t=0.05\,\omega_{pi}^{-1}$, and the simulation runs up to the time
$t=16500\,\omega_{pi}^{-1}$. With these parameters the initial relevant
frequency ratios are: $\omega_{pe}/\omega_{ce}=100$,
$\omega_{pi}/\omega_{ci}=4284.85$ and $\omega_{pi}/\omega_{ce}=2.33$, with
$\omega_{pe}=\sqrt{4\pi e^{2}n_{e}/m_{e}}$ being the electron plasma frequency
and $\omega_{cj}=eB_{0}/c\,m_{j}$ is the gyrofrequency of the species $j$. For
the spatial scales one should consider that: $d_{i}/\rho_{e}=56$, and
$d_{i}/d_{e}=42.85$ with $\rho_{e}=w_{e}/\omega_{ce}$ being the electron
gyroradius and $d_{e}=c/\omega_{pe}$ is the electron skin depth.

An important choice in our simulation is the expansion time $\tau_{exp}$, and
its ratio with respect to the e-folding time ($\tau$) of the instabilities we
aim to study. On the one hand, this ratio must be large enough to keep these
two timescales well separated. On the other hand, the effects of the expansion
must be fast enough to make the simulation feasible within reasonable
computation resources. The selection of $\tau_{exp}$ was done on the basis of
the typical growth rates of whistler heat flux instabilities for
characteristic solar wind parameters. To get an estimate of this value, we
relied on our previous work (Micera et al., 2020b), and on the results
obtained by López et al. (2020) who solved the linear dispersion relation for
a large range of plasma parameters. This was then followed by a convergence
study where simulations performed with different resolutions and expansion
times confirmed the results of the present study. We thus chose
$\tau_{exp}=10000\,\omega_{pi}^{-1}$, which gives
$\tau_{exp}/\tau_{~{}OWHFI}\approx 30$, where
$\tau_{~{}OWHFI}=1/\gamma_{~{}OWHFI}$ is the inverse of the maximum growth
rate $\gamma_{~{}OWHFI}$ of the oblique WHFI studied in Micera et al. (2020b).
$\tau_{exp}$ is lower than the realistic solar wind expansion time, but still
large enough to keep the timescales of the instability and of the expansion
well separated, hence our simulation results are physically significant. We
obtained this expansion time by maintaining the radial solar wind velocity
consistent with PSP observations during Encounter 1 ($u_{0}=540$ km/s =
$0.0018\,c$), but assuming a very small initial heliocentric distance of the
solar wind co-moving simulation box ($R_{0}=18~{}d_{i}$). It should be kept in
mind that $u_{0}$ and $R_{0}$ are free parameters and only their ratio
regulates how fast our plasma expands. This small $R_{0}$ is merely a
convenient value that gives a characteristic expansion time that fits the
criteria outlined.

The energy evolution in the simulation follows the expected adiabatic trend
with heliocentric distance. For example, when the initial magnetic field is
purely radial, the magnetic energy evolves as $E_{B}\sim B^{2}\sim R^{-4}$.
Magnetic energy fluctuations constitute only a small fraction of the total
magnetic energy (see Section 3.1).

##  3 Results

###  3.1 PiC simulation of wave-particle resonant interactions

We first verified, both through simulation and linear theory, that our plasma
system with realistic initial conditions described in Section 2 is stable with
respect to any kinetic instability if it does not expand. For the expanding
system, Figure 1(a) shows that the total magnetic fluctuation energy
normalized to the mean magnetic field energy $\delta B^{2}/B^{2}$ peaks at
$t\approx 3200\,\omega_{pi}^{-1}$ after a fast growth phase. Then it decreases
significantly until secondary instabilities are triggered. A second peak is
reached at $t\approx 13000\,\omega_{pi}^{-1}$, and finally the energy ratio
relaxes again until the end of the simulation. Due to the expansion, the
plasma parcel comoving with the solar wind becomes first unstable to the
oblique WHFI (Vasko et al., 2019; Verscharen et al., 2019; López et al., 2020;
Micera et al., 2020b), which is responsible for the first exponential rise of
the magnetic energy fluctuations. After these fastest growing modes are
saturated and the whistler instability completely relaxed, the energy of the
magnetic field fluctuations monotonically decreases. It then grows again when
new modes are destabilized due to the variation of plasma parameters caused by
the further plasma expansion.

The adiabatic expansion of the plasma produces a decrease of the core-strahl
differential velocity (Innocenti et al., 2020). This decrease is however less
steep than that of the local Alfvén speed, adiabatically evolving as $R^{-1}$
in the undisturbed solar wind with predominantly radial interplanetary
magnetic field. Consequently, the ratio $u_{s\parallel}/v_{A}$ as function of
the time, displayed in Figure 1(b), increases in the first stage of the
simulation, when the plasma expands following a double-adiabatic evolution,
until it reaches the oblique WHFI threshold. The instability onset produces a
further, sudden decrease of the strahl drift velocity and inverts the trend of
the ratio evolution. Thus, two concurrent mechanisms take place: on the one
hand, continuous expansion and the resulting excitation of whistler-mode
waves, and, on the other hand, significant regulation of the strahl drift
velocity by the generated waves and the subsequent saturation of the
instability. The oblique WHFI goes through cycles of stabilization and
destabilization, with effects clearly reflected in the evolution of
$u_{s\parallel}/v_{A}$.

In Figures 1(c), (d) we display, respectively, the tracks of the strahl and
core electron populations in the $\beta_{j\parallel}$ vs
$T_{j\perp}/T_{j\parallel}$ plane, with $\beta_{j\parallel}=8\pi
n_{j}k_{B}T_{j\parallel}/B^{2}$. Figure 1(c) confirms that the strahl starts
from a quasi-stable state, far from microinstability thresholds. After a first
phase of adiabatic expansion, the strahl enters in the cycle that leads to the
redistribution of the its kinetic energy from the parallel to the
perpendicular direction and vice versa, depending on whether the oblique WHFI
is in a growth or relaxation phase. The core also deviates from the double
adiabatic evolution when the first WHFI occurs (Figure 1(d)). However, its
largest contribution to the system evolution occurs in the last stages of the
simulation, when its track crosses the threshold of the oblique electron
firehose instability (EFI) (Li & Habbal, 2000; Gary & Nishimura, 2003;
Camporeale & Burgess, 2008; Shaaban et al., 2019a; Innocenti et al., 2019a)
and is successively bounced back to the stable region. In this way, the core
significantly participates in the generation of the late peak in the evolution
of the fluctuating magnetic energy of Figure 1(a). Thereby, the late increase
in energy can be traced back to an almost simultaneous stimulation of the
oblique EFI by the core and a secondary oblique WHFI by the strahl (this is
confirmed by the results of the linear dispersion analysis, see Section 3.2).
We note that in both Figures 1(c) and (d) there is no spread around the
average values of temperature anisotropy and $\beta_{\parallel}$ during the
adiabatic phase. The oscillations around the mean values increase when
instabilities develop, especially in the last stages of the simulation, when
oblique WHFI and EFI are triggered together.

Figure 1: Total magnetic energy fluctuations normalized to the the mean
magnetic field energy as a function of time (in the units of
$\omega_{pi}^{-1}$ and $\omega_{ce}^{-1}$) and heliocentric distance $R/R_{0}$
(panel (a)). Parallel strahl drift velocity normalized to the local Alfvén
velocity as a function of time and heliocentric distance (panel (b)). Strahl
(panel (c)) and core (panel (d)) trajectories in the $\beta_{\parallel}$ vs
$T_{\perp}/T_{\parallel}$ plane. Solid, dashed, and dotted black curves
represent isocontours of the oblique EFI for $\gamma=0.2,0.1$ and
$0.01~{}\omega_{ce}$, respectively (Gary & Nishimura, 2003). The dashed and
dotted blue curves are the isocontours of growth rates $\gamma=0.1$ and
$0.01~{}\omega_{ce}$ for the whistler temperature anisotropy instability,
respectively (Gary & Wang, 1996). Red and gray filled dots in panels (c) and
(d) correspond to temperature and $\beta$ values averaged over the whole
simulation box every $625\,\omega_{pi}^{-1}$, while the background shadings
depict the spread of the simulated data values. In all the four panels the
initial time is marked with a blue dot, the first onset of the oblique WHFI
($t=2500\,\omega_{pi}^{-1}$) with a yellow dot, the oblique WHFI + EFI second
growth stage ($t=11250\,\omega_{pi}^{-1}$) with a green dot, and the final
simulation time with a cyan dot.

In Figure 2 we show the total electron VDF as a function of $v_{\parallel}$
and $v_{\perp}$ at important stages of its evolution. The cuts of the total
electron VDF $f_{e}(v_{\parallel},v_{\perp}=0)$ and
$f_{e}(v_{\parallel}=0,v_{\perp})$ are shown in Figure 3. The excitation of
microinstabilities has noticeable effects on the electron distribution
function. The generated waves resonantly interact with the electron VDF,
modifying its bulk parameters as well as its shape. We depict in Figure 2(a)
the initial, stable electron VDF. Between $t=0$ and $t=2500\,\omega_{pi}^{-1}$
adiabatic expansions reduces the perpendicular velocity of the distribution.
For the strahl, this can be interpreted as adiabatic focusing in the expanding
solar wind.

At $t=2500\,\omega_{pi}^{-1}$, the VDF becomes unstable to WHFI, Figure 2(b).
At this point the triggered whistler waves start to resonate with the
electrons producing the “horns” that are visible at $v_{\parallel}>0$, which
correspond to a deviation of the suprathermal electron population from its
original bi-Maxwellian shape. The interaction of the electrons with the
whistler waves, is accentuated at $t=3000\,\omega_{pi}^{-1}$ (Figure 2(c)),
almost at the peak of the linear growth phase of the oblique WHFI, producing
multiple signatures of resonant interaction. A significant transfer of
electron momentum from the parallel to the perpendicular direction and a
considerable reduction of the strahl drift velocity can be observed, which
lead to a simultaneous broadening of the strahl pitch angle distribution.
Electrons with $v_{\parallel}$ satisfying the condition of resonance with the
growing whistler waves are scattered towards higher values of $v_{\perp}$,
forming a nascent halo that characterises the distribution function during the
whole propagation of the solar wind through the heliosphere.

In Figure 2(c) vertical lines mark the $v_{\parallel}$ values at which
electrons are expected to interact via cyclotron, Landau and anomalous
cyclotron resonances with the whistler waves which at that time permeate the
simulated plasma. The whistler resonance condition with electrons is given by
$v_{\parallel}=(n~{}\omega_{ce}+\omega_{r})/k_{\parallel}$, with $\omega_{r}$
and $k_{\parallel}$ being respectively the real wave frequency and the
parallel wave number of the fastest growing mode at a given time, with $n=-1$
for the cyclotron resonance, $n=0$ for the Landau resonance, and $n\geq 1$ for
the anomalous cyclotron resonances (Krall & Trivelpiece, 1973). The circles,
centered at $v_{\parallel}=\omega_{r}/k_{\parallel}$ and $v_{\perp}=0$,
illustrate the diffusion paths of the electrons that undergo anomalous
cyclotron resonant interactions, as constant energy surfaces in the reference
frame of the wave (Roberg-Clark et al., 2019; Verscharen et al., 2019). We
compute $\omega_{r}$ and $k_{\parallel}$, with their ratio giving the parallel
phase velocity of the waves ($v_{ph}=\omega_{r}/k_{\parallel}=0.00233~{}c$,
with $\omega_{r}=0.014~{}\omega_{pi}$ and $k_{\parallel}=6~{}\omega_{pi}/c$),
via the linear dispersion relation (see Figure 5(b) below) and fast Fourier
transform (FFT) in space of the out-of-plane component of the fluctuating
magnetic field $\delta B_{z}$ (see Figure 4(b) below), respectively.

One can see that the $n=1$ and $n=2$ resonant surfaces intersect quite
accurately the two large horn-like structures that appear in the frontal/anti-
sunward region of the electron VDF during the evolution of the oblique WHFI
(Figure 2(c)). This implies that a portion of the strahl electron population,
characterized by a $v_{\parallel}$ that fulfills the anomalous cyclotron
resonance conditions is scattered by whistler waves towards higher
perpendicular velocities. This leads to the isotropization of the suprathermal
electrons and to the formation of a nascent electron halo. Due to their
elliptical polarization, the oblique whistlers modes can trigger cyclotron
resonances with both positive and negative $n$ (Komarov et al., 2018). Thus,
electrons in the $v_{\parallel}<0$ region of the VDF are at this stage
efficiently scattered by the $n=-1$ cyclotron resonance (Figure 2(c)).

Shortly after the end of the linear growth phase of the oblique WHFI, at
$t=3500\,\omega_{pi}^{-1}$, the halo formation is still restricted mainly to
the portion of strahl with velocities satisfying the anomalous cyclotron
resonant conditions. The nascent halo appears as a marked deformation of the
distribution function, as depicted in Figure 2(d). However, this does not last
long as the electron VDF relaxes and the halo becomes more symmetric. In
Figure 2(e), at $t=5000\,\omega_{pi}^{-1}$, it can be seen that the overall
anisotropy of the electron VDF has been considerably reduced and a significant
amount of particles is present in the high $v_{\perp}$ \- negative
$v_{\parallel}$ region of the phase space. The combined effect of Landau
resonance, late stage of the anomalous cyclotron resonance, and normal
cyclotron resonance contributes to the isotropization of the total electron
VDF and the formation of a halo population distributed at all pitch angles.
These resonant wave-particle interaction mechanisms become predominant during
the relaxation phase of the oblique WHFI (Micera et al., 2020b), as the
whistler modes shift towards near-parallel propagation angles (see Figures
4(e) and (f)).

We calculated the wave frequency ($\omega_{r}=0.00622~{}\omega_{pi}$), used
for the evaluation of the resonant conditions at this time, by solving the
linear dispersion relation for the full spectrum of whistler heat flux
instabilities (see Section 3.2). We obtained the parallel wave number
($k_{\parallel}=9~{}c/\omega_{pi}$) from the FFT of $\delta B_{z}$ (see Figure
4(e)). We emphasize that as the solar wind propagates away from the Sun it
finds more favorable conditions for scattering processes to occur (Liewer et
al., 2001). Electrons resonate with the whistler waves at gradually higher $n$
and the values of $v_{\parallel}$ at which the anomalous cyclotron resonances
occur are becoming closer.

The competition between the perpendicular cooling of the plasma produced by
the expansion, and the perpendicular heating and pitch-angle scattering of the
strahl carried out by the WHFIs during the first part of the simulation,
results in the electron VDF displayed in Figure 2(f). At
$t=7000\,\omega_{pi}^{-1}$, the oscillating magnetic energy of the plasma is
close to local minimum (see Figure 1(a)) and the electron distribution
function appears stable with respect to kinetic instabilities. The VDF is
characterized by a field-aligned strahl, whose extent, however, is less
pronounced than that at the initial time of the simulation, and a suprathermal
halo that is present in all directions. The quasi-stable phase persists until
the moment when the second growth of whistler waves, combined with the
excitation of the EFI by the core, starts again to reshape the electron VDF.
In Figures 2(g),(h) and (i), where the electron VDF is displayed at $t=10500$,
$14000$ and $16000\,\omega_{pi}^{-1}$, respectively, it is possible to see the
reemergence of non-Maxwellian horn-like structures in the high $v_{\perp}$ \-
positive $v_{\parallel}$ region of the phase space. These features are
indicators of ongoing anomalous cyclotron resonant scattering processes.

Figure 2: Total electron distribution functions in the
$v_{\parallel}-v_{\perp}$ phase space at $t=0$ (panel (a)), $t=2500$ (oblique
WHFI first onset, panel (b)), $t=3000$ (oblique WHFI linear growth phase,
panel (c)), $t=3500$ (oblique WHFI non-linear relaxation phase, panel (d)),
$t=5000$ (oblique WHFI non-linear relaxation phase, panel (e)), $t=7000$
(WHFIs saturation, panel (f)), $t=10500$ (EFI onset, panel (g)), $t=14000$
(oblique WHFI + EFI non-linear relaxation phase, panel (h)), $t=16000$ (final
distribution, panel (i)). The time is in units of $\omega_{pi}^{-1}$. Black
and blue vertical lines in panels (c) and (e) indicate $v_{\parallel}$ values
at which gyroresonances of whistler waves with electrons are expected. They
are plotted at $v_{\parallel}=-0.041~{}c$ in panel (c) and
$v_{\parallel}=-0.02~{}c$ in panel (e) for $n=-1$, at
$v_{\parallel}=0.0457~{}c$ in panel (c) and $v_{\parallel}=0.0218~{}c$ in
panel (e) for $n=1$, at $v_{\parallel}=0.089~{}c$ in panel (c) and
$v_{\parallel}=0.0429~{}c$ in panel (e) for $n=2$, at
$v_{\parallel}=0.064~{}c$ in panel (e) for $n=3$. Red vertical lines, drawn at
$v_{\parallel}=0.00233~{}c$ in panel (c) and $v_{\parallel}=0.0007~{}c$ in
panel (e), respectively delineate the electron population that fulfills the
$n=0$ Landau resonant condition with whistler waves. The circles show
electrons diffusion paths due to the anomalous cyclotron resonance interaction
as constant energy surfaces in the frame moving with the parallel phase
velocity of the wave $v_{ph}=\omega_{r}/k_{\parallel}$. They are centered at
$v_{\parallel}=v_{ph}$, $v_{\perp}=0$, and their radius equals
$n~{}\omega_{ce}/k_{\parallel}$.

In Figure 3 the cuts of the total electron VDF along the parallel
($f_{e}(v_{\parallel},v_{\perp}=0)$) and perpendicular
($f_{e}(v_{\parallel}=0,v_{\perp})$) directions are shown at four different
stages of the simulation. The initial electron VDF in Figure 3(a) is
characterized by a pronounced field-aligned, anti-sunward directed “shoulder”,
which is signature of the presence of a strahl with limited angular extent (at
$t=0$, $T_{s\perp}/T_{s\parallel}=0.5$). In Figure 3(b), at
$t=3000\,\omega_{pi}^{-1}$, the effects of the wave-particle interactions are
already visible. The electron VDF starts to deviate from the Maxwellian shape.
A portion of the strahl electrons starts to resonate, and then the VDF
exhibits bumps corresponding to the “horns” that represent the seed for the
halo formation. At $t=5000\,\omega_{pi}^{-1}$ (Figure 3(c)), during the non-
linear stage of the oblique WHFI, the diffusion of the seed population towards
higher values of $v_{\perp}$ has already passed the peak of its activity. The
scattered electrons start to relax and are reorganized to form a suprathermal
halo distributed at all pitch angles at the expense of the strahl. The final
VDF cuts are depicted in Figure 3(d) for $t=16000\,\omega_{pi}^{-1}$. They
show the final state of the VDF after the electron non-thermal features are
reshaped by the destabilization and relaxation of the firehose and whistler
instabilities in the second part of the simulation as shown in Figure 2. In
comparison with Figure 3(a), one can see in Figure 3(d) that the number of
non-thermal electrons with high $v_{\parallel}$ significantly decreased.

Figure 3: Total electron VDF cuts along the parallel
($f_{e}(v_{\parallel},v_{\perp}=0$), orange solid line) and perpendicular
($f_{e}(v_{\parallel}=0,v_{\perp}$), purple solid line) directions at $t=0$
(panel (a)), $t=3000\,\omega_{pi}^{-1}$ (panel (b)),
$t=5000\,\omega_{pi}^{-1}$ (panel (c)) and $t=16000\,\omega_{pi}^{-1}$ (panel
(d)). Vertical lines in panels (b) and (c) are drawn at the intersections of
the resonant surfaces of $n=1,2,3$ with the $v_{\parallel}=0$ axis (see
Figures 2(c) and (e)).

To analyze the nature of the broad spectrum of the generated waves, in Figure
4 we display the FFT of the out-of-plane fluctuating magnetic field component
$\delta B_{z}$ at different times. In Figure 4(a), at $t=0$, the
$\mathit{k_{\parallel}~{}\mathchar 45\relax~{}k_{\perp}}$ plane shows only the
initial numerical noise. Due to the expansion and hence to the electron
kinetic instabilities that are triggered by it, the 2D power spectrum evolves
towards the situation shown in Figure 4(b), where it peaks at
$k_{\parallel}\approx 6\;\omega_{pi}/c$ and $k_{\perp}$ between $11$ and
$14\;\omega_{pi}/c$. This means that at $t=3000\,\omega_{pi}^{-1}$ the
triggered waves present highly oblique angles of propagation with respect to
the field, mainly from $61^{\circ}$ to $67^{\circ}$.

When the linear growth phase of the oblique WHFI has already come to the end,
and the subsequent relaxation phase of the fastest growing modes starts, the
power shifts towards higher $k_{\parallel}$ and smaller propagation angles.
Figures 4(c), (d) and (e) show a clear gradual transfer of the power towards
the background magnetic field direction. However, the wave magnetic energy at
these stages is still distributed at a range of propagation angles rather than
being exclusively concentrated along the parallel direction. This implies that
during the non-linear phase of the oblique WHFI, resonant scattering of
electrons by oblique waves can still occur, in agreement with Levinson &
Eichler (1992) and Komarov et al. (2018).

At $t=6000\,\omega_{pi}^{-1}$ (Figure 4(f)), the oblique WHFI is approaching
the saturated marginal state. The wave power has clearly faded, it is
concentrated in a region of the $k_{\parallel}>8\;\omega_{pi}/c$, and exhibits
also a purely field-aligned component that has become comparable to the
oblique one. It is during this phase of the simulation that most of the
redistribution of the scattered suprathermal electrons takes place to form a
symmetric halo.

After the quiet phase that follows the saturation of the WHFIs, the expansion
brings the electrons back to unstable conditions. At
$t=10500\,\omega_{pi}^{-1}$ the activity of electromagnetic waves starts again
and this time it is initially due to the core crossing the threshold of the
resonant EFI (see Figure 1(d)). The free energy associated to the electron
core temperature anisotropy is converted into the magnetic energy of wave
fluctuations. The maximum growth rate of the instability is again mainly in
the oblique direction, at an angle about $50^{\circ}$ with respect to
$\bm{B}_{0}$ (Figure 4(g)). The identification of these modes as EFI is
confirmed by the linear theory, which shows the non-propagating character of
this branch (see Section 3.2).

At $t=12500\,\omega_{pi}^{-1}$ in Figure 4(h) an additional contribution
appears in the $\delta B_{z}$ power spectrum. These modes are also oblique and
can be attributed to the excitation of a second surge of WHFIs (see Figures
1(b) and (c), and Figures 5(e) and (f) below). Their power is concentrated at
higher $k_{\parallel}$ and $k_{\perp}$ than the simultaneous modes produced by
the EFI. Finally, in Figure 4(i), where the FFT of $\delta B_{z}$ is computed
at $t=16000\,\omega_{pi}^{-1}$, we see the damping of both unstable mode waves
and their transition towards nearly parallel angles as expected for both
oblique EFI (Camporeale & Burgess, 2008; Micera et al., 2020a) and oblique
WHFI (Micera et al., 2020b).

Figure 4: Fast Fourier transform of the out-of-plane magnetic field
fluctuations ($\text{FFT}(\delta B_{z})$) at $t=0$ (initial power spectrum,
panel (a)), $t=3000$ (oblique WHFI linear growth phase, panel (b)), $t=3500$
(oblique WHFI non-linear relaxation phase, panel (c)), $t=4000$ (oblique WHFI
non-linear relaxation phase, panel (d)), $t=5000$ (oblique WHFI non-linear
relaxation phase, panel (e)), $t=6000$ (WHFI saturated marginal state, panel
(f)), $t=10500$ (EFI onset, panel (g)), $t=12500$ (oblique WHFI + EFI linear
growth phase, panel (h)), $t=16000$ (final power spectrum, panel (i)). The
time is in units of $\omega_{pi}^{-1}$.

###  3.2 Kinetic linear theory

Kinetic linear theory can help to understand the nature of the various wave
fluctuations and their impact on the electrons through resonant interactions.
We find numerically the unstable solutions by solving a Vlasov–Maxwell
dispersion formalism (Stix, 1992), using a root finder based on the Müller’s
method. This has been done by using the dispersion solver DIS-K in the
Maxwellian limit (López et al., 2021) and hence approximating the true
electron distribution with a drifting bi-Maxwellian model. The liner
dispersion relation confirms that initially our system is not subject to any
unstable mode. We then obtain the general dispersion and stability analysis
for the entire wave spectrum at three significant instants of the simulation.

In Figures 5(a) and (b) we show the instability growth rate and the real wave
frequency provided by the linear theory at $t=3000\,\omega_{pi}^{-1}$, for the
entire range of propagation angles. We confirm that at this time the system is
unstable to the right-hand elliptically polarized oblique WHFI (Vasko et al.,
2019; Verscharen et al., 2019; López et al., 2020): the fastest growing
whistler mode is located at $k_{\parallel}\approx 6\,\omega_{pi}/c$ and
$k_{\perp}\approx 16\,\omega_{pi}/c$, has a growth rate
$\gamma_{\text{max}}\approx 0.023\,\omega_{ce}=0.006\,\omega_{pi}$, and a real
wave frequency, $\omega_{r}\approx 0.055\,\omega_{ce}=0.014\,\omega_{pi}$. We
also notice the presence of unstable modes at much lower angles, around
$k_{\parallel}\approx 3\;\omega_{pi}/c$, $k_{\perp}\approx 0$. These modes
have a significantly lower growth rate than the oblique whistler fluctuations
and correspond to the simultaneous excitation of the left-hand circularly
polarized firehose heat flux instability (Shaaban et al., 2018a, b; López et
al., 2020). The latter modes exhibit negligible wave frequencies (non-
propagating modes with $\omega_{r}\approx 0$).

At $t=5000\,\omega_{pi}^{-1}$, in Figures 5(c) and (d), the strongest
instability for this regime remains the purely oblique whistler heat flux
instability, although in its relaxation phase, with maximum growth rate
$\gamma_{\text{max}}\approx 0.005\;\omega_{ce}=0.001\;\omega_{pi}$ at
$k_{\parallel}\approx 9\;\omega_{pi}/c$, $k_{\perp}\approx 11\;\omega_{pi}/c$,
and wave frequency for the fastest growing mode $\omega_{r}\approx
0.033\,\omega_{ce}=0.006\,\omega_{pi}$. Nevertheless, propagating whistler
modes aligned to the magnetic field (around $k_{\parallel}=9\;\omega_{pi}/c$)
are present at this stage, when the oblique WHFI is in its relaxation phase.
When these results are compared with the simulation, it should be noted that
at $t=5000\,\omega_{pi}^{-1}$ in Figure 4(e) the parallel modes are hardly
recognisable, but start to appear a bit later and only become comparable to
the oblique modes at $t=6000\,\omega_{pi}^{-1}$ (see Figure 4(f)).

When we solve the linear dispersion relation for the plasma parameters at
$t=11250\,\omega_{pi}^{-1}$ we observe a combination of oblique left-hand non-
propagating modes with $\omega_{r}\approx 0$, typical of the resonant electron
firehose instability (Li & Habbal, 2000; Camporeale & Burgess, 2008), and once
again propagating whistler modes with $\omega_{r}\neq 0$ (Figures 5(e) and
(f)). The EFI, which starts as soon as the core electrons approach the
instability threshold, as shown in Figure 1(d), excites modes characterized by
wave numbers that are slightly lower than those generated at the same time by
the oblique WHFI. As a result, the region of the wave spectrum characterized
by low $k_{\parallel}$ and $k_{\perp}$ also has the frequency values of
$\omega_{r}\approx 0$ (bottom-left corner of Figures 5(f)). The modes with
$k_{\parallel}>6\,\omega_{pi}/c$ present the frequency range typical of the
whistler waves that we have simulated, which is also the whistler frequency
range commonly observed in the near-Sun solar wind (below $0.1\;\omega_{ce}$)
(e.g. Agapitov et al., 2020; Jagarlamudi et al., 2021).

The results obtained from the kinetic linear theory are in good quantitative
agreement with those inferred from the simulation. The linear theory clearly
captures the oblique whistler modes responsible for the substantial
deformation of the electron VDF, identifies their transition to reduced
propagation angles, and finally, manages to distinguish the quasi-overlapping
whistler and firehose modes that are generated in the last phases of the
simulation. However, there are some minor differences between the linear
theory results and the FFT obtained from our simulation. They are essentially
due to the use of the Maxwellian distribution functions to describe the
electrons in the theory, while in the simulations they show clear departures
from this model.

Figure 5: Growth rates $\gamma$ (left column) and corresponding real wave
frequencies $\omega_{r}$ (right column) in the
$\mathit{k_{\parallel}~{}\mathchar 45\relax~{}k_{\perp}}$ plane obtained by
solving the linear dispersion relation. The results of the calculations are
shown at $t=3000\,\omega_{pi}^{-1}$ for the oblique WHFI linear growth phase
(panels (a) and (b)), at $t=5000\,\omega_{pi}^{-1}$ for the oblique WHFI non-
linear relaxation phase (panels (c) and (d)) and at
$t=11250\,\omega_{pi}^{-1}$ for the oblique WHFI + EFI linear growth stage
(panels (e) and (f)). Both growth rates and wave frequencies are in
$\omega_{ce}$, with
$\omega_{pi}(t=3000~{}\omega_{pi}^{-1})/\omega_{ce}(t=3000~{}\omega_{pi}^{-1})=3.86$,
$\omega_{pi}(t=5000~{}\omega_{pi}^{-1})/\omega_{ce}(t=5000~{}\omega_{pi}^{-1})=5.25$
and
$\omega_{pi}(t=11250~{}\omega_{pi}^{-1})/\omega_{ce}(t=11250~{}\omega_{pi}^{-1})=7.26$.

###  3.3 Solar wind heat flux

The evolution of the heat flux in the heliosphere can provide information
about how the solar wind is loaded with photospheric energy and how this
energy is redistributed as the wind propagates away from the Sun. The heat
flux is strongly regulated by kinetic phenomena and in particular by wave-
particle interaction processes (Bale et al., 2013; Halekas et al., 2020;
Cattell et al., 2021b). As the solar wind expands, it encounters different
plasma conditions and hence regimes characterized by different heat flux
regulation processes.

The solar wind heat flux is almost exclusively carried by electrons, and in
particular by the suprathermal electron populations (Feldman et al., 1975;
Pilipp et al., 1987). In the inner heliosphere, where the suprathermal portion
of the electron VDF is clearly dominated by the strahl (Halekas et al., 2020;
Berčič et al., 2020), it is this population of electrons that is the main
responsible for the heat flux in the solar wind. Given its marked anisotropy,
i.e. its being flattened along the direction of the magnetic field, the strahl
carries a heat flux with predominantly radial component. For these reasons,
the heat flux carried in the reference frame of the solar wind corresponds to
that carried by the strahl along the direction parallel to $\bm{B_{0}}$.

The heat flux can be defined as $Q_{s\parallel}=\frac{m_{s}}{2}\int
v_{\parallel}v^{2}f_{s}d^{3}v$, and split in three components (Feldman et al.,
1975; Scime et al., 1994): $Q_{s}=Q_{\text{enth},s}+Q_{\text{bulk},s}+q_{s}$,
with $Q_{\text{enth},s}=\frac{3}{2}~{}n_{s}m_{s}u_{s}w_{s}^{2}$ the strahl
electron enthalpy, $Q_{\text{bulk},s}=\frac{1}{2}~{}m_{s}n_{s}u_{s}^{3}$ the
energy flux associated to the bulk motion of the electrons, and
$q_{s}=\frac{m_{s}}{2}\int(v_{\parallel}-u_{s})(v-u_{s})^{2}f_{s}d^{3}v$ the
heat flux carried by the strahl in its frame of reference (skewness of the
VDF).

In agreement with Innocenti et al. (2020) and Micera et al. (2020b), one of
the main factors that leads to the reduction of the heat flux is the rapid
decrease of $u_{s\parallel}$, which is in turn due to the combined effect of
the interconnected expansion and development of electron-scale instabilities.
In Figure 6(a) we show the temporal evolution of the parallel strahl drift
velocity normalized to the speed of the light in vacuum, $u_{s\parallel}/c$
(cf. Figure 1(b), where the strahl drift velocity is normalized to the local
Alfvén velocity). It can be seen that, after the first phase of undisturbed
expansion, $u_{s\parallel}/c$ has a sudden decrease in correspondence with the
onset of the oblique WHFI. When this first oblique WHFI is saturated, the
parallel strahl drift velocity decreases less rapidly, but this only lasts
until the time when the core EFI and the second stage of the oblique WHFI are
triggered and further considerably decrease its value.

The evolution of the parallel and perpendicular components of the core and
strahl thermal velocities can be seen in Figure 6(b), which shows the trends
due to the successive phases of excitation and relaxation of the electron
kinetic instabilities. The pitch-angle scattering of the electrons during the
development of the WHFIs leads to a substantial redistribution of the strahl
thermal energy from the parallel to the perpendicular direction and hence to a
broadening of the strahl pitch-angle distribution. The process acts as a
driving force for the halo formation. Due to the generation of whistler waves,
the perpendicular cooling produced by the expansion has a limited impact not
only on the strahl, but also, to a lesser extent, on the core. Already at the
first onset of the oblique WHFI, the core shows a change of slope in the
evolution of both parallel and perpendicular components of $w_{c}$. A further
reduction of the core temperature anisotropy occurs as expected in
correspondence with the EFI growth phase.

Both the kinetic energy redistribution and the drop of the strahl drift
velocity are reflected in the evolution of the electron heat flux. In Figure
6(c) we show the field-aligned heat flux associated with the strahl in the
rest frame of the solar wind and the three terms in which it can be split, as
functions of time. All the heat flux components are normalized to the
saturation heat flux
$q_{max}=\frac{5}{2}~{}m_{e}(n_{c}w_{c}^{3}+n_{s}w_{s}^{3})$. Both enthalpy
and bulk components of the heat flux exhibit a sudden and steep decrease as
soon as whistler waves start to be generated. The total heat flux presents the
characteristic multi-slope behavior seen also in other electron bulk
parameters. Figure 6(c) also shows that the heat flux is mainly carried by the
suprathermal electron enthalpy, but the component that is affected by the
collisionless regulation more significantly is the one related to the bulk
motion, in agreement with Innocenti et al. (2020) and Micera et al. (2020b).
The electron heat flux in the rest frame of the strahl, $q_{s}$, presents
first a rapid growth simultaneous with the oblique WHFI onset, when the VDF is
more asymmetric, and a second growth, more moderate, at the EFI onset.

In Figure 6(d) we quantify the heat flux reduction, by showing the evolution
of its percentage decrease. At the end of our simulation, the solar wind has
dissipated about $80\%$ of the heat flux that it was carrying at the initial
time step. Almost a half of it was dissipated during the first interaction
stage between the electron VDF and the whistler waves. The heat flux
percentage decrease shows the typical multi-trend behaviour corresponding to
the effectiveness of the repeated wave-particle interactions.

Figure 6: Temporal evolution of the parallel strahl drift velocity normalized
to the speed of the light in vacuum (panel (a)). Parallel (solid lines) and
perpendicular (dashed lines) thermal velocities of the strahl (red) and core
(black) electrons as a function of time (panel (b)). Evolution of the heat
flux components carried by the strahl along the magnetic field direction:
$Q_{s}$ (black), $Q_{\text{enth},s}$ (red), $Q_{\text{bulk},s}$ (blue), and
$q_{s}$ (purple). All the heat flux components are normalized to
$q_{max}=\frac{5}{2}~{}m_{e}(n_{c}w_{c}^{3}+n_{s}w_{s}^{3})$ (panel (c)).
Percentage decrease of the total electron heat flux as a function of time
((panel (d)). In all the four panels, the upper axis shows the heliocentric
distance $R/R_{0}$, while the two bottom axes show the time in units of
$\omega_{pi}^{-1}$ and $\omega_{ce}^{-1}$. The initial time is marked with a
blue dot, the first oblique WHFI onset ($t=2500\,\omega_{pi}^{-1}$) with a
yellow dot, the oblique WHFI + EFI second growth stage
($t=11250\,\omega_{pi}^{-1}$) with a green dot, and the final simulation time
with a cyan dot, similarly to Figure 1.

##  4 Discussion

###  4.1 Comparison with observations

We now compare our results to recent observations. The implications of our
simulation can help to clarify the dynamics of electrons in the solar wind. We
have confirmed the effects of the whistler heat flux instabilities on the
electron VDF (Micera et al., 2020b), namely the scattering of the strahl
electrons into the halo (Stverák et al., 2008) and the significant heat flux
suppression simultaneous to the instability development (Halekas et al.,
2021).

We have demonstrated how the expansion can act as a driver for the excitation
of whistler instabilities and hence as a self-sustaining mechanism for the
generation of whistler waves, which seem to have many aspects in common with
those detected in the near-Sun solar wind by PSP. The fastest growing
whistler-mode waves triggered in our simulation present angles of propagation
with respect to the background magnetic field from $61^{\circ}$ to
$67^{\circ}$ (Figure 4(b)). This is in qualitative agreement with the highly
oblique whistler waves detected during the PSP Encounter 1 and reported by
Agapitov et al. (2020) and Cattell et al. (2021a).

The magnetic amplitude of such highly oblique whistler waves reaches values up
to $\lvert\delta B\rvert\sim 0.5~{}B$ (Figure 1(a)). These values are
consistent with large-amplitude whistler waves reported by Agapitov et al.
(2020) and Cattell et al. (2021a) in the near-Sun solar wind and with those
observed in the STEREO electric field data (Breneman et al., 2010; Cattell et
al., 2020), but are significantly higher than the amplitudes of parallel
whistler waves observed in Cluster and THEMIS search coil data, which are
typically less than $0.02~{}B$ (Lacombe et al., 2014; Tong et al., 2019). In
addition, in line with the results by Cattell et al. (2021a), the peak
amplitude of the simulated whistlers does not change much with the
heliocentric distance (see Figure 1(a)). This suggests that the whistlers’
effect on the solar wind electron VDF may occur over a wide region of the
heliosphere.

The values of the whistler wave frequency in our simulation fall in the
observed range as reported by Agapitov et al. (2020), Cattell et al. (2021a),
Jagarlamudi et al. (2021). Indeed, according to these works, the whistler
perturbations in the inner heliosphere have frequencies ranging from about
$0.03$ up to $0.2~{}\omega_{ce}$ (cf. Figure 5 above).

The perpendicular electron heating driven by the whistler waves limits the
development of the temperature anisotropy produced by the expansion (Innocenti
et al., 2019b). The temperature of the suprathermal electron populations, as
well as that of the core, diverge from the double adiabatic expansion law
(Chew et al., 1956) already when the first oblique WHFI is triggered and way
before reaching the EFI threshold (see Figures 1(c) and (d) and Figure 6(b)).
This suggests that whistler heat flux instabilities are a key mechanism to
explain the limited values of electron temperature anisotropy observed in the
solar wind, on a par with the firehose and whistler temperature anisotropy
instabilities that are most commonly considered in this respect. Similarly to
observations of the halo reported by Stverák et al. (2008), our simulation
shows that the suprathermal portion of the electron VDF (strahl and halo),
despite the expansion, never really approaches the threshold of the firehose
instability. It is rather concentrated in the stable region at the center of
the electron $\beta_{\parallel}$ vs $T_{\perp}/T_{\parallel}$ plane (so-called
“Brazil plot”, see Figure 1(c)). This plausibly means that, when looking
exclusively at the classical “Brazil” plot for the core and for the
suprathermal electrons, one does not have the full information on which
instability limits the temperature anisotropies. This is particular relevant
for the suprathermal component of the electron VDF (Stverák et al., 2008;
Berčič et al., 2019).

The whistler waves also produce perpendicular heating of the core, which
causes the core temperature anisotropy evolution to deviate from that expected
from a double adiabatic expansion. This corresponds to observations reported
by Jagarlamudi et al. (2020) and Cattell et al. (2021b), according to which
during the intervals of whistler wave detection there seems to be a
concomitant increase of the core electron $\beta$.

###  4.2 Origin of field-aligned sunward whistler waves

In this Section, we propose a mechanism for the self-generation of whistler
waves aligned with the background magnetic field, propagating mainly in the
sunward direction. They may be the result of the non-linear evolution of the
oblique WHFI and of the shift of the oblique whistler waves towards small
propagation angles.

In Figure 7 we show the space-time Fourier power spectrum, in the
$k_{\parallel}-\omega_{r}$ plane, obtained exclusively for field-aligned
whistler-mode waves. It was computed as the FFT of the cut along the parallel
direction of the out-of-plane component of the fluctuating magnetic field,
$\delta B_{z}(x,y=0)$. We have selected a time interval ranging from
$t=5000\,\omega_{pi}^{-1}$ to $t=8000\,\omega_{pi}^{-1}$ (relaxation phase of
the oblique WHFI). As shown in Figure 4, this is the stage of the simulation
in which the unstable modes are no longer concentrated only at high
propagation angles, but instead the field-aligned component has become
significant. We see that during the selected interval there is a considerable
portion of whistler waves propagating in the sunward direction and hence anti-
parallel to the electron heat flux. The fastest growing parallel and anti-
parallel whistler waves are characterized by wave numbers between $8$ and
$12\;\omega_{pi}/c$, confirming those shown in Figure 4(f), and have
frequencies below $0.1\;\omega_{ce}$, in line with those depicted in Figure 5
and with in situ observations (e.g. Agapitov et al., 2020; Cattell et al.,
2021a; Jagarlamudi et al., 2021).

This result has a fundamental importance. Indeed, for a long time it was
considered that most of the whistler waves in the solar wind, and consequently
most of the collisionless heat flux regulation, were due to the parallel WHFI,
which has its maximum growth rate in the direction of the magnetic field (Gary
et al., 1975, 1994). In support of this, Tong et al. (2019b) made a
statistical analysis of the occurrence of whistler waves around 1 au and
showed that there is a very good similarity between the properties of the
observed waves and those expected from the waves generated by the parallel
WHFI.

Nevertheless, it remains unexplained how the whistler waves generated by the
parallel WHFI instability can regulate the electron heat flux. Indeed, it has
been demonstrated through both linear theory (Shaaban et al., 2019b; Vasko et
al., 2020) and simulations (Kuzichev et al., 2019; López et al., 2019) that
these waves travel only parallel to the heat flux and consequently are unable
to scatter the suprathermal component of the electron VDF.

Our simulation suggests that the relaxation of the oblique WHFI is a plausible
source of whistler waves propagating at limited angles with respect to the
magnetic field and, more importantly, mainly in the sunward direction. These
anti-parallel waves are able to interact with the filed-aligned strahl and
thus to explain the significant non-collisional contribution to the heat flux
regulation (Salem et al., 2003; Bale et al., 2013; Halekas et al., 2021).

We suggest that the fast (compared to the expansion time scales) and
continuous cycle of excitation and relaxation of the oblique WHFI in the inner
heliosphere produces the whistler waves with propagation angles ranging from
parallel to oblique, and the portion of parallel whistlers propagating
sunward, in agreement with observations (Agapitov et al., 2020; Cattell et
al., 2021a). This mechanism of whistler wave generation by the expansion can
essentially occur in the inner heliosphere. As the distance from the Sun
increases, on the one hand, the relative velocity between the core and the
suprathermal component of the electron VDF decreases due to the combined
effect of kinetic instabilities and solar wind expansion (e.g. Innocenti et
al., 2020). On the other hand, $v_{A}$ stops to decrease steeply when the
transverse component of the magnetic field becomes dominant (Chew et al.,
1956): when the field is radial, then $v_{A}\sim R^{-2}/R^{-1}\sim R^{-1}$;
when the azimuthal component dominates, then $v_{A}\sim R^{-1}/R^{-1}=const$.
The combination of these two factors significantly inhibits the development of
the oblique WHFI (López et al., 2020), which requires values of $u_{s}/v_{A}$
above a certain threshold to be triggered.

At larger heliocentric distances, the question seems to be more intricate.
Regarding the parallel small-amplitude and commonly anti-sunward whistler
waves (Lacombe et al., 2014; Stansby et al., 2016; Tong et al., 2019), they
can be explained by the fact that with the increase of the heliocentric
distance, the suprathermal component of the electron VDF becomes progressively
more isotropic (e.g. Berčič et al., 2019) so that the plasma satisfies the
conditions for the occurrence of the parallel WHFI instability (see Kuzichev
et al., 2019; López et al., 2019; Micera et al., 2020b). As far as the highly
oblique and large-amplitude whistler waves observed by (Breneman et al., 2010)
and (Cattell et al., 2020) are concerned, the question of their origin remains
open. Around 1 au, where the expansion does not increase the $u_{s}/v_{A}$
ratio anymore, there may be other phenomena that can significantly modify the
structure of the electron VDF and bring the plasma into unstable conditions.
According to Cattell et al. (2021a), at 1 au oblique whistler waves are often
associated with stream interaction regions or even coronal mass ejections that
may give rise to conditions similar to those reproduced by Micera et al.
(2020b), where the instability is independent of the presence of the azimuthal
magnetic field component, and the large value of $u_{s}$ is the main driver of
the instability.

It is important to emphasise a crucial difference between the presented
results and those obtained by Micera et al. (2020b) in a non-expanding case
study. The process of relaxation of oblique modes towards parallel angles can
be seen in both expanding and non-expanding simulations. However, in the
simulation reported by Micera et al. (2020b), the electron VDF, and in
particular its suprathermal component, became almost isotropic after the
scattering on the oblique whistler waves and remained almost unperturbed until
a new instability was triggered: the resulting VDF was indeed suitable for the
excitation of the parallel WHFI. In addition to the relaxation of oblique
modes towards parallel angles, long after the relaxation of the oblique WHFI,
new modes appeared, which were exclusively anti-sunward directed and unrelated
to the oblique WHFI modes. In the present work, the expansion of the solar
wind has the effect of continuously bringing the electrons into anisotropic
conditions. The presence of anisotropic suprathermal electron populations
inhibits the growth of the parallel WHFI, and the high $k_{\parallel}$,
exclusively anti-sunward directed modes, which are typical of parallel WHFI,
are not seen in the expanding simulation. Therefore, the only modes with a
significant growth and propagating along the magnetic field direction are
those resulting from the relaxation of the oblique WHFI, which manifest also a
significant sunward component. Our results confirm the need to take into
account the solar wind expansion in order to have a more complete picture of
the evolution of the electron VDF during its interaction with whistler waves.

###  4.3 Limitations of our approach

We now comment on the approximations adopted in our numerical model. In order
to reduce the computational cost of our simulation, we significantly decreased
the expansion time scale with respect to realistic ones. If one considers a
heliocentric distance $R=30~{}R_{\odot}$ and the wind speed $u=600$ km/s, the
resulting realistic expansion time will be $\tau_{r}=R/u\sim 10^{4}~{}$s. If
this value is compared with our $\tau_{exp}=0.4$ s, then
$\tau_{exp}/\tau_{r}\sim 10^{-5}$. This means that the expansion in our
simulation has been sped up by a factor $10^{5}$ with respect to the realistic
expansion times. However, this is acceptable as the expansion time scale is
longer than the time scales of triggered instabilities. With the selected
expansion time we have $\tau_{exp}/\tau_{~{}WHFI}\approx 10$ and
$\tau_{exp}/\tau_{~{}OWHFI}\approx 60$, with $\tau_{~{}WHFI}$ and
$\tau_{~{}OWHFI}$ evaluated as the inverse of the maximum growth rates of the
simulated parallel and oblique WHFIs, respectively (see Section 3.2). This
confirms that the selected $\tau_{exp}$ value is within the limits of validity
of the expanding box model and gives us confidence about the physical
significance of our results.

The speedup of the expansion effects leads to the electron core firehose
instability triggered almost simultaneously with the second cycle of
destabilization of the oblique WHFI. The compression of the temporal scales
that separate the occurrence of these two instabilities intensifies their
effects on the electron VDF and hence on the heat flux that the electrons
carry. The heat flux carried by the electrons is almost entirely dissipated
while the wind travels a heliocentric distance of only $R=2.65~{}R_{0}$, which
is clearly not the case in reality.

We also note that the background magnetic field in our simulation is purely
radial. This causes the Alfvén velocity to decrease as $v_{A}\sim R^{-1}$
throughout the entire simulated time frame (as $B_{r}\sim R^{-2}$ and $n\sim
R^{-2}$), causing in turn the oblique WHFI to be destabilized repeatedly
during the solar wind expansion. However, this behavior is only expected in
the inner heliosphere. As the solar wind moves away from the Sun, the
transverse component of the magnetic field $B_{t}\sim R^{-1}$ becomes dominant
and causes the Alfvén velocity to become almost constant. This means that the
ratio $u_{s}/v_{A}$ does not increase anymore and, together with the depletion
of the strahl caused by the oblique WHFI itself, implies that this instability
is no longer excited (Verscharen et al., 2019; López et al., 2020), and the
simulated cycles of stabilization and destabilization come to an end.

Finally, we argue that in the expansion phase that precedes the development of
the microinstabilities, our plasma follows a double adiabatic evolution,
which, however, does not appropriately describe the true evolution of the
solar wind where collisionless instabilities, turbulence and collisions are
ubiquitous (Landi et al., 2012). An approach involving all these processes
together goes beyond the scope of this work and will be investigated in future
simulations.

Figure 7: Space-time Fourier power spectrum ($k_{\parallel}-\omega_{r}$) of
both sunward ($k_{\parallel}<0$) and anti-sunward ($k_{\parallel}>0$) whistler
waves propagating along the background magnetic field direction for the
interval $5000<t\,\omega_{pi}<8000$.

##  5 Conclusions

In this work, we used an electron VDF which faithfully reproduces those
recently observed during the PSP Encounter 1 (Halekas et al., 2020; Berčič et
al., 2020), and we followed its evolution with heliocentric distance, in order
to shed light on the kinetic processes that influence the electron dynamics
during solar wind propagation from the inner heliosphere onward. This was
carried out by performing a two spatial dimensions, three velocity components
(2D3V) fully kinetic Expanding Box Model simulation through which we could
monitor the nonlinear development of kinetic instabilities and study their
effects on the solar wind electrons. The results of the simulations have been
compared with those obtained through kinetic linear theory and discussed
within the context of recent PSP observational studies (e.g. Agapitov et al.,
2020; Cattell et al., 2021a, b; Jagarlamudi et al., 2021).

It has been demonstrated that the expansion of the solar wind can be
considered as one of the main drivers of the oblique whistler heat flux
instability. This instability leads to the generation of whistler waves with
main properties being in qualitative agreement with those observed in the
inner heliosphere by PSP. It has also been confirmed that these waves, which
begin to propagate predominantly in the oblique direction with respect to the
ambient magnetic field, can significantly modify the shape of the electron
VDF, leading to the formation of the halo at the expense of the strahl, and
consequently regulate the electron heat flux. In addition, due to the
continuing expansion, the plasma crosses the threshold of electron kinetic
instabilities multiple times, triggering a cycle of excitation and saturation
of whistler heat flux and firehose instabilities. This leads to the continuous
destabilization of waves which interact with the electron VDF during the
entire simulation.

PiC simulations have recently shown that the whistler waves produced by
whistler heat flux instabilities can strongly scatter the electron
suprathermal populations (Roberg-Clark et al., 2018; Komarov et al., 2018;
Micera et al., 2020b). One of the crucial differences of our study with
respect to these earlier works is that our initial electron VDF is completely
stable to microinstabilities and naturally evolves towards unstable conditions
due to solar wind expansion.

As already shown by Innocenti et al. (2020), solar wind expansion can affect
the heat flux regulation via a two-step process: by modifying the evolution of
collisionless instabilities that in turn affect the heat flux regulation. In
the present work, we verify this assumption using a more appropriate
description for the strahl and 2D simulation geometry, which allow for oblique
whistler heat flux and firehose instabilities to develop. We demonstrated that
the evolution of heat flux regulating instabilities is modified by the
expansion (cf. the results reported by Micera et al. (2020b)) and that this
has consequences for the heat flux evolution.

We have also shown that sunward-directed parallel whistler waves can result
from the relaxation of oblique whistler waves. This can help to solve the
apparent controversy that the field-aligned whistler waves generated by WHFI
propagate only parallel to the heat flux direction (Kuzichev et al., 2019;
López et al., 2019; Vasko et al., 2020) and fail to interact with the strahl,
despite non-collisional phenomena shaping the non-thermal features of the
electron VDF being ubiquitous in the solar wind (e.g. Salem et al., 2003;
Halekas et al., 2021).

Our work represents an important step in clarifying the role of small-scale
electron kinetic processes in the broader context of solar wind physics. It
adds new insight into viable mechanisms for heat flux regulation by
collisionless processes, by including the effects of plasma expansion in fully
kinetic simulations of relevant instabilities. Our results are consistent with
several key trends observed in the heliosphere and may provide an explanation
for the variation of the average whistler wave propagation angle with
heliocentric distance. The parameters of the simulated whistler waves, such as
frequencies and amplitudes, are in remarkable agreement with observational
data. The further approach of PSP to the Sun and its observational campaigns
coordinated with those carried out by Solar Orbiter during the radial
alignment of the two spacecraft, will be able to provide us with information
on even more pristine electron VDFs and their radial evolution. This will
bring fundamental insight regarding the complementary evolution of kinetic
instabilities and electron energy fluxes.

A.M. acknowledge L. Matteini for useful discussions and L. Franci and E.
Papini for helpful suggestions on theoretical and numerical aspects. This work
was supported by a PhD grant awarded by the Royal Observatory of Belgium to
one of the authors (A. M.). These simulations were performed on the
supercomputers SuperMUC (LRZ) and Marconi (CINECA) under PRACE allocations. A.
N. Z. thanks the European Space Agengy (ESA) and the Belgian Federal Science
Policy Office (BELSPO) for their support in the framework of the PRODEX
Programme. R.A.L acknowledges the support of ANID Chile through FONDECyT grant
No. 11201048. This research was supported in part by the NASA DRIVE HERMES
project, grant No. 80NSSC20K0604.

## References

  * Agapitov et al. (2020) Agapitov, O. V., Wit, T. D. d., Mozer, F. S., et al. 2020, ApJ, 891, L20, doi: 10.3847/2041-8213/ab799c
  * Bale et al. (2013) Bale, S. D., Pulupa, M., Salem, C., Chen, C. H. K., & Quataert, E. 2013, ApJ, 769, L22, doi: 10.1088/2041-8205/769/2/L22
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49, doi: 10.1007/s11214-016-0244-5
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Nature, 576, 237, doi: 10.1038/s41586-019-1818-7
  * Berčič et al. (2019) Berčič, L., Maksimović, , M., Landi, S., & Matteini, L. 2019, MNRAS, 486, 3404, doi: 10.1093/mnras/stz1007
  * Berčič et al. (2020) Berčič, L., Larson, D., Whittlesey, P., et al. 2020, ApJ, 892, 88, doi: 10.3847/1538-4357/ab7b7a
  * Brackbill & Forslund (1982) Brackbill, J. U., & Forslund, D. W. 1982, J Comput Phys, 46, 271, doi: 10.1016/0021-9991(82)90016-X
  * Breneman et al. (2010) Breneman, A., Cattell, C., Schreiner, S., et al. 2010, Journal of Geophysical Research (Space Physics), 115, A08104, doi: 10.1029/2009JA014920
  * Camporeale & Burgess (2008) Camporeale, E., & Burgess, D. 2008, J. Geophys. Res., 113, A07107, doi: 10.1029/2008JA013043
  * Cattell et al. (2021a) Cattell, C., Short, B., Breneman, A., et al. 2021a, A&A, 650, A8, doi: 10.1051/0004-6361/202039550
  * Cattell et al. (2021b) Cattell, C., Breneman, A., Dombeck, J., et al. 2021b, ApJ, 911, L29, doi: 10.3847/2041-8213/abefdd
  * Cattell et al. (2020) Cattell, C. A., Short, B., Breneman, A. W., & Grul, P. 2020, ApJ, 897, 126, doi: 10.3847/1538-4357/ab961f
  * Chew et al. (1956) Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, Proceedings of the Royal Society of London Series A, 236, 112, doi: 10.1098/rspa.1956.0116
  * Crooker et al. (2003) Crooker, N. U., Larson, D. E., Kahler, S. W., Lamassa, S. M., & Spence, H. E. 2003, Geophys. Res. Lett., 30, 1619, doi: 10.1029/2003GL017036
  * Feldman et al. (1975) Feldman, W. C., Asbridge, J. R., Bame, S. J., Montgomery, M. D., & Gary, S. P. 1975, J. Geophys. Res., 80, 4181, doi: 10.1029/JA080i031p04181
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7, doi: 10.1007/s11214-015-0211-6
  * Gary (1978) Gary, S. P. 1978, PhPl, 20, 47, doi: 10.1017/S0022377800021358
  * Gary (1993) —. 1993, Theory of Space Plasma Microinstabilities
  * Gary & Feldman (1977) Gary, S. P., & Feldman, W. C. 1977, J. Geophys. Res., 82, 1087, doi: 10.1029/JA082i007p01087
  * Gary et al. (1975) Gary, S. P., Feldman, W. C., Forslund, D. W., & Montgomery, M. D. 1975, J. Geophys. Res., 80, 4197, doi: 10.1029/JA080i031p04197
  * Gary & Nishimura (2003) Gary, S. P., & Nishimura, K. 2003, Physics of Plasmas, 10, 3571, doi: 10.1063/1.1590982
  * Gary et al. (1994) Gary, S. P., Scime, E. E., Phillips, J. L., & Feldman, W. C. 1994, J. Geophys. Res., 99, 23391, doi: 10.1029/94JA02067
  * Gary et al. (1999) Gary, S. P., Skoug, R. M., & Daughton, W. 1999, Physics of Plasmas, 6, 2607, doi: 10.1063/1.873532
  * Gary & Wang (1996) Gary, S. P., & Wang, J. 1996, J. Geophys. Res., 101, 10749, doi: 10.1029/96JA00323
  * Gonzalez-Herrero et al. (2018) Gonzalez-Herrero, D., Boella, E., & Lapenta, G. 2018, Comput Phys Commun, 229, 162, doi: 10.1016/j.cpc.2018.03.020
  * Grappin & Velli (1996) Grappin, R., & Velli, M. 1996, J. Geophys. Res., 101, 425, doi: 10.1029/95JA02147
  * Gurgiolo et al. (2012) Gurgiolo, C., Goldstein, M. L., Viñas, A. F., & Fazakerley, A. N. 2012, Annales Geophysicae, 30, 163, doi: 10.5194/angeo-30-163-2012
  * Halekas et al. (2020) Halekas, J. S., Whittlesey, P., Larson, D. E., et al. 2020, ApJS, 246, 22, doi: 10.3847/1538-4365/ab4cec
  * Halekas et al. (2021) Halekas, J. S., Whittlesey, P. L., Larson, D. E., et al. 2021, A&A, 650, A15, doi: 10.1051/0004-6361/202039256
  * Hellinger & Trávníček (2008) Hellinger, P., & Trávníček, P. M. 2008, Journal of Geophysical Research (Space Physics), 113, A10109, doi: 10.1029/2008JA013416
  * Horaites et al. (2018) Horaites, K., Astfalk, P., Boldyrev, S., & Jenko, F. 2018, MNRAS, 480, 1499, doi: 10.1093/mnras/sty1808
  * Innocenti et al. (2020) Innocenti, M. E., Boella, E., Tenerani, A., & Velli, M. 2020, ApJ, 898, L41, doi: 10.3847/2041-8213/aba591
  * Innocenti et al. (2017) Innocenti, M. E., Johnson, A., Markidis, S., et al. 2017, Advances in Engineering Software, 111, 3, doi: https://doi.org/10.1016/j.advengsoft.2016.06.011
  * Innocenti et al. (2019a) Innocenti, M. E., Tenerani, A., Boella, E., & Velli, M. 2019a, ApJ, 883, 146, doi: 10.3847/1538-4357/ab3e40
  * Innocenti et al. (2019b) Innocenti, M. E., Tenerani, A., & Velli, M. 2019b, ApJ, 870, 66, doi: 10.3847/1538-4357/aaf1be
  * Jagarlamudi et al. (2020) Jagarlamudi, V. K., Alexandrova, O., Berčič, L., et al. 2020, ApJ, 897, 118, doi: 10.3847/1538-4357/ab94a1
  * Jagarlamudi et al. (2021) Jagarlamudi, V. K., Dudok de Wit, T., Froment, C., et al. 2021, A&A, 650, A9, doi: 10.1051/0004-6361/202039808
  * Jeong et al. (2020) Jeong, S.-Y., Verscharen, D., Wicks, R. T., & Fazakerley, A. N. 2020, ApJ, 902, 128, doi: 10.3847/1538-4357/abb099
  * Kasper et al. (2016) Kasper, J. C., Abiad, R., Austin, G., et al. 2016, Space Sci. Rev., 204, 131, doi: 10.1007/s11214-015-0206-3
  * Ko et al. (1996) Ko, Y. K., Fisk, L. A., Gloeckler, G., & Geiss, J. 1996, Geophys. Res. Lett., 23, 2785, doi: 10.1029/96GL02449
  * Komarov et al. (2018) Komarov, S., Schekochihin, A. A., Churazov, E., & Spitkovsky, A. 2018, Journal of Plasma Physics, 84, 905840305, doi: 10.1017/S0022377818000399
  * Krall & Trivelpiece (1973) Krall, N. A., & Trivelpiece, A. W. 1973, Principles of plasma physics
  * Kuzichev et al. (2019) Kuzichev, I. V., Vasko, I. Y., Rualdo Soto-Chavez, A., et al. 2019, ApJ, 882, 81, doi: 10.3847/1538-4357/ab3290
  * Lacombe et al. (2014) Lacombe, C., Alexandrova, O., Matteini, L., et al. 2014, ApJ, 796, 5, doi: 10.1088/0004-637X/796/1/5
  * Landi et al. (2012) Landi, S., Matteini, L., & Pantellini, F. 2012, ApJ, 760, 143, doi: 10.1088/0004-637X/760/2/143
  * Lapenta et al. (2006) Lapenta, G., Brackbill, J. U., & Ricci, P. 2006, Physics of Plasmas, 13, 055904, doi: 10.1063/1.2173623
  * Lapenta et al. (2017) Lapenta, G., Gonzalez-Herrero, D., & Boella, E. 2017, J. of Plasma Phys., 83, 705830205, doi: 10.1017/S0022377817000137
  * Levinson & Eichler (1992) Levinson, A., & Eichler, D. 1992, ApJ, 387, 212, doi: 10.1086/171072
  * Li & Habbal (2000) Li, X., & Habbal, S. R. 2000, J. Geophys. Res., 105, 27377, doi: 10.1029/2000JA000063
  * Liewer et al. (2001) Liewer, P. C., Velli, M., & Goldstein, B. E. 2001, J. Geophys. Res., 106, 29261, doi: 10.1029/2001JA000086
  * López et al. (2020) López, R. A., Lazar, M., Shaaban, S. M., Poedts, S., & Moya, P. S. 2020, ApJ, 900, L25, doi: 10.3847/2041-8213/abaf56
  * López et al. (2019) López, R. A., Lazar, M., Shaaban, S. M., et al. 2019, ApJ, 873, L20, doi: 10.3847/2041-8213/ab0c95
  * López et al. (2021) López, R., Shaaban, S., & Lazar, M. 2021, Journal of Plasma Physics, 87, 905870310, doi: 10.1017/S0022377821000593
  * Maksimovic et al. (1997) Maksimovic, M., Pierrard, V., & Lemaire, J. F. 1997, A&A, 324, 725
  * Maksimovic et al. (2005) Maksimovic, M., Zouganelis, I., Chaufray, J. Y., et al. 2005, J. Geophys. Res., 110, A09104, doi: 10.1029/2005JA011119
  * Markidis et al. (2010) Markidis, S., Lapenta, G., & Rizwan-uddin. 2010, Math. Comput. Simul., 80, 1509, doi: 10.1016/j.matcom.2009.08.038
  * Marsch (2006) Marsch, E. 2006, Living Reviews in Solar Physics, 3, 1, doi: 10.12942/lrsp-2006-1
  * Matteini et al. (2012) Matteini, L., Hellinger, P., Landi, S., Trávníček, P. M., & Velli, M. 2012, Space Sci. Rev., 172, 373, doi: 10.1007/s11214-011-9774-z
  * Micera et al. (2020a) Micera, A., Boella, E., Zhukov, A. N., et al. 2020a, ApJ, 893, 130, doi: 10.3847/1538-4357/ab7faa
  * Micera et al. (2020b) Micera, A., Zhukov, A. N., López, R. A., et al. 2020b, ApJ, 903, L23, doi: 10.3847/2041-8213/abc0e8
  * Pavan et al. (2013) Pavan, J., Viñas, A. F., Yoon, P. H., Ziebell, L. F., & Gaelzer, R. 2013, ApJ, 769, L30, doi: 10.1088/2041-8205/769/2/L30
  * Pierrard & Lemaire (1996) Pierrard, V., & Lemaire, J. 1996, J. Geophys. Res., 101, 7923, doi: 10.1029/95JA03802
  * Pilipp et al. (1987) Pilipp, W. G., Miggenrieder, H., Montgomery, M. D., et al. 1987, J. Geophys. Res., 92, 1075, doi: 10.1029/JA092iA02p01075
  * Roberg-Clark et al. (2019) Roberg-Clark, G. T., Agapitov, O., Drake, J. F., & Swisdak, M. 2019, ApJ, 887, 190, doi: 10.3847/1538-4357/ab5114
  * Roberg-Clark et al. (2018) Roberg-Clark, G. T., Drake, J. F., Swisdak, M., & Reynolds, C. S. 2018, ApJ, 867, 154, doi: 10.3847/1538-4357/aae393
  * Saito & Gary (2007) Saito, S., & Gary, S. P. 2007, Geophys. Res. Lett., 34, L01102, doi: 10.1029/2006GL028173
  * Salem et al. (2003) Salem, C., Hubert, D., Lacombe, C., et al. 2003, ApJ, 585, 1147, doi: 10.1086/346185
  * Scime et al. (1994) Scime, E. E., Bame, S. J., Feldman, W. C., et al. 1994, J. Geophys. Res., 99, 23401, doi: 10.1029/94JA02068
  * Shaaban et al. (2019a) Shaaban, S. M., Lazar, M., López, R. A., Fichtner, H., & Poedts, S. 2019a, MNRAS, 483, 5642, doi: 10.1093/mnras/sty3377
  * Shaaban et al. (2018a) Shaaban, S. M., Lazar, M., & Poedts, S. 2018a, MNRAS, 480, 310, doi: 10.1093/mnras/sty1567
  * Shaaban et al. (2018b) Shaaban, S. M., Lazar, M., Yoon, P. H., & Poedts, S. 2018b, PhPl, 25, 082105, doi: 10.1063/1.5042481
  * Shaaban et al. (2019b) Shaaban, S. M., Lazar, M., Yoon, P. H., Poedts, S., & López, R. A. 2019b, Monthly Notices of the Royal Astronomical Society, 486, 4498, doi: 10.1093/mnras/stz830
  * Shevchenko & Galinsky (2010) Shevchenko, V. I., & Galinsky, V. L. 2010, Nonlinear Processes in Geophysics, 17, 593, doi: 10.5194/npg-17-593-2010
  * Smith et al. (2006) Smith, C. W., Hamilton, K., Vasquez, B. J., & Leamon, R. J. 2006, ApJ, 645, L85, doi: 10.1086/506151
  * Spitzer & Härm (1953) Spitzer, L., & Härm, R. 1953, Phys. Rev., 89, 977, doi: 10.1103/PhysRev.89.977
  * Stansby et al. (2016) Stansby, D., Horbury, T. S., Chen, C. H. K., & Matteini, L. 2016, ApJ, 829, L16, doi: 10.3847/2041-8205/829/1/L16
  * Stawicki et al. (2001) Stawicki, O., Gary, S. P., & Li, H. 2001, J. Geophys. Res., 106, 8273, doi: 10.1029/2000JA000446
  * Stix (1992) Stix, T. H. 1992, Waves in plasmas
  * Stverák et al. (2008) Stverák, Š., Trávníček, P., Maksimovic, M., et al. 2008, Journal of Geophysical Research (Space Physics), 113, A03103, doi: 10.1029/2007JA012733
  * Stverák et al. (2009) Stverák, v., Maksimovic, M., Trávníček, P. M., et al. 2009, J. Geophys. Res., 114, A05104, doi: 10.1029/2008JA013883
  * Tang et al. (2020) Tang, B., Zank, G. P., & Kolobov, V. I. 2020, ApJ, 892, 95, doi: 10.3847/1538-4357/ab7a93
  * Tenerani & Velli (2017) Tenerani, A., & Velli, M. 2017, ApJ, 843, 26, doi: 10.3847/1538-4357/aa71b9
  * Tong et al. (2019b) Tong, Y., Vasko, I. Y., Artemyev, A. V., Bale, S. D., & Mozer, F. S. 2019b, ApJ, 878, 41, doi: 10.3847/1538-4357/ab1f05
  * Tong et al. (2019) Tong, Y., Vasko, I. Y., Pulupa, M., et al. 2019, ApJ, 870, L6, doi: 10.3847/2041-8213/aaf734
  * Vasko et al. (2019) Vasko, I. Y., Krasnoselskikh, V., Tong, Y., et al. 2019, ApJ, 871, L29, doi: 10.3847/2041-8213/ab01bd
  * Vasko et al. (2020) Vasko, I. Y., Kuzichev, I. V., Artemyev, A. V., et al. 2020, Physics of Plasmas, 27, 082902, doi: 10.1063/5.0003401
  * Velli et al. (1992) Velli, M., Grappin, R., & Mangeney, A. 1992, in American Institute of Physics Conference Series, Vol. 267, Electromechanical Coupling of the Solar Atmosphere, ed. D. S. Spicer & P. MacNeice, 154–159, doi: 10.1063/1.42861
  * Verscharen et al. (2019) Verscharen, D., Chandran, B. D. G., Jeong, S.-Y., et al. 2019, ApJ, 886, 136, doi: 10.3847/1538-4357/ab4c30
  * Verscharen et al. (2019b) Verscharen, D., Klein, K. G., & Maruca, B. A. 2019b, Living Reviews in Solar Physics, 16, 5, doi: 10.1007/s41116-019-0021-0
  * Vocks et al. (2005) Vocks, C., Salem, C., Lin, R. P., & Mann, G. 2005, ApJ, 627, 540, doi: 10.1086/430119