Particle-in-cell simulations of the whistler heat-flux instability in the solar wind conditions

TL;DR

This study presents the first particle-in-cell simulations of the whistler heat-flux instability in solar wind conditions, revealing how it saturates and affects electron distributions without initial temperature anisotropies.

Contribution

It provides the first realistic PIC simulation of WHFI in solar wind conditions, demonstrating its saturation mechanism and impact on electron velocity distributions.

Findings

WHFI saturates at moderate amplitudes in simulations.

Enhanced whistler fluctuations cause electron pitch-angle scattering.

Strahl electron distribution becomes skewed with increasing energy.

Abstract

In collision-poor plasmas from space, e.g., solar wind or stellar outflows, the heat-flux carried by the strahl or beaming electrons is expected to be regulated by the self-generated instabilities. Recently, simultaneous field and particle observations have indeed revealed enhanced whistler-like fluctuations in the presence of counter-beaming populations of electrons, connecting these fluctuations to the whistler heat-flux instability (WHFI). This instability is predicted only for limited conditions of electron beam-plasmas, and was not captured in numerical simulations yet. In this letter we report the first simulations of WHFI in particle-in-cell (PIC) setups, realistic for the solar wind conditions, and without temperature gradients or anisotropies to trigger the instability in the initiation phase. The velocity distributions have a complex reaction to the enhanced whistler fluctuations conditioning the instability saturation by a decrease of the relative drifts combined with induced (effective) temperature anisotropies (heating the core electrons and pitch-angle and energy scattering the strahl). These results are in good agreement with a recent quasilinear approach, and support therefore a largely accepted belief that WHFI saturates at moderate amplitudes. In anti-sunward direction the strahl becomes skewed with a pitch-angle distribution decreasing in width as electron energy increases, that seems to be characteristic to self-generated whistlers and not to small-scale turbulence.

Tables

Table 1. Table 1: Initial plasma parameters for the simulation.

Parameter	Beam ($b$)	Core ($c$)	Protons ($p$)
$n_j/n_0$	0.05	0.95	1.0
$T_{j,\parallel }/T_{c,\parallel }$	6.0	1.0	1.0
${\beta }_{j,\parallel }$	18.0	3.0	3.0
$m_p/m_j$	1836	1836	1.0
$T_{j,⟂}/T_{j,\parallel }$	1.0	1.0	1.0
$U_j/v_A$	$40.0$	$-2.1$	0.0

Equations

$$f_e\left(v_{\perp },v_{\parallel }\right)=\frac{n_c}{n_e}{f}_c\left(v_{\perp },v_{\parallel }\right)+\frac{n_b}{n_e}{f}_b\left(v_{\perp },v_{\parallel }\right)$$

Full text

methods: numerical – plasmas – solar wind – waves – instabilities – interplanetary medium

††journal: ApJL

1 Motivations

Among the kinetic instabilities invoked in the self-regulation of solar wind properties the heat-flux instabilities, and in particular the whistler heat-flux instability (WHFI), are still the most controversial, though in the last decade an increased effort has been devoted to understanding their fundamental properties (Saito & Gary, 2007a; Pavan et al., 2013; Seough et al., 2015; Saeed et al., 2017a, b; Shaaban et al., 2018a, b, 2019; Lee et al., 2019) and find their signatures in the observation (Breneman et al., 2010; Gurgiolo et al., 2012; Wilson et al., 2013; Landi et al., 2014; Lacombe et al., 2014; Stansby et al., 2016; Tong et al., 2019b, a). The WHFI is triggered by the relative drift, $U=|U_{c}|+U_{b}$, of the counter-beaming electrons, a central population (summing up the core and halo electrons) here called generically core and denoted by subscript $c$, and the beam or strahl population (with subscript $b$) satisfying the zero net-current condition $n_{c}|U_{c}|=n_{b}U_{b}$, see Gary (1985) and refs therein. Conditions for the whistlers to be excited (resonantly) by the beaming electrons are however very restrained, namely, to a beaming velocity limited between two threshold values, roughly given by $\theta_{c}<U_{b}<\theta_{b}$, where $\theta_{c,b}$ are thermal velocities (Gary, 1985; Shaaban et al., 2018a, b). The quasi-stable states are expected in this case only for low drifts $U_{b}$ (or $U_{c}$), below the lower threshold (Gary et al., 1999a; Shaaban et al., 2018a) that seems to be confirmed by the observations (Gary et al., 1999a, b; Tong et al., 2018). Theoretically, whistlers may also satisfy resonance conditions with both electron populations, especially, for more energetic beams ($U_{b}>\theta_{b}$), but never develop, being heavily competed by the other faster growing modes, e.g., the electrostatic beam-plasma instabilities or the oblique instabilities (Gary & Saito, 2007; Saito & Gary, 2007a; Seough et al., 2015; Saeed et al., 2017a; Horaites et al., 2018; Lee et al., 2019; Vasko et al., 2019; Verscharen et al., 2019a). If the core electrons exhibit an important temperature anisotropy $T_{c,\perp}>T_{c,\parallel}$ the regime of WHFI may be significantly altered becoming specific to a standard whistler instability driven by temperature anisotropy, with lower thresholds and higher growth rates (Seough et al., 2015; Shaaban et al., 2018b).

The first investigations of WHFI have been stimulated by the observations suggesting a potential implication of whistlers in the regulation of suprathermal populations. If binary collisions are rare, in the solar wind the electron heat-flux is less than a conventional Spitzer-Härm level (Spitzer & Härm, 1953), and such constraint is attributed mainly to the wave-particle interactions (Bale et al., 2013). Moreover, with the expansion of the solar wind the electron halo shows a continuous build-up on the expense of strahl that lowers in intensity and undergoes a pitch-angle scattering (Maksimovic et al., 2005; Pagel et al., 2007; Gurgiolo et al., 2012; Berčič et al., 2019). In the absence of collisions an immediate explanation for these evolutions is offered by the small scale wave turbulence and/or the fluctuations self-generated by the instabilities. Higher plasma beta conditions stimulate the implication of self-generated instabilities in the regulation of suprathermal populations, in particular of the electron strahls (Pilipp et al., 1987; Crooker et al., 2003).

Theory and simulations have confirmed that whistler fluctuations, either predefined by a power spectrum decreasing monotonically with increasing frequency or self-generated by kinetic instabilities, can pitch-angle and energy scatter the suprathermal electrons and lead to asymmetric beaming-like distributions, broader (bulge) or decreasing (skewness) in pitch-angle at larger electron energies (Vocks & Mann, 2003; Vocks et al., 2005; Saito & Gary, 2007a, b; Seough et al., 2015). In particular, for the WHFI, quasilinear studies have also suggested a potential role in the limitation of the electron heat-flux, probably by the same mechanisms, which reduce the relative drift and induce effective anisotropies of electron populations (Gary & Feldman, 1977; Shaaban et al., 2019). However, a confirmation of these effects in simulations have not been reported yet. To our knowledge, numerical experiments have provided extended descriptions only for other different branches of heat-flux instabilities, e.g., electrostatic beam-plasma, firehose-like (Gary & Saito, 2007; Lee et al., 2019), or for the temperature anisotropy driven instabilities (Saito & Gary, 2007a; Seough et al., 2015).

Limiting conditions predicted for the WHFI (Gary, 1985; Shaaban et al., 2018a) and small amplitudes of the resulting fluctuations (Shaaban et al., 2019; Tong et al., 2019a) might have also prevented a direct detection in the observations, and leaded sometimes to contradictory correlations between plasma states and fluctuations (Scime et al., 2001). Clear evidences of WHFI in the solar wind have recently been provided by simultaneous electron and field measurements with a well established connection to the electron counter-beaming populations and their temperature anisotropy (Tong et al., 2019b, a). These observations confirm recent predictions that WHFI must be quenched by a slight anisotropy $T_{b,\parallel}\gtrsim T_{b,\perp}$ of the beam, but growth rates may significantly be increased by an opposite anisotropy of the core $T_{c,\perp}\gtrsim T_{c,\parallel}$ (Shaaban et al., 2018b).

This letter reports the first particle-in-cell (PIC) simulations of the WHFI, realistic for solar wind conditions. The characteristics of this instability (see above) impose serious limitations to describe it using simulations (with realistic parameters), requiring an immense amount of numerical resources, which are practically impossibles with standard approaches. Here we make use of an implicit PIC code developed by Markidis et al. (2010), able to resolve multiple temporal and spatial scales characteristic to the solar wind plasma dynamics (Verscharen et al., 2019b), without the strict limitations in time step and grid spacing imposed typically in explicit codes. Simulations capture the energy transfer between the electron core and beam populations, and correctly describe the saturation of WHFI via the relaxation of the velocity distributions. A quasilinear approach allows time variations of the moments of the distribution (e.g., drifts, temperatures), but implies only a single wave mode in the energy and momentum transfers (Shaaban et al., 2019). Instead the simulations enable quasilinear and nonlinear effects of multiple (concurrent, coupled) wave modes.

2 Particle-in-cell simulations

Our initial setup in Table 1 is intended to the solar wind plasma conditions (Maksimovic et al., 2005; Tong et al., 2018, 2019b), and, for simplicity, both counter-beaming electron populations (in a frame fixed to protons)

[TABLE]

are assumed Maxwellian distributed. Here $n_{c}$ and $n_{b}$ are the core and beam number densities, respectively, and $n_{e}\equiv n_{0}$ is the total number density of electrons, in a neutral plasma with zero charge $n_{e}\approx n_{i}$ and zero net current $n_{c}U_{c}+n_{b}U_{b}=0$, where $U_{b,c}$ are the corresponding drift velocities, and here the ions (subscript $i$) are assumed to be only protons.

We use an implicit one-dimensional PIC code, i.e., iPic3D (Markidis et al., 2010) with a high enough resolution to resolve the electron inertial length and the electron gyromotion. The spatial grid is composed of $n_{x}=1024$ cells, with $5000$ particles per species per grid. The box size is $L_{x}=16\,d_{i}$, then the cell size is $\Delta x=0.0156\,d_{i}$. Here $d_{i}=c/\omega_{pi}$ is the ion inertial length, with $\omega_{pi}=(4\pi n_{0}e^{2}/m_{p})^{1/2}$ the ion plasma frequency. The mass ratio is $m_{p}/m_{e}=1836$, and the plasma to gyro frequency ratio for ions is $\omega_{pi}/\Omega_{ci}=4390.07$, which implies that the Alfvén speed is $v_{A}=B_{0}/\sqrt{4\pi n_{p}m_{p}}=0.00023\,c$ and the plasma to gyro frequency of the electrons, $\omega_{pe}/\Omega_{ce}=102.48$, which are typical values encounter under solar wind conditions. The background magnetic field is set in the $x$ direction, $\mathbf{B}_{0}=B_{0}\hat{\mathbf{x}}$. The time step is $\Delta t=0.0375/\omega_{pi}$ and the simulations ran until $t_{\text{max}}=17560.265/\omega_{pi}$ or equivalently $t_{\text{max}}=4.0/\Omega_{ci}$. In terms of electron quantities, the time step used correspond to $\Delta t=0.016/\Omega_{ce}$ and the cell size is $\Delta x=0.7\,d_{e}$.

Fig. 1 shows (normalized) time variations, with the increase and saturation of the magnetic power ($W_{B}=\int\delta B^{2}/B_{0}^{2}\,dx$) of the enhanced fluctuations, and the relaxation of the main moments of electron velocity distributions which continues after the instability saturation. The rapid growth of $W_{B}$ corresponding to the excitation of WHFI in the early stage of the simulations slows down close to $\Omega_{i}t\approx 3$, and shows then a slower increase up to the end of the simulation, $\Omega_{p}t=4.0$. The entire period of the simulation can be identified in this case as characteristic to a pure WHFI. In order to identify this interval of pure whistler-like fluctuations we have used the fast Fourier transforms in space of the transverse magnetic fluctuations ($|\text{FFT}(B_{y}-iB_{z})|^{2}$), which are displayed in Fig. 2. In this interval only the intense power of the WHFI corresponding to positive wave numbers are present, see also Fig. 3.

Plasma beta parameters $\beta_{c,b}\equiv 8\pi n_{0}T_{c,b}/B_{0}^{2}$ plotted in Fig. 1 are defined with total number density $n_{0}$ and reflect therefore the variations of the corresponding temperatures $T_{c,b}$, in parallel (red) and perpendicular (blue) directions with respect to the background magnetic field. Initially isotropic, i.e., $\beta_{c,\parallel}(0)=\beta_{c,\perp}(0)$, the core temperatures are subjected to parallel cooling (red) and small perpendicular heating (blue) by a resonant cyclotron interaction of whistlers with the cooler electrons from the core. Beaming electrons have also isotropic temperatures at the beginning, i.e., $\beta_{b,\parallel}(0)=\beta_{b,\perp}(0)$, but their pitch-angle and energy scattering induces an opposite anisotropy. Consequently, at the saturation the core exhibits an excess of perpendicular temperature, i.e. $\beta_{c,\perp}(t^{*})>\beta_{c,\parallel}(t^{*})$, while the beam shows an excess of parallel temperature, i.e. $\beta_{b,\perp}(t^{*})<\beta_{b,\parallel}(t^{*})$. Both the linear theory of WHFI (Shaaban et al., 2018b) and the observations (Tong et al., 2019b) suggest indeed that this instability is inhibited by such a temperature anisotropy of the beam. These results are also in good agreement, at least at qualitative level and for the same time scales, with the QL evolutions predicted by theory, see, for instance, Fig. 9 (middle panels) and Fig. 10 in Shaaban et al. (2019). However, for such comparison we have to keep in mind that initial conditions, like the number of particles used in PIC simulations, or the initial level of the electromagnetic fluctuations in the QL theory, or both are crucial for the onset time of the instability (López & Yoon, 2018).

Right panels in Fig. 1 show the time relaxation of the electron heat-flux $q_{e}=m_{e}/2\int dv\,v_{x}vf_{e}$ (normalized by $q_{\text{max}}=3n_{0}T_{\parallel c}\alpha_{c,\parallel}/\sqrt{2}$, where $\alpha_{c,\parallel}=\sqrt{k_{B}T_{\parallel c}/m_{e}}$ is the thermal speed), and the core and beam (normalized) drift velocities ($U_{c,b}/v_{A}$). In the time interval relevant for the WHFI the heat flux and (counter-)drifts are only partially relaxed, showing similar reductions of about $25\%$ or $30\%$ of initial magnitude. We can state that the relaxation of relative drift velocities, i.e. $U_{c,b}(t_{\text{max}})\approx~{}0.67~{}U_{c,b}(0)$, is slowed down by a concurrent effect of the enhanced fluctuations, which interact with the electrons and induce opposite temperature anisotropies in the core and beam populations.

In Fig. 3 we plot the normalized power spectra for the initial stage of the simulation, i.e. $0<\Omega_{i}t<2.24$, to guarantee we are capturing the linear stage of the WHF instability and to have a fair comparison with the linear dispersion relation at $\Omega_{i}t=0$. The spectra is obtained from $|\text{FFT}(B_{y}-iB_{z})|^{2}$, where here the FFT is computed in space and time (then normalized to the maximum value of the spectra). By doing so, we are able to separate the contribution of left-handed (LH) and RH circularly polarized modes. Thus, for $\omega>0$ we observe the RH contribution with positive and negative helicity, $k>0$ and $k<0$, respectively, see Saeed et al. (2017a) for details. Most of the magnetic power is concentrated in the part with $\omega>0$ and $k>0$, corresponding to the RH unstable modes with positive helicity, confirming the linear theory predictions (top panel) for a RH WHFI. Moreover, we observe a very low intensity in the negative wave-numbers part of the spectrum, but those are modes with negative helicity and are damped, according to linear calculations. The other combinations do not show any significant power (not shown here). Here we can state that our PIC simulations are capable to capture the low intensity whistler fluctuations associated with the WHFI that can develop only for $\omega>0$ and $k>0$. The dispersion shown by the simulated fluctuations is not an instantaneous picture, but rather a cumulative contribution of fluctuations in the entire period under consideration, when macroscopic plasma values evolve from the initial condition. However, the unstable wavenumber interval does not change much with the beaming speed (see, Figs. 1 and 2 in Shaaban et al. (2019)) to explain the broad wavenumber spectra, which may probably result from the small error in the energy conservation in this simulation. Reducing the time step or increasing the number of particles per grid cell would help to improve the energy conservation and therefore obtain more accurate results, but more computational resources will be needed.

Fig. 4 presents the velocity distribution $f(v_{x},v_{y})$ at different relevant stages of the simulation $\Omega_{i}t=0.0$, $2.24$, and $3.92$, for the total electron population (upper panels) and the beam component (lower panels), as well as the reduced distribution (integrated along $v_{y}$) $f_{e}(v_{x})$ at the initial and (almost) final stages of the simulation, i.e. $\Omega_{i}t=0.0$, and $3.92$ (right panel). In order to highlight the deformation of the electron components in the distribution, we have carefully selected three particular contours, as indicated with dotted lines, at $2\times 10^{-4},\,3\times 10^{-4},\,8\times 10^{-4}$. It is clear that the highest contour of level $8\times 10^{-4}$ (black dotted) becomes more symmetric at the end of the simulation, showing also a slight increase of temperature anisotropy of the core population $T_{c,\perp}>T_{c,\parallel}$ (upper panels) and giving an indication for the relaxation of the drift velocities. In the case of the beam (lower panels) this contour shows the behavior observed in Fig. 1, a generation of parallel anisotropy. At later stages of the simulation contours of lower level (white), e.g., $3\times 10^{-4}$ and $\,2\times 10^{-4}$, are slightly different than those at the initial state, and specifically show an asymmetric skewness of less scattered particles (pitch-angle scattering of the beam decreasing in parallel direction as electron energy increases). Moreover, a lower (relaxed) but still finite drift velocity is more obviously shown by the reduced distributions in the right panel. The reduced distributions confirm the previous description that the initial drift velocities (black line) are regulated by the enhanced WHF fluctuations and the electron components ended up with small but finite relative drift velocities (blue line). Moreover, at final stage, i.e. $\Omega_{i}t=3.92$, the reduced eVDF $f_{e}(v_{x})$ shows the formation of a small, but still noticeable, “shoulder” in the parallel direction for the beam component, suggesting already that not all beaming electrons are scattered by the enhanced fluctuations, an hypothesis confirmed by the results in Fig. 5.

Finally, Fig. 5 shows the departures of the distributions from the initial condition $\delta f_{j}(t)=f_{j}(t)-f_{j}(0)$, for core (top panels) and beam (bottom panels) electrons, and for the two relevant moments $\Omega_{i}t=2.24$ (left) and $\Omega_{i}t=3.92$ (right). Red contours show $\delta f_{j}>0$ with an abundance of scattered electrons, while blue contours $\delta f_{j}<0$ mark the electron loss. Here we can see how different electron components are scattered (or not) by the whistler waves. Correlating with Fig. 1, the diffusion of core electrons occurs under the effect of whistlers which interact resonantly with the electrons with $v_{x}<|U_{c}|$, while the instability itself is (resonantly) triggered by the beaming electrons with $v_{x}<U_{b}$, cooling them down in perpendicular direction and increasing their effective temperature (or kinetic energy) in parallel direction. The lighter blue color population at higher energies in Fig. 5, right-lower panel, indicates those electrons less scattered by whistlers, and corresponds to the small shoulder (or small plateau) showed in Fig. 4. In time this population is naturally reduced leading to a lower pitch-angular width that becomes however prominent due to a concomitant decrease of the drift.

3 Summary

In this Letter we have provided a detailed description of the whistler heat-flux instability (WHFI) using an implicit one-dimensional PIC simulation. The instability is triggered by the relative drift of the counter-beaming electron populations, without temperature gradients or temperature anisotropies. The initial stage of the simulation is characterized by a rapid growth of the magnetic energy density, corresponding to the excitation of the WHFI, then corroborated by the spectral analysis with a good agreement with linear theory. The enhanced whistler fluctuations interact with both electron components, reducing the relative drift ($\sim 30\%$) and inducing (effective) temperature anisotropies, i.e., an excess of perpendicular temperature for the core and excess of parallel temperature for the beam. The interplay of temperature anisotropies is in good agreement with a recent QL approach (Shaaban et al., 2019) (although the drift relaxation is less significant under the effect of a single mode in QL theory), and support therefore a largely accepted belief that WHFI saturates at typically small amplitudes.

In anti-sunward direction the strahl becomes skewed with a pitch-angle distribution decreasing in width as electron energy increases, that seems to be characteristic to self-generated whistlers and not to small-scale turbulence. However, this skewness (a decreasing pitch-angle distribution with increasing energy) is shown only by the lower levels (white dashed contours) with a lower contribution to the moments of the distribution, and implicitly to the effective temperature anisotropy, which is reduced.

Future refinements to clarify the nonlinear evolution of this instability need to be taken with caution and eventually using new codes that conserve much better the energy.

These results were obtained in the framework of the projects SCHL 201/35-1 (DFG-German Research Foundation), GOA/2015-014 (KU Leuven), G0A2316N (FWO-Vlaanderen), and C 90347 (ESA Prodex 9). S.M.S. acknowledges support by a FWO Postdoctoral Fellowship, Grant No. 12Z6218N. P.H.Y. acknowledge BK21 Plus program from NRF Korea. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI. We acknowledge fruitful discussions at the meeting of international team on Kappa Distributions hosted by ISSI-Bern.