The Effect of Solar Wind Turbulence on Parallel and Oblique Firehose Instabilities

We consider the firehose instability coexisting with the omnipresent ambient solar wind turbulence. The characteristic temporal and spatial scales of the turbulence are comparable to those of the instability. Therefore, turbulence may violate the common assumption of a uniform and stationary background used to describe instabilities and make the properties of the instabilities different. To investigate this effect, we perform three-dimensional hybrid simulations with particle-in-cell ions and a quasi-neutralizing electron fluid. We find that the turbulence significantly reduces the growth rates and saturation levels of both instabilities. Comparing the cases with and without turbulence, the former results in a higher temperature anisotropy in the asymptotic marginally stable state at large times. In the former case, the distribution function averaged over the simulation box is also closer to the initial one.


Introduction
The firehose instability is thought to constrain the thermal ion temperature anisotropy in the solar wind when the temperature in the direction parallel to the mean magnetic field becomes larger than the one in the perpendicular direction. At 1 au, this has been shown by Kasper et al. (2002) using the Wind spacecraft data.  further argued that the oblique firehose instability acts more efficiently in the slow solar wind, while the parallel one is more important in the fast wind. Closer to the Sun, Huang et al. (2020) have found with the help of the Parker Solar Probe data that the proton temperature anisotropy is consistent with the limits imposed by the firehose instability. Despite their common name, the parallel and oblique firehose instabilities have distinctly different properties . In particular, the former is resonant and the latter is nonresonant, and the oblique modes may be growing even when the parallel ones are marginally stable. When both instabilities are present, the protons are isotropized mainly by the oblique mode.
In the solar wind, any instability has to operate in the ambient turbulent environment and the relationship between them is not clear. The turbulence is nearly omnipresent as it transports the energy of very-large-scale fluctuations generated by the Sun to much smaller kinetic scales (Coleman 1968;Jokipii 1973;Matthaeus & Goldstein 1982;Tu et al. 1984;Matthaeus et al. 1994;Goldstein et al. 1995;Tu & Marsch 1995;Zank et al. 1996). The description of instabilities is usually based on the assumption of a uniform and stationary background. However, this assumption is violated by the ambient turbulence. The nonuniform medium is an essential factor that needs to be taken into account (Hellinger et al. 2017;Ozak et al. 2015;Ofman et al. 2017;Markovskii et al. 2019a). In particular, an inhomogeneity across the background magnetic field makes the unstable waves more oblique. This can change the character of the wave-particle interaction and the resulting particle heating (Markovskii et al. 2019a).
The turbulent cascade also makes the medium in which the instabilities grow time dependent (Markovskii et al. 2006). The characteristic time, τ, can be estimated in the solar wind as the nonlinear eddy turnover time of a strong turbulence. Based on the rms amplitude of the magnetic field fluctuations at the proton kinetic scales, τ is usually around 30 p 1 Wat 1 au, where Ω p is the proton gyrofrequency. If τ is comparable to or shorter than the inverse growth rate of an instability, γ −1 , then the variation of the background over the time of the instability growth has to be taken into account (Mikhailenko & Stepanov 1984;Maslennikov et al. 1995;Markovskii & Hollweg 2002).
The temperature anisotropy instability coexisting with the turbulence has been analyzed by Markovskii et al. (2019bMarkovskii et al. ( , 2020 with the opposite sense of the anisotropy, where the perpendicular temperature is larger than the parallel one. To some extent, this setup is similar to the case considered here because it also produces two separate competing instabilities, a parallel proton-cyclotron and an oblique mirror mode. However, differences arise even before the turbulence is taken into account. In particular, when the particle distribution functions are bi-Maxwellian, the presence of alpha particles and other minor ions reduces the linear growth rate of the proton-cyclotron but not the mirror mode (Price et al. 1986;McKean et al. 1994). At the same time, both parallel and oblique firehose instabilities are affected by the minor ions . Furthermore, the protoncyclotron mode is more efficient at isotropizing the protons than the mirror mode, whereas the protons are mostly isotropized by the oblique, but not the parallel, firehose fluctuations.
Adding a turbulent ambient medium results in a considerable modification of the proton-cyclotron and mirror instabilities. The initial growth rates of the mirror mode are close, but the saturation level is significantly reduced when the turbulence is present. The saturation level of the proton-cyclotron mode is not affected as strongly. However, the turbulence extends the spectrum of these unstable waves to higher perpendicular wavenumbers making them more oblique, in a way similar to what would occur in the presence of a static nonuniform background. In the present paper, we investigate the effect of an approximately two-dimensional (2D) turbulence, with wavevectors highly oblique to the mean magnetic field, on the parallel and oblique firehose instabilities.

Numerical Setup
We perform three-dimensional (3D) hybrid simulations with particle-in-cell protons and a quasi-neutralizing electron fluid. The numerical code is described by Terasawa et al. (1986) and Vasquez (1995Vasquez ( , 2015. The code solves the following equations: where n e = n p + 2n α . The quantities x p,α and v p,α are positions and velocities of individual protons and alpha particles, q p,α is the ion charge, E is the electric field, c is the speed of light, e is the elementary charge, n e is the electron number density, and V e is the electron fluid velocity. The proton and alpha number density n p,α and bulk velocity V p,α for each spatial cell are calculated as moments of the distribution. The mean magnetic field B 0 is in the positive x direction. Equations (1)-(5) neglect the global evolution of the solar wind because all the characteristic spatial and temporal scales in our study are much smaller than the global-expansion scales. Furthermore, the instabilities that we are interested in are often seen on the fringes of the observed parameter distribution. Therefore, they are not necessarily caused by the quasi-static evolution associated with the expanding wind, but rather by some transient disturbances. To simplify the numerical model, the electron temperature is set to zero. Equations (1)-(5) also imply that the electron fluid is massless. The simulation grid is 256 cells in the x direction and 128 cells in the y and z directions. The simulation box size is L ∥ = 256 in the x direction and L ⊥ = 64 in the y and z directions in units of the proton inertial length V .
A p 1 W -The time step is 0.01 in units of the inverse proton gyrofrequency .
p 1 W -Here, Ω p and the Alfvén speed, V A , and thereby the spatial and temporal scales, are defined with the initial mean values of the magnetic field, B 0 , and the proton number density, n 0p . The boundary conditions are periodic.
We initiate a turbulent cascade as follows. The magnetic field fluctuations ΔB(t, x) at t = 0 are given by the formula where x is the Cartesian spatial position vector, k is the wavevector, and f(k) is a random phase. The seed spectrum is confined to the modes (k x = 0, k y = ±2π/L ⊥ , k z = ±2π/L ⊥ ) with all modes having the same amplitude. The proton and alpha bulk velocity fluctuations ΔV p,α are defined in the same way as ΔB in Equation (6). The components of the vectors δB (0, k) and δV p,α (0, k) obey the polarization relations of linear Alfvén waves in the MHD limit. This is merely a convenient way to control the initial cross-helicity and does not mean that such waves exist in our configuration. In the present case, the cross-helicity is zero. The relative total rms amplitude of the magnetic field fluctuations is set to The ion distribution functions are loaded as a spatially uniform drifting bi-Maxwellian: is the total ion temperature, and the subscripts "∥" and "⊥" refer to the direction of the background magnetic field. The initial total proton beta β p = 1.5. The proton temperature anisotropy T ⊥p /T ∥p = 0.38. The temperature of the alpha particles T α = 2T p is isotropic and their number density n 0α = 0.05n 0p . The initial ion number densities and temperatures are uniform. The relative drift speed between the ion species averaged over the simulation box is zero. The number of proton and alpha particles per cell is 1500 and 500, respectively. The system is then allowed to evolve freely.

Simulation Results
The ambient turbulent medium used in this paper has been described by Markovskii et al. (2019bMarkovskii et al. ( , 2020. Strictly speaking, a decaying turbulence is not in a steady state. However, it goes into a quasi-steady phase after a rapid initial relaxation. The turbulence is then consistent with the observed properties of the solar wind fluctuations at 1 au, in particular the rms amplitude at the ion kinetic scale. It obeys the Kolmogorov law and is weakly compressive at larger scales. At smaller scales, the spectrum steepens and the fluctuations become more compressive. The average temperature anisotropy 〈T ⊥p 〉/〈T ∥p 〉 resulting from the coexisting turbulence and instability is displayed by a black line in Figure 1. Here the angle brackets denote averaging over the simulation box. The relaxation phase manifests itself as oscillations at early times ( t 60 p W  ). The turbulence develops and enters the quasi-steady phase around t 180 .
p 1 = W -To separate the contributions from the turbulence and instability, we ran a 2D simulation of the turbulence without the instability. The initial parameters were the same as in the 3D setup, and the simulation plane was perpendicular to the background magnetic field. Since the parallel wavenumbers were cut off, the instability was suppressed. The temperature anisotropy in this case is plotted in Figure 1 as a blue line. At first, it follows closely the black curve describing the 3D configuration but then the latter deviates far away.
As the system evolves, the distribution function averaged over the simulation box becomes non-bi-Maxwellian. However, more importantly, the turbulence makes the distribution function inhomogeneous and time dependent. This means that the turbulence can affect the instability not just by interacting with the unstable waves: it also changes the environment where the instability develops before it is even excited. As a result, a direct comparison of the instabilities with and without turbulence becomes impossible even if the averaged distributions are the same. This is true regardless of the structure of the distribution in the phase space. Therefore, the instability in a uniform and stationary background can be used only for reference.
To have a point of reference, we ran a simulation of the instability without turbulence. This is achieved by putting A = 0 in Equation (7). The initial proton distribution function was bi-Maxwellian with T ⊥p /T ∥p = 0.4 and β p = 1.5. For a better match, the temperature anisotropy was set to the time average of the rapid oscillations in the turbulent case. The same total beta as in the turbulent case was chosen because β p does not undergo rapid oscillations and remains almost constant throughout the relaxation phase. It is also instructive to evaluate the contribution of the alpha particles to the instability and compare it to the effect of the turbulence. Therefore, we performed a separate simulation run with n 0α = 0.
The time evolution of the average temperature anisotropy is shown by a green line in Figure 1 for A = 0 and n 0α = 0.05n 0p and by a magenta line for A = 0 and n 0α = 0. As can be seen from the figure, the alpha-particle concentration does not have a strong effect on the temperature anisotropy. At the same time, the turbulence produces a significant difference. It is also clear that the initial bump during the relaxation phase is due to the turbulence because there is very little contribution from the instability at those times.
The intensity of the magnetic field fluctuations associated with the instability is displayed in Figure 2 as a function of time. As will be demonstrated below, the firehose modes are    well separated from the turbulence in the wavenumber space.
To eliminate the turbulence and single out the instability, we sum the intensity of the modes over a certain fraction of the wavenumbers defined below. We notice that the alpha-particle concentration does affect the wave intensity somewhat, unlike the anisotropy. It suppresses the oblique modes and boosts the parallel ones. Nevertheless, the particle isotropization still occurs during the oblique-mode phase, while the parallel mode plays a minor role. Compared to the alpha-particle concentration, the effect of the turbulence is more significant. It reduces both growth rate and the saturation level of the two instabilities. Note that the instability in the turbulent case does not diminish in time. Rather, it becomes weaker from the very beginning. The reason, as pointed out above, is the distribution function becomes spatially nonuniform. A 2D spectrum, P, of the magnetic field fluctuations for the instability coexisting with turbulence is plotted at the peak intensity of the parallel and oblique firehose modes in Figures 3  and 4, respectively. The spectrum has been summed over the azimuthal angle and normalized to B . 0 2 The separation of the instability from the turbulence in the wavenumber space is visible in the figures, although some overlap is present for the oblique modes. To calculate the curves in Figure 2, we summed the modes in two regions with 0.15 < |k ∥ |V A /Ω p < 0.27 and |k ⊥ |V A /Ω p < 0.39, and 0.27 < |k ∥ |V A /Ω p < 0.98 and |k ⊥ |V A /Ω p < 1.08.
The spatial inhomogeneity associated with the turbulent fluctuations extends the spectrum of the unstable waves to higher perpendicular wavenumbers making them more oblique (Markovskii et al. 2019b). However, for the oblique modes, this effect is not as strong as for the parallel modes. To evaluate the spectral anisotropy from the simulation results, we use a   A p 1 W -The turbulence depends predominantly on the y and z coordinates and produces considerable deviations of the anisotropy from the average value (seen in the x = 0 plane). At t = 200, only parallel-propagating unstable modes are present and their intensity is substantially lower than that of the turbulence (as follows from Figure 3). Furthermore, even though the rapid oscillations of the average anisotropy stop in the quasi-steady phase, local time variations continue. These variations are associated with the turbulent eddies. The turnover time, τ, of the largest eddy in the system can be estimated as ( ) L A 2 9 3 .
that the Kolmogorov scaling of the eddy amplitude continues to the kinetic scales, τ at the proton inertial length becomes 20 .
p 1 W -As follows from the discussion in Section 1, this is smaller than the typical value observed by Markovskii et al. (2006) in the solar wind but not by much.
At later times, the local variations of the anisotropy associated with the decaying turbulence subside and it is useful to examine the asymptotic particle distribution at the end of the simulation run (t 800 p 1 = W -). The proton distribution summed over v z is displayed in arbitrary units in Figure 8 as a function of v x and v y in the turbulent case. The velocity is in units of the Alfvén speed, V A . The distortion of the initial bi-Maxwellian distribution shows up as a flat top at |v x | ≈ 0 and "horns" at |v x | ≈ 1.25. Without the turbulence, the distortion is qualitatively the same but noticeably stronger (Figure 9), consistent with the stronger isotropization (Figure 1). To provide a measure of how much the distribution changes over time, we plot an initial uniform bi-Maxwellian distribution in Figure 10.
The asymptotic distribution and anisotropy are closer to their unstable counterparts in the turbulent case than without turbulence, and yet the instability is saturated. One possible reason is that the ambient medium is still not uniform and stationary at the end of the simulation run, so conventional instability criteria do not apply. The spatial 3D structure of the anisotropy at this time is displayed in Figure 11. We see that the inhomogeneity is not as strong as it was earlier (in Figure 7), but is present nonetheless. Without the turbulence, the corresponding structure is much more uniform and appears as a monochrome green color (not shown). Another explanation of the different saturation efficiency is the mechanism by which it happens. As shown by Hellinger & Matsumoto   (2001), the oblique firehose instability saturates by nonlinear interactions and transformation of the unstable waves. The turbulent fluctuations are likely to contribute to these interactions and, thereby, alter the saturation level.

Conclusion
We have carried out 3D hybrid simulations of the parallel and oblique firehose instabilities in a turbulent ambient plasma with alpha particles. The instability was excited by the protons with temperature perpendicular to the mean magnetic field smaller than the parallel temperature. It has a complex structure, with waves interacting and being reabsorbed by the particles. The turbulence was initiated with wavevectors perpendicular to the mean magnetic field and remained approximately 2D during the simulations. Initially, the turbulence spectrum spreads rapidly from lower to higher wavenumbers. Later, it enters a quasi-steady phase where the variation is much slower but the turbulence is still decaying.
We used the parameters of the turbulence close to those observed in the solar wind. The properties of the instability were compared to the case of a uniform and stationary background. In this case, we also evaluated the contribution of the typical alpha-particle concentration to the instability. We found that the proton isotropization is caused mostly by the oblique firehose modes regardless of the turbulence. However, the turbulence significantly reduces the growth rates and saturation levels of the parallel and oblique firehose instabilities and its effect is stronger than that of the alpha particles.
The asymptotic proton distribution at large times gets distorted from its initial bi-Maxwellian structure. The distortion is weaker when the turbulence is present than when it is not, consistent with the smaller change of the anisotropy. Therefore, the instability is saturated more efficiently. In part, this may be due to the fact that the distribution is still spatially inhomogeneous and time dependent even at the end of our simulation run, despite the decaying turbulence. Also, the saturation of the oblique firehose instability occurs through nonlinear interactions of the unstable modes with each other. The turbulent fluctuations are likely to modify this process by interacting with the oblique unstable modes, especially because the two are close in the wavenumber space.