Collisionless Heat Flux Regulation via the Electron Firehose Instability in the Presence of a Core and Suprathermal Population in the Expanding Solar Wind

Electrons play a fundamental role in solar wind dynamics, because they drive coronal expansion and carry the greatest share of the heat flux due to their low mass (Feldman et al. 1975). A comprehensive solar wind theory then requires an understanding of the mechanisms behind electron heat flux regulation in the expanding solar wind. Outstanding open questions are whether collisional, Spitzer–Härm (Spitzer & Härm 1953; Salem et al. 2003), or collisionless  (Scime et al. 1994; Crooker et al. 2003; Landi et al. 2012) processes (the transition between regimes is discussed, e.g., in Bale et al. 2013 and references therein) act as the main heat flux regulator and, in the latter case, which specific collisionless process (i.e., instability) plays the most prominent role. Candidate instabilities are proposed (or ruled out) in Gary et al. (1975, 1994), Scime et al. (1994, 2001), Roberg-Clark et al. (2018), Komarov et al. (2018), Tong et al. (2019), Vasko et al. (2019), Verscharen et al. (2019), López et al. (2019), and Kuzichev et al. (2019).

The closer one gets to the Sun, the more the electron energy flux problem intersects the coronal heating problem, and the fundamental issue of solar wind acceleration (Scudder & Olbert 1983; Dorelli & Scudder 1999; Meyer-Vernet 1999; Landi & Pantellini 2001; Lie-Svendsen et al. 2001; Dorelli & Scudder 2003).

The role of collisions in heat flux regulation falls within the larger topic of their role in solar wind evolution  (Lie-Svendsen et al. 1997, 2001; Landi et al. 2012; Yoon et al. 2019). In approaching the discussion, one should keep in mind that collisions and wave–particle interactions ultimately act in the same direction, i.e., reducing the anisotropy of velocity distributions driven by the solar wind expansion, and that higher than observed collisionality levels may in fact approximate the effects of wave–particle interactions in models that, "on paper," should not account for them.

Solar wind expansion plays a fundamental role in heat flux evolution. Even in the simplified (and unrealistic) double-adiabatic (DA) framework (Chew et al. 1956), where wave–particle interactions and dissipative processes are neglected, one can expect a drop of the electron energy flux with heliocentric distance  (Scime et al. 1994). Abandoning the DA framework in favor of a more realistic solar wind description, where wave–particle interaction occurs, the possibility arises that expansion modifies solar wind bulk parameters in a way that facilitates the onset of wave–particle resonances and kinetic instabilities that can, in turn, regulate the heat flux.

Solar wind expansion increases the parallel beta and reduces the perpendicular-to-parallel temperature ratio: this drives the system toward ion (Hellinger et al. 2003; Matteini et al. 2006) and electron (Innocenti et al. 2019a) firehose instabilities, which in fact are observed to constrain the ion and electron populations in the respective "Brazil" plots (Matteini et al. 2007, 2013; Štverák et al. 2008; Berčič et al. 2019).

The electron firehose instability (EFI) is an electromagnetic kinetic instability that develops in the presence of a background magnetic field and of a T_⊥ < T_∥ thermal anisotropy over spatial and temporal scales that are relatively large and slow for electrons, and exhibits lower threshold and higher growth rates at oblique, rather than parallel, propagation (Paesold & Benz 1999; Li & Habbal 2000; Gary & Nishimura 2003; Camporeale & Burgess 2008).

Studies of the EFI are usually limited (with few exceptions, e.g., Shaaban et al. 2018, 2019) to a single, non-drifting, electron population, which therefore does not carry significant heat flux in the mean velocity frame. To understand heat flux dynamics, instead one has to take into account the multiple electron populations which compose the electron velocity distribution function (eVDF) in the solar wind.

In this paper, we discuss global electron heat flux regulation by the EFI triggered as a result of solar wind expansion within a purely collisionless description of the solar wind. We simulate plasma expansion self-consistently with the semi-implicit, fully kinetic expanding box model (EBM) code EB-iPic3D (Innocenti et al. 2019a, 2019b). The EBM framework maps the evolution of a solar wind plasma parcel that moves radially away from the Sun with constant velocity U₀, while expanding in the transverse directions with characteristic expansion timescale τ = R/U₀ (R is the heliocentric distance), to a Cartesian, non-expanding, co-moving grid where the secular evolution with distance appears through time-dependent terms and coefficients. In EB-iPic3D, this grid is resolved with a fully kinetic, semi-implicit algorithm (Markidis et al. 2010; Innocenti et al. 2017). Other EBM implementations rely on hybrid (Liewer et al. 2001; Matteini et al. 2006; Hellinger & Trávníček 2013) or magnetohydrodynamic (Tenerani & Velli 2017) descriptions.

The solar wind eVDF is composed of three populations: a colder, nearly Maxwellian core, a tenuous, suprathermal halo, and a field-aligned strahl  (Feldman et al. 1975; Pilipp et al. 1987; Maksimovic et al. 2005; Horaites et al. 2018). Since the aim of this work is to study heat flux evolution, this eVDF can be simplified with a two-population eVDF composed of a core component and a tenuous, suprathermal population (representing both the halo and the strahl) which we call "suprathermal electrons." The heat flux resulting from such an eVDF presents characteristics similar enough to the observed heat flux (see Section 3) to justify the approximation.

In Section 2, we investigate how plasma expansion drives our two-electron-component ("2E") simulation into the EFI-unstable area. For comparison, we evolve also a one-electron-component ("1E") simulation. In Section 3 we then focus on the role of the EFI in regulating the heat flux in the presence of multiple electron populations. The formula of Feldman et al. (1975) is used to distinguish among enthalpy, bulk, and heat flux in the frame of reference of the species components within the total electron energy flux.

We compare the development of the oblique EFI in the 2E and 1E simulations. Our simulation setup is quite similar to that of López et al. (2019). We simulate a 1D box with length L_x/d_i = 16, with d_i the ion skin depth, resolved with 1024 cells. The mass ratio is m_r = 1836 and the Alfvén speed is v_A/c = 0.00023, with c the speed of light: both electron skin depth and the electron gyroradius are resolved. The time step is ${\omega }_{{pi}}t=0.0375$, with ω_pi the ion plasma frequency. The ratio between the plasma frequency and the gyrofrequency is ω_pi/Ω_ci = 4390 and ω_pe/Ω_ce = 102 respectively for ions and electrons, compatible with values observed in the solar wind. The magnetic field at initialization is inclined with an angle ${\theta }_{0}=26^\circ$ with respect to the radial, x direction. In the 2E simulation, the core ("c") and suprathermal ("s") populations are modeled with Maxwellian distributions with thermal velocity v_th,c/c = 0.0119512, v_th,s/c = 0.0292744 and density n_c/n_e = 0.95, n_s/n_e = 0.05 respectively, with n_e the total electron density. Their drift velocity in the radial, x direction is V_x,c/c = −0.000479551 and V_x,s/c = 0.00911147, which satisfies the zero-current relation (Feldman et al. 1975; Scime et al. 1994). The ions have the same temperature as the core electrons, v_th,i/c = 0.00027. In the 1E simulation, the entire electron population has v_th,e/c = 0.0119512 and is not drifting with respect to the ions. A total of 5000 particles per species per cell are used in all cases.

The EB expansion time, τ, is ω_piτ = R₀/U₀ ∼ 5500, with R₀/d_i = 20 the initial distance of the box from the Sun, and ${U}_{0}/c=0.00358536$ the solar wind velocity. In EB simulations R₀ should be regarded as a sort of free parameter, chosen in order to have expansion dynamics fast enough to impact the evolution of the simulation, but slow with respect to the characteristic timescales of the processes under investigation (in this case, EFI evolution), as assumed in the formulation of the EB method. In our choice of simulation parameters we respect the solar wind ordering ${v}_{A}\sim {v}_{\mathrm{th},i}\lt \lt {U}_{0}\lt \,{v}_{\mathrm{th},e}$.

In Figure 1, panel (a), we show the oscillating magnetic field energy for the 2E (blue) and 1E (red) simulations. The global evolution of the magnetic and kinetic energy in the two simulations matches remarkably well the DA expectations. Observing the evolution of the oscillating magnetic field energy (which constitutes a small fraction of the total magnetic energy), we notice, in both simulations, the onset of an oblique EFI with growth rate γ/ω_pi = 0.0025 and γ/ω_pi = 0.0037 respectively, for the case with two and one electron species. We notice that these growth rates give a ratio between the expansion timescale τ and the EFI e-folding time t_e–f of τ/t_e–f ∼ 14 and τ/t_e–f ∼ 20, well in line with the τ > t_e–f requirement of EBMs.

The identification of the developing instability as an oblique EFI is supported by Figure 1, panel (b), where we plot the traces of the two simulations in the electron β_∥ versus T_⊥/T_∥ plane. The blue and red traces are the core and suprathermal electrons in the 2E simulation, and the cyan trace corresponds to the electrons in the 1E simulation. The traces "bounce back" in the stable area of the plot after briefly entering the EFI-unstable area. This is a departure from the "ideal expansion-driven traces," depicted with dotted lines, which show how the simulation would evolve due to expansion alone, i.e., in the absence of instability development. The differences between the 2E and 1E simulation traces are minimal, since the density of the suprathermal electrons is very low when compared to the core density in the 2E simulation. Figure 1, panel (c), shows the thermal velocity evolution of the core and suprathermal electrons, parallel and perpendicular components, as a function of R/R₀ compared with the expansion-driven expected behavior. We see that the observed evolution (solid lines) matches well the expansion-driven one (circles) before the onset of the instability. After that, energy is redistributed from the parallel to the perpendicular direction.

We depict in Figure 1, panel (d), the eVDF at R/R₀ = 1.88 (third dot in panels (a) and (b)), immediately after the end of the linear growth phase of the EFI. In the online animated version of Figure 1 the eVFD evolution during the entire simulation for the 2E and 1E simulations is provided.

In the video, we see a first phase of DA-like cooling. Since the initial magnetic field is not purely radial, the cooling affects both the parallel and the perpendicular direction. At R/R₀ = 1.88, consistently with the oscillating magnetic field energy evolution shown in Figure 1, panel (a), the 1E eVDF starts showing traces of resonant wave–particle interaction, similar to those observed in Innocenti et al. (2019a) in 2D3V EBM simulations. The core electrons of the 2E simulation exhibit similar traces. When the entire electron population of the 2E simulation is plotted, two features are quite evident: the reduction of the initial drift between the core and the suprathermal component, and the emergence of an asymmetry between the v_∥ < 0 and v_∥ > 0 sides of the eVDF.

The former feature is associated with a reduction of the observed heat flux, as discussed in the next section.

In the 2E simulation, the asymmetry of the initial eVDF in the v_x < 0 and v_x > 0 semi-planes originates an electron heat flux, which evolves with the eVDF.

In our simulation, the rest frame of the entire electron population, core plus suprathermal electrons, corresponds, due to the zero-current relation, to the co-moving frame, i.e., the rest frame of the solar wind. The total electron heat flux (which is also the total electron energy flux in this frame; Feldman et al. 1975) is composed of a core and suprathermal contribution. Labeling the electron populations, core and suprathermal, with j, the heat flux associated with population j, Q_j, is composed of three terms (Feldman et al. 1975):

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{Q}}}_{j} & = & {m}_{e}/2\displaystyle \int {\boldsymbol{v}}{v}^{2}{f}_{j}{d}^{3}v={n}_{j}{{\boldsymbol{V}}}_{d,j}(3/2{T}_{\parallel ,j}+{T}_{\perp ,j})\\ & & +\,{m}_{e}/2\,{n}_{j}{{\boldsymbol{V}}}_{d,j}{V}_{d,j}^{2}\\ & & +\,{{\boldsymbol{q}}}_{j}={{\boldsymbol{Q}}}_{\mathrm{enth},j}+{{\boldsymbol{Q}}}_{\mathrm{bulk},j}+{{\boldsymbol{q}}}_{j},\end{array}\end{eqnarray} \tag{ 1 }$$

where v is the velocity and V_d,j the species drift velocity in the co-moving frame. The heat flux in the frame of reference of each population, q_j, is

$$\begin{eqnarray}&&{{\boldsymbol{q}}}_{j}={m}_{e}/2\int {{\boldsymbol{w}}}_{j}{w}_{j}^{2}{f}_{j}{d}^{3}v\end{eqnarray} \tag{ 2 }$$

with ${{\boldsymbol{w}}}_{j}={\boldsymbol{v}}-{{\boldsymbol{V}}}_{d,j}$. Q_enth,j and Q_bulk,j are the enthalpy and bulk velocity contributions. While Q_enth,j and Q_bulk,j depend on "large-scale" solar wind parameters, such as parallel and perpendicular temperatures and drift velocities, q_j is a direct measure of the skewness of the VDF. In Scime et al. (1994), it is noted that Equation (1) excellently reproduces Ulysses energy flux observations between 1 and 5 au.

We observe that the suprathermal electrons carry the largest contribution to the energy flux, and that the two population contributions are oppositely directed, with the suprathermal contribution being anti-Sunward (Feldman et al. 1975; Scime et al. 1994). Since the velocity drifts are initialized in the radial, rather than parallel, direction, smaller perpendicular energy flux components (not depicted here, and generally negligible in the solar wind) are observed together with the dominant parallel ones. In Figure 2, we plot in red the suprathermal enthalpy, bulk and heat flux q_s (panels (a), (b), and (d)) contributions to the parallel textcolorredsuprathermal energy flux Q_s, customarily normalized to the free streaming value ${q}_{\max }=3/2{m}_{e}{\sum }_{j}{n}_{j}{v}_{\mathrm{th},j,0}^{3}$ (where we have considered that the two electron species have different thermal velocities).

Zoom In Zoom Out Reset image size

In the simulations as in the solar wind, the convection of the suprathermal electron enthalpy, ${Q}_{\parallel ,\mathrm{enth},s}$, largely dominates the energy flux (Feldman et al. 1975; Scime et al. 1994) (compare the y-axis scales in panels (a), (b), (d)).

The electron energy flux decreases during the simulation due to a combination of three processes: plasma expansion, reduction of the relative drift velocity of the populations with respect to the solar wind rest frame, and onset and development of the EFI.

To disentangle the role of these three processes, we plot in Figure 2, panels (a) and (b), two fits for the parallel enthalpy and bulk components. The first fit, in cyan, is calculated from Equation (1) using as thermal velocities the DA evolution (depicted as circles in Figure 1, panel (c)). In this fit, the module of the drift velocity is kept constant, and the parallel component is calculated projecting the initial drift velocity value, in the x direction, in the direction of the background magnetic field at each heliocentric distance, as calculated from DA evolution (since the radial magnetic field component drops with R faster than the transverse component, the angle of the magnetic field vector with the radial direction increases with R). For the second fits, in black, we use expansion-driven thermal velocity evolution (as before) and, for the drift velocity term, the values calculated at each heliocentric distance from the simulation.

We observe in both cases that the "black" fit excellently reproduces the observed values. The "cyan" fit is poorer. In the particular case of the enthalpy plot, panel (a), we notice that the cyan fit is adequate only up to EFI onset. This means that solar wind expansion alone (cyan fit), while contributing to the observed enthalpy and bulk component evolution, does not completely explain it. Collisionless processes, reproduced within our fully kinetic, EB simulation, reduce the heat flux to values lower than expected from expansion effects alone, as observed in the solar wind already in Scime et al. (1994).

In Figure 2, panel (c), we depict the parallel component of the drift velocity for the core (blue, left axis) and the suprathermal (red, right axis) electrons, as a function of R/R₀ and mark with vertical lines the linear growth phase of the EFI. The dotted lines, blue and red, depict the parallel drift component as expected from DA evolution (the dotted lines were used to calculate the cyan fits in panels (a) and (b)) for the two electron components. We notice that the sudden drops in the drift velocities are associated to the EFI. The EFI thus affects enthalpy and bulk components chiefly through its effect on the species drift velocities, rather than by redistributing energy from the parallel to the perpendicular direction (the expansion-driven thermal velocity evolution, which does not account for EFI energy redistribution, still gives a good fit of the observed values; black fits in Figure 2, panels (a) and (b)). Finally, in Figure 2, panel (d), we plot the parallel suprathermal electron heat flux in the rest frame of the suprathermal electrons, ${q}_{\parallel ,s}$. We notice that its value increases, as a consequence of EFI onset, due to the eVDF asymmetry in the parallel direction highlighted in Figure 1, panel (d). After saturation of the EFI, the heat flux starts decreasing, seemingly following the trend of the magnetic field oscillations.

In this paper, we analyze the role of solar wind expansion and selected collisionless processes (namely, the EFI in presence of two electron populations) in electron heat flux regulation. Our key finding is that, at least in our simplified scenario, the drift velocity between the electron populations plays a fundamental role in solar wind electron energy flux regulation. In our simulation, we observe a thermal velocity evolution that differs from DA expectations due to EFI onset. However, this does not seem to crucially affect heat flux regulation. The electron energy flux is instead regulated by the sudden drop of the electron relative drift velocity due to EFI onset. Solar wind plasma expansion indirectly contributes to the heat flux regulation by triggering the EFI instability. This work supports an indirect role of solar wind expansion in electron heat flux regulation, where expansion drives or modifies the evolution of heat flux instabilities.

We notice that the suprathermal (halo plus strahl) population could have been more realistically modeled with a non-Maxwellian distribution, composed, e.g., of a kappa distribution for the halo and of the realistic shape for the strahl described in Horaites et al. (2018). We expect that simulations initialized with such a VDF would not differ significantly from the one presented here, since EFI onset is driven by the core (see Figure 1, panel (b)), and therefore is not affected by the shape of the suprathermal eVDF. Also, drift velocity evolution of the suprathermal electrons is driven by the halo, since it has higher fractional density than the strahl (Maksimovic et al. 2005). The halo can be modeled with a kappa distribution, which is expected to behave similarly to a Maxwellian for processes that do not depend on the high-energy tails of the distribution.

Scime et al. (1994) study how the electron energy flux, obtained from Ulysses observations, varies between 1 and 5 au. They obtain a dependence with the heliocentric distance r of $\propto {r}^{-2.7}$, inconsistent with the expected variation from collisional processes alone and steeper than what is expected from expansion-driven evolution alone; they deduce that collisionless processes should contribute to heat flux regulation.

They confirm the validity of the phenomenological formula of Feldman et al. (1975), rewritten here as Equation (1), in explaining energy flux evolution with heliocentric distance, and confirm the Feldman et al. (1975) intuition that the collisionless mechanisms (i.e., instabilities) responsible for heat flux regulation are those that also control drift velocity evolution.

An obvious limitation of the present work is the fact that the global temperature profiles before EFI onset follow a DA evolution, rather than profiles with heliocentric distance consistent with solar wind observations—as shown, for example, in Landi et al. (2012). This is due to our simplified initial conditions, which do not allow for processes such as turbulence that we believe affect the observed radial profile of the temperature at least as much as collisions. Incidentally, we observe that introducing turbulent fluctuations in the initialization fields may result in a more accurate evolution of the eVDF. In Tang et al. (2018), for example, it is observed that introducing whistler wave turbulence into the kinetic Fokker–Planck transport equation allows one to obtain a core–halo–strahl eVDF and realistic radial dependences for key quantities starting from a core and a suprathermal electron component, with the latter being scattered by whistler fluctuations.

In future work, we will address this issue and analyze other scenarios, prone to the development of other instabilities potentially controlling energy flux evolution through drift velocity regulation, for a better understanding of the outer corona energy balance. The aim is to provide a reliable interpretative framework for observations of heat flux at different heliocentric distances along the same magnetic field line, which will become available with coordinated Parker Solar Probe (Fox et al. 2016) and Solar Orbiter (Mueller et al. 2013) observational campaigns.

M.E.I. thanks Rodrigo López for useful discussions. M.E.I.'s work is supported by an FWO postdoctoral fellowship. The simulations were performed on the supercomputer Marconi-Broadwell (Cineca, Italy) under a PRACE allocation. A.T. acknowledges support from NASA HSR grant No. 80NSSC18K1211. This research was supported in part by the NASA DRIVE HERMES project, grant No. 80NSSC20K0604.