Nonlinear evolution of the whistler heat flux instability

TL;DR

This study uses particle-in-cell simulations to analyze the nonlinear evolution of the whistler heat flux instability in solar wind conditions, revealing that it produces waves that do not significantly suppress electron heat flux.

Contribution

The paper demonstrates through simulations that the whistler heat flux instability generates waves consistent with observations but does not effectively regulate electron heat flux in the solar wind.

Findings

Simulations produce whistler waves matching spacecraft amplitudes.

Wave amplitude correlates with electron heat flux and plasma beta.

Whistler waves do not significantly suppress electron heat flux.

Abstract

We use the one-dimensional TRISTAN-MP particle-in-cell code to model the nonlinear evolution of the whistler heat flux instability that was proposed by Gary et al. (1999, 2000) to regulate the electron heat flux in the solar wind and astrophysical plasmas. The simulations are initialized with electron velocity distribution functions typical for the solar wind. We perform a set of simulations at various initial values of the electron heat flux and $\beta_{e}$. The simulations show that parallel whistler waves produced by the whistler heat flux instability saturate at amplitudes consistent with the spacecraft measurements. The simulations also reproduce the correlations of the saturated whistler wave amplitude with the electron heat flux and $\beta_{e}$ revealed in the spacecraft measurements. The major result is that parallel whistler waves produced by the whistler heat flux instability do not significantly suppress the electron heat flux. The presented simulations indicate that coherent parallel whistler waves observed in the solar wind are unlikely to regulate the heat flux of solar wind electrons.

Tables

Table 1. Table 1: Parameters for the first set of simulations performed at β c = 1 subscript 𝛽 𝑐 1 \beta_{c}=1 . Indicated are initial values of u c / v A subscript 𝑢 𝑐 subscript 𝑣 𝐴 u_{c}/v_{A} and q e / q 0 subscript 𝑞 𝑒 subscript 𝑞 0 q_{e}/q_{0} . In all these simulations n c / n e = 0.85 subscript 𝑛 𝑐 subscript 𝑛 𝑒 0.85 n_{c}/n_{e}=0.85 , T h / T c = 10 subscript 𝑇 ℎ subscript 𝑇 𝑐 10 T_{h}/T_{c}=10 and ω p ​ e / Ω e ≈ 12.3 subscript 𝜔 𝑝 𝑒 subscript Ω 𝑒 12.3 \omega_{pe}/\Omega_{e}\approx 12.3 . In Figure 2 the initial conditions for the first set of simulations are indicated by the blue squares.

$\#$	$-u_c/v_A$	$q_e/q_0$
1	3.5	0.3
2	5.5	0.45
3	7.5	0.65
4	9.0	0.8
5	14.0	1.3
6	18.5	1.9

Table 2. Table 2: Parameters for the second set of simulations performed at the same initial heat flux value, q e / q 0 = 0.45 subscript 𝑞 𝑒 subscript 𝑞 0 0.45 q_{e}/q_{0}=0.45 . Indicated are initial values of β c subscript 𝛽 𝑐 \beta_{c} , u c / v A subscript 𝑢 𝑐 subscript 𝑣 𝐴 u_{c}/v_{A} and ω p ​ e / Ω e subscript 𝜔 𝑝 𝑒 subscript Ω 𝑒 \omega_{pe}/\Omega_{e} . In all simulations n c / n e = 0.85 subscript 𝑛 𝑐 subscript 𝑛 𝑒 0.85 n_{c}/n_{e}=0.85 and T h / T c = 10 subscript 𝑇 ℎ subscript 𝑇 𝑐 10 T_{h}/T_{c}=10 . In Figure 2 the initial conditions for the second set of simulations are indicated by four squares (red, green, blue and black) at q e / q 0 = 0.45 subscript 𝑞 𝑒 subscript 𝑞 0 0.45 q_{e}/q_{0}=0.45 .

$\#$	${\beta }_c$	$-u_c/v_A$	${\omega }_{pe}/{\Omega }_e$
1	0.4	3.5	7.8
2	1	5.5	12.3
3	2	7.8	17.3
4	3	9.5	21.2

Equations

$$f_{\alpha }$$

$$\frac{k^2{c}^2+{\omega }_{pi}^2}{{\omega }^2}\approx \frac{{\omega }_{pe}^2}{{\omega }^2}\sum _{\alpha =c,h}\frac{n_{\alpha }}{n_e}\frac{(\omega -k{u}_{\alpha })}{k{v}_{\alpha }}Z\left(\frac{\omega -k{u}_{\alpha }-{\Omega }_e}{k{v}_{\alpha }}\right)$$

$$\max (B_w)/B_0\sim A({\gamma }_{\max }/{\Omega }_e{)}^{\alpha }$$

Full text

solar wind — heat flux — whistler waves

††journal: The Astrophysical Journal

1 Introduction

The electron heat flux in collisionless or weakly collisional plasma is generally not described by the collisional Spitzer-Härm law (Spitzer & Härm, 1953). The heat flux suppression below the collisional value was demonstrated by direct in-situ measurements in the solar wind (Feldman et al., 1975; Scime et al., 1994; Crooker et al., 2003; Bale et al., 2013) and remote observations of the temperature distribution of a hot gas in galaxy clusters (Cowie & McKee, 1977; Bertschinger & Meiksin, 1986; Fabian, 1994; Zakamska & Narayan, 2003; Wagh et al., 2014; Fang et al., 2018). The Spitzer-Härm law was hypothesized to be inadequate for describing the heat conduction in the solar corona (Scudder, 1992; Landi & Pantellini, 2001; Dorelli & Scudder, 2003) and that was recently confirmed in the analysis of coronal loop oscillations (Wang et al., 2015). The observations in the solar wind were interpreted in terms of the heat flux regulation by wave-particle interactions (Feldman et al., 1975, 1976; Gary & Feldman, 1977; Scime et al., 1994; Crooker et al., 2003), but the wave activity potentially regulating the electron heat flux is still under debate (e.g., Gary et al., 1999; Scime et al., 2001; Pagel et al., 2007; Roberg-Clark et al., 2016; Tong et al., 2018; Roberg-Clark et al., 2018b, a; Komarov et al., 2018; Vasko et al., 2019; Verscharen et al., 2019). The alternative view is that wave-particle interactions may be not necessary to explain the observed heat flux values in the solar wind (Landi et al., 2012, 2014; Horaites et al., 2015).

The mechanism of the electron heat flux regulation is intimately related to kinetic features of the electron velocity distribution function (VDF). The electron VDF in a slow solar wind is often adequately described by a combination of bi-Maxwellian thermal dense core and suprathermal tenuous halo populations (Feldman et al., 1975, 1976; Maksimovic et al., 1997; Tong et al., 2019b). In the plasma rest frame, the core and halo populations are streaming along the background magnetic field, the net electron current is approximately zero, while the electron heat flux is anti-sunward and carried predominantly by halo electrons (Feldman et al., 1975, 1976; Scime et al., 1994; Tong et al., 2019b). In what follows we focus on a slow solar wind and do not consider effects of the strahl population (anti-sunward beam-like population) frequently observed in a fast solar wind and carrying a significant portion of the electron heat flux (e.g., Rosenbauer et al., 1977; Pilipp et al., 1987; Štverák et al., 2009).

At sufficiently large heat flux values electrons are capable of generating whistler waves propagating quasi-parallel to the background magnetic field (Gary et al., 1975, 1994). Spacecraft measurements of the electron heat flux values below a bound dependent on $\beta_{e}$ (electron beta parameter) were interpreted in terms of the heat flux regulation by this so-called whistler heat flux instability (WHFI) (Feldman et al., 1976; Gary & Feldman, 1977; Gary et al., 1999). The argument behind that hypothesis was that the heat flux bound dependent on $\beta_{e}$ is similar to the marginal stability threshold of the WHFI. Gary & Li (2000) extrapolated the hypothesis and proposed that the WHFI may regulate the electron heat flux in high-$\beta_{e}$ astrophysical plasma. Pistinner & Eichler (1998) called into question that hypothesis, but no detailed analysis of the nonlinear stage of the classical WHFI was performed to support their arguments. Scime et al. (2001) questioned that the WHFI regulates the electron heat flux in the solar wind based on Ulysses measurements.

The simultaneous wave and particle measurements have recently shown that, though intermittently, whistler waves are indeed present in the solar wind and they are predominantly quasi-parallel that is consistent with the WHFI scenario (Lacombe et al., 2014; Stansby et al., 2016; Kajdič et al., 2016; Tong et al., 2019b, a). Tong et al. (2019b) have demonstrated for several events that the whistler waves were indeed generated locally by the WHFI. The extensive statistical analysis by Tong et al. (2019a) has shown that in the solar wind whistler waves have amplitudes typically less than a few hundredths of the background magnetic field. The amplitude is positively correlated with $\beta_{e}$ and generally with the electron heat flux. Though the WHFI is indeed operating in the solar wind, there has been no analysis that would establish whistler wave properties in the nonlinear stage of the WHFI and clarify whistler wave effects on the electron heat flux evolution.

In this paper, we present Particle-In-Cell (PIC) simulations of the WHFI using the fully relativistic massively parallelized TRISTAN-MP code (Spitkovsky, 2008; Park et al., 2015). The simulations are restricted to parallel whistler waves, because they have the largest growth rates according to the linear theory of the WHFI (Gary et al., 1975, 1994), and initialized with proton and electron VDFs typical for the solar wind. The simulations demonstrate that in a uniform plasma parallel whistler waves produced by the WHFI are incapable of regulating the electron heat flux. The whistler waves are shown to saturate at amplitudes consistent with the recent spacecraft measurements. The paper is organized as follows. We summarize the results of the WHFI linear theory in Section 2, then present results of the PIC simulations in Section 3, discuss and summarize our conclusions in Sections 4 and 5.

2 WHFI linear theory

In the linear theory of the original WHFI (Gary et al., 1975), the electron VDF consists of Maxwellian core and halo populations streaming along the background magnetic field, $f_{e}(v_{||},v_{\perp})=f_{c}(v_{||},v_{\perp})+f_{h}(v_{||},v_{\perp})$, where

[TABLE]

and $n_{\alpha}$, $u_{\alpha}$ and $T_{\alpha}$ denote density, bulk velocity and temperature of the core ($\alpha=c$) and halo ($\alpha=h$) populations, $m_{e}$ is the electron mass. The electron current in the plasma rest frame is assumed to be zero, $n_{c}u_{c}+n_{h}u_{h}=0$. The electron heat flux associated with the core plus halo electron VDF, $q_{e}=-0.5\;n_{c}u_{c}\;(5T_{h}+m_{e}u_{h}^{2}-5T_{c}-m_{e}u_{c}^{2})$, is often normalized to the free-streaming heat flux value $q_{0}=1.5\;n_{e}T_{e}(2T_{e}/m_{e})^{1/2}$, where $n_{e}=n_{c}+n_{h}$ is the total electron density, $T_{e}=(n_{c}T_{c}+n_{h}T_{h})/n_{e}$ is typical macroscopic electron temperature (e.g., Gary et al., 1999). The normalized electron heat flux $q_{e}/q_{0}$ depends on $n_{c}/n_{e}$, $T_{h}/T_{c}$ and $u_{c}\sqrt{m_{e}/T_{c}}$. At realistic parameters in the solar wind, thermal protons do not interact resonantly with whistler waves unstable to the WHFI (e.g., Gary et al., 1994; Gary & Li, 2000), therefore the contribution of protons to the whistler wave linear dispersion relation is considered in the frame of the cold fluid approach.

At $\omega_{pe}\gg\Omega_{e}$ the frequency and the linear growth rate of a parallel whistler wave are determined by the following dispersion relation (e.g., Mikhailovskii, 1974)

[TABLE]

where $\omega$ and $k$ are the whistler wave frequency and wavenumber, $\omega_{pe}$ and $\omega_{pi}$ are electron and proton plasma frequencies, $\Omega_{e}$ is the electron cyclotron frequency, $Z(\xi)=(2\pi)^{-1/2}\int_{-\infty}^{+\infty}dx\;e^{-x^{2}/2}/(x-\xi)$ is the plasma dispersion function, $v_{\alpha}=(T_{\alpha}/m_{e})^{1/2}$ denotes thermal velocity of the core and halo populations. We solve the dispersion relation in the limit of a weak instability via the standard techniques (e.g., Mikhailovskii, 1974). The linear growth rate $\gamma$ normalized to $\Omega_{e}$ depends on $kc/\omega_{pe}$, $n_{c}/n_{e}$, $T_{h}/T_{c}$, $\beta_{c}=8\pi n_{c}T_{c}/B_{0}^{2}$ and $u_{c}/v_{A}$, where $v_{A}=B_{0}/(4\pi n_{e}m_{p})^{1/2}$ is the Alfvén velocity, $B_{0}$ is the background magnetic field, $m_{p}$ is the proton mass. In what follows, we present results of the linear stability analysis of the WHFI for $n_{c}/n_{e}=0.85$ and $T_{h}=10\,T_{c}$ if not stated otherwise.

Figure 1 presents results of the linear stability analysis at $\beta_{c}=1$ and various $u_{c}/v_{A}$ in the range from $-1$ to $-19$. Panel (a) shows that the whistler wave dispersion curves are dependent on $u_{c}/v_{A}$ and generally different from the whistler dispersion curve in a cold plasma, $\omega=\Omega_{e}k^{2}c^{2}/(k^{2}c^{2}+\omega_{pe}^{2})$ (e.g., Mikhailovskii, 1974). Panel (b) presents the whistler wave growth rates at various $u_{c}/v_{A}$ and demonstrates that the growth rates reach maximum values of about 0.01 $\Omega_{e}$ at $kc/\omega_{pe}\approx 0.3$. The wavelength of the fastest growing whistler waves is about 20 $c/\omega_{pe}$. We have performed similar growth rate computations at $\beta_{c}=0.4,2$ and $3$ and determined parameters of the fastest growing whistler waves.

Figure 2 presents parameters of the fastest growing whistler waves in dependence on $u_{c}/v_{A}$ at various $\beta_{c}$. Panel (a) shows that whistler waves at lower frequencies are preferentially unstable at larger $\beta_{c}$. Panel (b) shows that the growth rate of the fastest growing whistler wave is larger at larger $\beta_{c}$. Panel (c) presents the normalized heat flux $q_{e}/q_{0}$ as a function of $u_{c}/v_{A}$ for the further references. The normalized heat flux $q_{e}/q_{0}$ varies with $u_{c}/v_{A}$ and $\beta_{c}$, because it depends on $m_{e}u_{c}^{2}/T_{c}\sim u_{c}^{2}/v_{A}^{2}\beta_{c}$. Panels (b) and (c) show that at any given $\beta_{c}$, the growth rate of the fastest growing whistler wave is a non-monotonous function of $u_{c}/v_{A}$ or, equivalently, $q_{e}/q_{0}$ (in accordance with Gary, 1985).

In addition to the linear theory results for $T_{h}/T_{c}=10$, Figure 2 also presents parameters of the fastest growing whistler waves in a plasma with $\beta_{c}=1$ and $T_{h}/T_{c}=4$. The comparison of the growth rates computed at $T_{h}/T_{c}=4$ and 10 shows that the smaller halo temperature results in more than five times smaller maximum growth rate. The specific values of $u_{c}/v_{A}$ and $q_{e}/q_{0}$ at which the growth rate reaches maximum depend on $T_{h}/T_{c}$ and $n_{c}/n_{e}$. For example, panels (b) and (c) show that at $T_{h}/T_{c}=10$, the growth rate reaches maximum at $u_{c}\approx-10\,v_{A}$ and $q_{e}/q_{0}\approx 0.95$, while at $T_{h}/T_{c}=4$, it reaches maximum at $u_{c}\approx-6\,v_{A}$ and $q_{e}/q_{0}\approx 0.5$.

3 PIC simulations

We perform PIC simulations of the WHFI to establish the effect of the whistler waves on the electron heat flux evolution and determine the amplitude of whistler waves in the nonlinear stage. We use the TRISTAN-MP code (Spitkovsky, 2008; Park et al., 2015) with both ions and electrons treated as particles at realistic mass ratio. We restrict the analysis to 1D simulations with a single spatial coordinate $x$ along the background magnetic field, hence only strictly parallel whistler waves are modelled. The code is initialized with core and halo electron VDFs given by Maxwell-Jüttner distribution functions, which reduce to Maxwell distributions in a non-relativistic limit (that is the case in our simulations). The initial proton VDF is a Maxwell-Jüttner distribution function with the proton temperature equal to the core electron temperature (that is realistic in the solar wind at 1 AU, see Artemyev et al., 2018). Due to computational limitations, the simulations are performed at $\omega_{pe}/\Omega_{e}$ of the order of 10, while in the realistic solar wind this parameter is about ten times larger. At $\beta_{c}\sim 1$, this implies the temperature of the core population about 100 times larger than in the realistic solar wind (a few keV instead of 10 eV), because $\omega_{pe}/\Omega_{e}=\beta_{c}^{1/2}(m_{e}c^{2}/2T_{c})^{1/2}$. In fact, in all simulations $T_{c}=2$ keV, while the specific value of $\omega_{pe}/\Omega_{e}$ is determined by $\beta_{c}$. The length of the simulation box is $L_{x}=1296\,c/\omega_{pe}$ that is much larger than the typical wavelength $20\,c/\omega_{pe}$ of the fastest growing whistler waves (Figure 1). The temporal and spatial integration steps are 0.09 $\omega_{pe}^{-1}$ and 0.2 $c/\omega_{pe}$ that is adequate to reproduce the expected unstable whistler waves. The number of spatial cells is about 6480, the number of particles per cell is $4\cdot 10^{4}$, giving the total number of particles $\approx 2.6\cdot 10^{8}$.

In all simulations, we assume $n_{c}/n_{e}=0.85$ that is realistic for the solar wind (e.g., Maksimovic et al., 2005; Tong et al., 2019b). Typically in the solar wind $T_{h}/T_{c}\sim 4$ (e.g., Maksimovic et al., 2005; Tong et al., 2019b), but in all simulations we assume $T_{h}/T_{c}=10$ to have larger initial growth rates and thereby reduce the saturation time of the instability. We perform two set of simulations. In the first set, $\beta_{c}=1$ and $\omega_{pe}/\Omega_{e}\approx 12.3$, while initial values of $u_{c}/v_{A}$ and $q_{e}/q_{0}$ are given in Table 1 and indicated in Figure 2. In the second set, the initial electron heat flux is $q_{e}/q_{0}=0.45$ and simulations are performed at $\beta_{c}=0.4-3$. The corresponding initial values of $u_{c}/v_{A}$ and $\omega_{pe}/\Omega_{e}$ are given in Table 2 and indicated in Figure 2. The realistic parameters in the solar wind are $\beta_{c}\sim 1$, $|u_{c}|/v_{A}\lesssim 10$ and $q_{e}/q_{0}\lesssim 1$ (e.g., Tong et al., 2018, 2019b). Several simulations of the first set are performed at $u_{c}/v_{A}$ larger than typically observed in the solar wind (Table 1).

Figure 3 presents results of a simulation from the first set with $u_{c}=-9\,v_{A}$. Panel (a) presents temporal evolution of the magnetic field $B_{\perp}(x,t)=\left(B_{z}^{2}+B_{y}^{2}\right)^{1/2}$ perpendicular to the background magnetic field, demonstrating the growth of electromagnetic waves propagating parallel to the electron heat flux. Panel (b) presents the spectrum of the magnetic field $B_{y}$ around the center of the simulation box, demonstrating that the waves are at the central frequency of about 0.05 $\Omega_{e}$ and have a bandwidth comparable to the central frequency. Panel (c) presents waveforms of $B_{y}$ and $B_{z}$ around the center of the simulation box and shows that the amplitude of the waves is less than 0.1 $B_{0}$. The 90*∘* phase-shift between $B_{y}$ and $B_{z}$, meaning the right-hand polarization, and the wave propagation parallel the electron heat flux prove that the waves are the whistler waves produced by the WHFI. The central frequency of 0.05 $\Omega_{e}$ is consistent with the linear stability theory (Figure 2a). Panel (d) presents the temporal evolution of the whistler wave amplitude averaged over the simulation box, $B_{w}\equiv\left(L_{x}^{-1}\int B^{2}_{\perp}(x,t)\;dx\right)^{1/2}$. The instability saturates after about 2000 $\Omega_{e}^{-1}$ at the averaged whistler wave amplitudes of $0.03\,B_{0}$. Panel (e) presents the temporal evolution of the growth rate computed as $\gamma=B_{w}^{-1}dB_{w}/dt$ and shows that the initial growth rate $\gamma/\Omega_{e}\sim 0.01$ is consistent with the linear stability theory (Figure 2b).

Figure 4 clarifies the mechanism of the saturation of the WHFI. Panel (a) presents the initial electron VDF, while panel (b) presents the electron VDF after 2000 $\Omega_{e}^{-1}$ that is after the saturation of the instability. Panel (c) presents the difference between electron VDFs at the initial moment and after the saturation of the instability. The scattering by the whistler waves results in redistribution of the phase space density of electrons with $v_{||}\sim-0.3\,c$. Namely, electrons with initial parallel velocities $v_{||}\gtrsim-0.3\,c$ are transported to parallel velocities $v_{||}\lesssim-0.3\,c$ and vice versa. These are actually electrons in the first normal cyclotron resonance, $v_{||}=v_{R}\approx(\omega-\Omega_{e})/k$ (e.g., Shklyar & Matsumoto, 2009; Soto?Chavez et al., 2014), because at $\omega_{pe}/\Omega_{e}\approx 9$, $\omega\sim 0.05\,\Omega_{e}$ and $kc/\omega_{pe}\approx 0.3$, we have $v_{R}\approx(\omega-\Omega_{e})/k\approx-0.3\,c$. The resonant electrons provide the energy for the whistler wave growth, while the gradual scattering of these electrons by the whistler waves results in saturation of the instability via formation of a “plateau” on the electron VDF (e.g., Kennel & Engelmann, 1966; Sagdeev & Galeev, 1969; Karpman, 1974). Panel (c) shows that the scattering is efficient for electrons with $|v_{||}-v_{R}|\lesssim\Delta v_{R}$, where $\Delta v_{R}\sim 0.2\,c$. The finite resonance width $\Delta v_{R}$ is due to the finite bandwidth of the whistler waves, $\Delta v_{R}\approx\Omega_{e}\Delta k/k^{2}$, where $\Delta k$ is the bandwidth of the whistler waves in the wavenumber space. Adopting the cold dispersion relation for estimates, $\omega\approx\Omega_{e}k^{2}c^{2}/\omega_{pe}^{2}$, we derive $\Delta v_{R}\approx|v_{R}|\;\Delta\omega/2\omega\sim 0.2\,c$, because the bandwidth $\Delta\omega$ is comparable to the central frequency $\omega$ (Figure 3b).

Figure 5 summarizes results of the first set of simulations performed at $\beta_{c}=1$ and various initial $u_{c}/v_{A}$ and $q_{e}/q_{0}$ (Table 1). Panel (a) shows the temporal evolution of the averaged whistler wave amplitudes and demonstrates that the amplitude in the saturation stage depends on the initial heat flux value. The whistler waves saturate at amplitudes of about $0.02-0.04\,B_{0}$. Panels (b) and (c) show that the scattering of electrons by the whistler waves result in variations of the electron heat flux and $\beta_{c}$ by less than a few percent. The variation of $\beta_{c}$ is mostly due to the numerical heating of cold electrons that is recognized by the monotonous growth of $\beta_{c}$. Though the initial heat flux values in the first set of simulations are in the range from 0.3 to 1.9 (Table 1), the heat flux values in the saturation stage differ from the initial values by less than a few percent. These results do not support the scenario that parallel whistler waves produced by the WHFI suppress the electron heat flux below the marginal stability threshold dependent on $\beta_{c}$ (Gary et al., 1999; Gary & Li, 2000). According to that scenario, the heat flux values in the saturation stage would decrease down to approximately the same value in all the simulations, because the marginal stability threshold depends only on $\beta_{c}$. Yet, Figure 5 shows that the heat flux does not decrease significantly, remaining almost constant for all the simulations in this set with different initial heat flux values (Table 1).

Figure 6 summarizes results of the second set of simulations performed at the same initial heat flux $q_{e}/q_{0}=0.45$ and various $\beta_{c}$ (Table 2). Panel (a) presents the temporal evolution of the averaged whistler wave amplitude and demonstrates that the amplitude in the saturation stage depends on $\beta_{c}$. For $\beta_{c}=0.4-3$, the whistler waves saturate at amplitudes of $0.01-0.05\,B_{0}$. We should mention that in the simulation at $\beta_{c}=0.4$ the computations were performed up to $3500\,\Omega_{e}^{-1}$, because it takes longer for the instability to saturate at low $\beta_{c}$ due to smaller initial growth rates (Figure 2b). Panels (b) and (c) show that the scattering of electrons by the whistler waves results in variation of the heat flux and $\beta_{c}$ by less than a few percent. Again, these results do not support the scenario that parallel whistler waves produced by the WHFI suppress the electron heat flux below the marginal stability threshold (Gary et al., 1999; Gary & Li, 2000). The heat flux values in the simulations with various $\beta_{c}$ differ by less than a few percent, while according to that scenario the heat flux values in the saturation stage would have been different, because the marginal stability threshold depends on $\beta_{c}$.

Figure 7 presents correlations between the whistler wave amplitude in the saturation stage ${\rm max}(B_{w})/B_{0}$ and initial parameters of the electron VDF. Panel (a) presents ${\rm max}(B_{w})/B_{0}$ found in the first set of simulations versus initial $q_{e}/q_{0}$ and demonstrates that these quantities are positively correlated, though the correlation switches to negative at $q_{e}/q_{0}\gtrsim 0.8$. Panel (b) presents ${\rm max}(B_{w})/B_{0}$ found in the second set of simulations versus $\beta_{c}$ and demonstrates a positive correlation between these quantities. Panels (a) and (b) also present the linear growth rates $\gamma_{\rm max}/\Omega_{e}$ of the fastest growing whistler waves corresponding to the performed simulations (Figure 2). There is a clear positive correlation between the whistler wave amplitude in the saturation stage and the initial maximum growth rate. The results of the simulations are well fitted to

[TABLE]

with $A\sim 0.57$ and $\alpha\sim 0.64$. The specific values of $A$ and $\alpha$ are likely dependent on $n_{c}/n_{e}$ and $T_{h}/T_{c}$ that were fixed in our simulations, but the positive correlation between ${\rm max}(B_{w})/B_{0}$ and $\gamma_{\rm max}/\Omega_{e}$ is expected to hold.

4 Discussion

The electron heat flux is one of the leading terms in the electron energy balance in the expanding solar wind (e.g., Cranmer et al., 2009; Stveraḱ et al., 2015), that is why kinetic processes regulating the electron heat flux are of fundamental importance in the solar wind physics. In this paper, we have presented results of 1D PIC simulations demonstrating that in a uniform plasma, the WHFI is incapable of regulating the electron heat flux. In other words, the presented simulations do not support the hypothesis proposed by Gary et al. (1999) and Gary & Li (2000) that whistler waves produced by the WHFI suppress the electron heat flux below the $\beta_{e}-$dependent marginal stability threshold. Though that hypothesis was previously questioned based on theoretical estimates (Pistinner & Eichler, 1998) and Ulysses measurements (Scime et al., 2001), no detailed analysis of the nonlinear stage of the WHFI was provided to support those statements. We underline that our conclusions are based on the simulations restricted to parallel whistler waves. Though parallel whistler waves have the largest growth rates (Gary et al., 1975, 1994), 2D PIC simulations should be performed in the future to clarify effects of slightly oblique whistler waves produced by the WHFI.

The recent 1D PIC simulations by Roberg-Clark et al. (2016) have also demonstrated that strictly parallel whistler waves are incapable of regulating the electron heat flux in a collisionless plasma. However, the initial electron VDF set by Roberg-Clark et al. (2016) is an asymptotic solution of the collisional Fokker-Plank equation (e.g., Pistinner & Eichler, 1998), that is very different from electron VDFs typical for the solar wind (Feldman et al., 1975; Maksimovic et al., 1997; Tong et al., 2019b). In contrast to Roberg-Clark et al. (2016), we have presented 1D PIC simulations of the classical WHFI with initial electron VDFs more relevant for the solar wind, so that the results of the simulations can be compared to in-situ measurements. We underline that in the realistic solar wind, the halo population is better fitted to $\kappa-$distribution with $\kappa\sim 4$, and typically $T_{h}/T_{c}\sim 4$ that is lower than used in our simulations (e.g., Maksimovic et al., 2005; Tong et al., 2019b). Such a choice of distribution function and temperature ratio would affect the quantitative results for linear growth rates, and, consequently, the saturated wave amplitudes (see Eq. 1)), as we discuss below. But it does not affect the qualitative results for the WHFI, still producing low-amplitude quasi-parallel whistler waves.

The presented simulations show that whistler waves saturate at amplitudes of about $0.01-0.05\,B_{0}$. These amplitudes are comparable to the largest amplitudes of whistler waves measured in the solar wind at 1 AU (Tong et al., 2019a). Namely, the spacecraft measurements showed that whistler waves amplitudes are generally below $0.02\,B_{0}$, while typical amplitudes are of about $3\cdot 10^{-3}\,B_{0}$ (Tong et al., 2019a). In the simulations with smaller density and temperature of the halo population ($n_{c}/n_{e}>0.85$ and $T_{h}/T_{c}<10$), whistler waves would certainly saturate at amplitudes smaller than $0.01-0.05\,B_{0}$ that would be consistent with the typical amplitudes found in the spacecraft measurements. The simulations have shown that the saturated whistler wave amplitude is positively correlated with $\beta_{c}$ that is consistent with the correlation between the whistler wave amplitude and $\beta_{e}$ revealed in the spacecraft measurements (Tong et al., 2019a).

Our results demonstrate that the saturated whistler wave amplitude ${\rm max}(B_{w})/B_{0}$ depends non-monotonously on $q_{e}/q_{0}$. It is related to the fact that the initial (linear) maximum growth rate $\gamma_{\rm max}/\Omega_{e}$ depends on $q_{e}/q_{0}$ non-monotonously (Figure 2), while ${\rm max}(B_{w})/B_{0}\sim A(\gamma_{\rm max}/\Omega_{e})^{\alpha}$. Interestingly, similar correlation between the whistler wave amplitude and $q_{e}/q_{0}$ was revealed in the spacecraft measurements, though the correlation switches from positive to negative at $q_{e}/q_{0}\sim 0.3$ (Tong et al., 2019a). Figure 2 shows that at $T_{h}/T_{c}\sim 4$ that is more typical in the solar wind (e.g., Maksimovic et al., 2005; Tong et al., 2019b), $\gamma_{\rm max}/\Omega_{e}$ reaches maximum at $q_{e}/q_{0}\sim 0.5$. Therefore, in simulations with $T_{h}/T_{c}=4$ the correlation between ${\rm max}(B_{w})/B_{0}$ and $q_{e}/q_{0}$ would switch from positive to negative at $q_{e}/q_{0}\sim 0.5$ that would be in better agreement with the spacecraft measurements. Tong et al. (2019a) have considered whistler waves in the solar wind not disturbed by shocks, but we note that whistler waves produced by the WHFI were also identified around interplanetary shocks (Wilson et al., 2013; Liu et al., 2018) and a positive correlation between the whistler wave amplitude and electron heat flux was reported (Wilson et al., 2013).

The dependence of the saturated whistler wave amplitude on the initial maximum growth rate is similar to the one established by Tao et al. (2017) in the quasi-linear analysis of the classical temperature-anisotropy instability (Sagdeev & Shafranov, 1960; Kennel & Petschek, 1966). That is not surprising, because the nonlinear evolution of the WHFI seems to be quasi-linear in nature. Whistler waves have rather broad spectrum and saturate at relatively small amplitudes (Figure 3), so that the electron scattering is determined by the finite bandwidth, rather than by the finite amplitude of the whistler waves (Figure 4c). Tong et al. (2019a) have recently demonstrated that the quasi-linear theory is indeed applicable to describe effects of whistler waves observed in the solar wind, thereby demonstrating an adequacy of the previous quasi-linear simulations of whistler wave effects in the solar wind (e.g., Vocks et al., 2005; Vocks, 2012).

The presented results indicate that there is a fundamental problem concerning processes controlling the electron heat conduction in the solar wind. We have shown that coherent parallel whistler waves observed in the solar wind (Lacombe et al., 2014; Stansby et al., 2016; Kajdič et al., 2016; Tong et al., 2019a, b) are unlikely to regulate the electron heat flux. We underline that we have addressed the evolution of the WHFI in a uniform plasma. Because the saturation time of the WHFI in physical units is from a few seconds to several minutes (see Figures 5a and 6a in this paper and Figure 16 in Tong et al. (2019a)), whistler waves with typical group velocities of thousands km/s may cover a significant spatial region and can be influenced by low-frequency density and magnetic field fluctuations. The analysis of effects of low-frequency fluctuations on the nonlinear evolution of the WHFI deserves a separate study (see, e.g., Breizman & Ruytov, 1970; Voshchepynets et al., 2015, for analysis of non-uniform plasma density effects on a beam instability evolution). In addition, the evolution of the WHFI may be influenced by the solar wind expansion, which provides the evolution of the background plasma parameters and may continuously support instability of a particular expanding plasma parcel (e.g., Innocenti et al., 2019, for microscopic effects of the expansion). The effects of a slightly oblique propagation (within 20*∘*, see, e.g., Lacombe et al., 2014; Kajdič et al., 2016) of the observed whistler waves should be addressed in the future, though we note that the slight obliqueness may be also due to uncertainty of measurements. Incoherent whistler waves produced due to a turbulence cascade (see, e.g., Gary et al., 2012) could generally regulate the electron heat flux, but spacecraft measurements indicate that the magnetic field turbulence spectrum is likely dominated by kinetic Alfveń waves (e.g., Bale et al., 2005; Salem et al., 2012; Chen et al., 2013) and the contribution of whistler waves to the turbulence spectrum has not been estimated yet (e.g., Narita et al., 2016; Kellogg et al., 2018). Recent PIC simulations showed that rather oblique and large-amplitude electromagnetic whistler waves (wave normal angle of about $45^{\circ}$ and saturated amplitude $B_{w}\sim 0.3B_{0}$) can be generated in a system with a sustained electron temperature gradient and suppress the electron heat flux (Komarov et al., 2018; Roberg-Clark et al., 2018b, a). The analysis of magnetic field spectra showed no evidence for such oblique whistler waves in the solar wind at 1 AU (e.g., Tong et al., 2019a). The likely reason is that the electron VDF developing in those simulations is very different from what is observed in the solar wind. In particular, Landau resonant electrons damp too oblique whistler waves in the solar wind (that is why the WHFI of electron VDFs typical for the solar wind produces only quasi-parallel whistler waves) (Gary et al., 1994), while Landau resonant electrons seem to drive oblique whistler waves in the simulations by Roberg-Clark et al. (2018a). Thus, it is currently not clear what processes regulate the electron heat flux in the solar wind (see Vasko et al., 2019; Verscharen et al., 2019, for a recent discussion), while simulations presented in this paper exclude coherent parallel whistler waves in a uniform plasma from the list of potential processes.

5 conclusion

We have presented results of 1D PIC simulations of the whistler heat flux instability. We summarize the results as follows:

In a uniform plasma, parallel whistler waves produced by the whistler heat flux instability are incapable of regulating the electron heat flux. 2. 2.

Whistler waves saturate at amplitudes consistent with the recent spacecraft measurements. The simulations reproduce the correlations of the saturated whistler wave amplitude with the electron heat flux $q_{e}/q_{0}$ and $\beta_{e}$ revealed in the spacecraft measurements. 3. 3.

The saturated whistler wave amplitudes are dependent on the maximum growth rate computed from the linear theory, ${\rm max}(B_{w})/B_{0}\sim A(\gamma_{\rm max}/\Omega_{e})^{\alpha}$.

The presented simulations exclude coherent parallel whistler waves observed in the solar wind from the list of potential processes regulating the electron heat flux. The effects of non-uniform background plasma density and magnetic field and slightly oblique propagation of the observed whistler waves remain to be addressed.

I. Kuzichev and A. R. Soto-Chavez were supported by the National Science Foundation (NSF) grant No. 1502923. They would like to acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by NSF. The NJIT Kong cluster at New Jersey Institute of Technology (Newark, NJ, USA) was also used. I. Kuzichev would also like to acknowledge the support of the RBSPICE Instrument project by JHU/APL subcontract 937836 to the New Jersey Institute of Technology under NASA Prime contract NAS5-01072. Y. Tong and S. D. Bale were supported in part by NASA contract NNN06AA01C.