Particle-In-Cell simulations of the parallel proton firehose instability influenced by the electron temperature anisotropy in solar wind conditions

TL;DR

This study uses Particle-In-Cell simulations to investigate how electron temperature anisotropy influences the development of the parallel proton firehose instability in solar wind conditions, revealing that anisotropy accelerates instability onset and promotes isotropization.

Contribution

It provides the first non-linear kinetic simulation analysis of electron influence on the proton firehose instability in solar wind plasma.

Findings

Electron anisotropy accelerates firehose instability onset.

Wave fluctuations reduce particle temperature anisotropy.

Simulation results agree with linear theory.

Abstract

In situ observations of the solar wind show a limited level of particle temperature anisotropy with respect to the interplanetary magnetic field direction. Kinetic electromagnetic instabilities are efficient to prevent the excessive growth of the anisotropy of particle velocity distribution functions. Among them, the firehose instabilities are often considered to prevent the increase of the parallel temperature and hence to shape the velocity distribution functions of electrons and protons in the solar wind. We present a non-linear modeling of the parallel firehose instability, retaining a kinetic description for both the electrons and protons. One-dimensional (1D) fully kinetic Particle-In-Cell simulations using the Energy Conserving semi-implicit method (ECsim) are performed to clarify the role of the electron temperature anisotropy in the development of the parallel proton firehose instability. We found that in the presence of an electron temperature anisotropy, such that the temperature parallel to the background magnetic field is higher than the temperature in the perpendicular direction, the onset of the parallel proton firehose instability occurs earlier and its growth rate is faster. The enhanced wave fluctuations contribute to the particle scattering reducing the temperature anisotropy to a stable, nearly isotropic state. The simulation results compare well with linear theory. A test case of 1D simulations at oblique angles with respect to the magnetic field is also considered, as a first step to study the cumulative effect of protons and electrons on the full spectrum of instabilities.

Figures

Tables

Table 1. Table 1: Input parameters

Case	${\beta }_{p,\parallel }$	${\beta }_{e,\parallel }$	$A_p$	$A_e$
1 (PFHI)	4.0	4.0	0.1	1.0
2 (PFHI)	4.0	0.1	0.1	0.1
3 (EFHI)	4.0	4.0	1.0	0.1
4 (PFHI + EFHI)	4.0	4.0	0.1	0.1

Equations

$$F_j(v_{\parallel },v_{\perp })=$$

$${\tilde{k}}^2=$$

$$+A_p-1+\frac{A_p\left(\tilde{\omega }+1\right)-1}{\tilde{k}\sqrt{{\beta }_{p,\parallel }}}{Z}_p\left(\frac{\tilde{\omega }+1}{\tilde{k}\sqrt{{\beta }_{p,\parallel }}}\right),$$

$$Z_j\left({\xi }_j^{\pm }\right)=\frac{1}{{\pi }^{1/2}}{\int }_{-\infty }^{\infty }\frac{\exp \left(-x^2\right)}{x-{\xi }_j^{\pm }}dt$$

Full text

plasmas — solar wind — instabilities — waves

1 Introduction

Solar wind is a natural plasma laboratory, enabling in situ measurements and investigations of a rich variety of physical processes operating there. More than 50 years of observations (made e.g. by Helios, Ulysses, and Wind space missions) revealed that proton and electron velocity distribution functions (VDFs) may be markedly non-Maxwellian, and therefore non-thermal. The observed distributions are often anisotropic with respect to the magnetic field direction (see Marsch et al., 1982; Marsch, 2006; Štverák et al., 2008; Kasper et al., 2016, and references therein). Deviations from isotropy are limited by such mechanisms as particle-particle (binary) collisions and wave-particle interactions, resulting in fluctuations of electric and magnetic fields which scatter plasma particles and reduce their temperature anisotropy. The fluctuations may be either pre-existing (e.g. turbulence) or generated by plasma instabilities. These mechanisms are usually invoked in the models to explain the observed VDFs (Hellinger & Trávníček, 2006; Štverák et al., 2008; Bale et al., 2009; Chen et al., 2016; Shaaban et al., 2017; Lazar et al., 2017; Shaaban et al., 2019a, b; Ofman, 2019).

Solar wind is a hot and low-density plasma, such that binary collisions are not efficient to limit the temperature anisotropies of plasma particles. Despite significant efforts in fluid modeling (see e.g. the review by Hunana et al., 2019), a kinetic approach is needed to describe other mechanisms potentially constraining the anisotropies, e.g., the plasma kinetic instabilities and wave-particle interactions, and to explain the low level of anisotropy recorded in situ for both electrons and protons at 1 AU (Matteini et al., 2007; Štverák et al., 2008). A kinetic approach is also important to take into account wave-like processes operating at frequencies ranging between the ion gyrofrequency and the electron plasma frequency (Bruno & Carbone, 2013). Spectra of wave-like fluctuations of the solar wind electric and magnetic fields are well established by the observations, including those at kinetic scales. These waves can resonantly interact with particles resulting in the exchange of energy and momentum between them (Alexandrova et al., 2013). These are key aspects to understand the kinetic processes in the solar corona and the solar wind.

Kinetic simulations enable detailed studies of the plasma instabilities and their consequences, e.g., enhancement of the observed wave fluctuations and their reaction on particles that may lead to their energization. They allow for modelling the plasma from first principles and can correctly describe non-thermal behaviours, anisotropic distribution functions and non-symmetric pressure tensors in the solar wind plasma. Among the different numerical techniques, the Particle-In-Cell (PIC) method (Birdsall & Langdon, 2004; Hockney & Eastwood, 1988) represents our choice because the algorithm is scalable and can take advantage of modern high-performance computing systems (Gonzalez-Herrero et al., 2018).

Firehose instabilities are thought to be a crucial mechanism for constraining the electron and proton temperature anisotropies observed in expanding astrophysical plasmas, such as the solar wind (Bale et al., 2009; Chen et al., 2016; Lazar et al., 2017; Shaaban et al., 2019a). These electromagnetic kinetic instabilities are thought to develop around 1 AU and beyond (Matteini et al., 2006), but perhaps also closer to the Sun (van der Holst et al., 2019). Firehose instabilities can develop in a plasma with sufficiently large parallel plasma beta $\beta_{j,\parallel}=8\pi n_{j}k_{B}T_{j,\parallel}/B_{0}^{2}>~{}1$ and anisotropy $A_{j}=T_{j,\perp}/T_{j,\parallel}<1$, where $j$ refers to a given plasma species ($j=e,p$ for electrons and protons), $n$ and $T$ are the particle number density and temperature respectively, $k_{B}$ the Boltzmann constant, while $\perp$ and $\parallel$ denote directions perpendicular and parallel to the background magnetic field, $\bm{B}_{0}$.

A number of studies have been conducted to demonstrate the co-existence of the firehose instabilities with evolving turbulent spectra (Hellinger et al., 2015), to confirm the evolution of Alfvénic fluctuations in the firehose regime (Tenerani & Velli, 2018) and to test the properties of the instabilities in the expanding solar wind (Hellinger et al., 2003; Matteini et al., 2006; Innocenti et al., 2019a, b; Hellinger et al., 2019). The proton firehose instability has been largely examined via hybrid simulations, e.g. (Matteini et al., 2006; Hellinger, 2017), where the interest is focused only on the ion-scale processes, since the electrons are considered as a massless, charge neutralizing fluid. The proton firehose instability evolves out of the whistler wave and has right-hand circular polarization (Gary, 1993).

Two branches of the firehose instability can be driven by the anisotropic proton distributions: the periodic proton firehose with the finite real part of the wave frequency $\omega_{r}\neq 0$ and the most unstable modes in the direction parallel to the magnetic field, and the aperiodic modes ($\omega_{r}=0$) that are present only for oblique angles of propagation, (e.g. Yoon et al., 1993; Hellinger & Matsumoto, 2000; Hunana & Zank, 2017). The properties of these two modes, when predicted using simplified models neglecting the effects of electrons, do not always agree with observations and contradictory conclusions were sometimes obtained. For example, Kasper et al. (2002) compared observations of proton velocity distribution functions with the limits for parallel and oblique firehose instabilities derived from the linear theory and simulations as reported by Gary et al. (1998). The agreement with the parallel mode was better. Hellinger & Matsumoto (2000) found that the oblique proton firehose instability may have the growth rate comparable or even grater than that of the parallel proton firehose. Hellinger & Matsumoto (2001) extended the work by Hellinger & Matsumoto (2000) and determined that the parallel firehose instability initially dominates and saturates faster than the oblique one, which then can grow even when the parallel modes saturate. Hellinger et al. (2006) and Bale et al. (2009) state that the periodic (parallel) mode should develop faster, but the observed limits of temperature anisotropy are better shaped by the thresholds of the aperiodic (oblique) firehose. However, Hellinger & Trávníček (2006) analyzed proton firehose instabilities in the presence of alpha particles and concluded that the parallel and oblique modes have comparable growth rates for a vast range of plasma parameters (see also Hunana & Zank, 2017; Hunana et al., 2019) and that both modes are relevant in the solar wind context.

On the other hand, the observations confirm the (co-)existence of solar wind electrons with $T_{e\parallel}>T_{e\perp}$ (Marsch, 2006; Štverák et al., 2008), which can trigger two similar branches of the electron firehose instability. There is no unique definition for the two branches of the electron firehose instability. Li & Habbal (2000) define the quasi-parallel branch as periodic modes in directions within 30∘ with respect to the magnetic field, while the oblique branch represents aperiodic modes at larger angles. Gary & Nishimura (2003) define the parallel electron firehose as non-resonant while the oblique mode is resonant. Camporeale & Burgess (2008) distinguish propagating ($\omega_{r}\neq 0$) and non-propagating ($\omega_{r}=0$) modes. Even though Camporeale & Burgess (2008), in agreement with Li & Habbal (2000), did not find any non-propagating perturbation at small angles, they demonstrated the extension of propagating perturbations to relatively high oblique angles. Both branches have been investigated but generally this has been done neglecting their interplay with the proton firehose instability (Quest & Shapiro, 1996; Paesold & Benz, 1999; Li & Habbal, 2000; Gary & Nishimura, 2003; Paesold & Benz, 2003; Camporeale & Burgess, 2008; Hellinger et al., 2014; Shaaban et al., 2019c; López et al., 2019). Neglecting the interplay of different species may prevent understanding of the instability implications and the anisotropy relaxation mechanisms (Yoon & Sarfraz, 2017; Yoon et al., 2019).

Recently, more realistic approaches have been proposed showing that linear properties of both branches can be markedly altered by the interplay of anisotropic electrons and protons (Michno et al., 2014; Maneva et al., 2016; Shaaban et al., 2017) and may thus lead to a better agreement with observations. However, these approaches do not allow analyzing such phenomena as the effect of the wave growth on the scattering of particles and hence on the rate of their isotropization. In an attempt to check and complete these predictions from linear theory, in this work we perform one-dimensional (1D) PIC simulations to investigate the periodic parallel branch of the proton firehose instability under the effect of anisotropic electrons. We are aware that reducing the scope of our investigation to perturbations propagating only in the parallel direction may be a limitation. The oblique electron firehose instability has generally a lower threshold and higher growth rate than the parallel electron firehose instability. Nevertheless, quasi-parallel and exactly parallel modes may be dominant after the saturation of the linear stage (Camporeale & Burgess, 2008; Innocenti et al., 2019a). Even if our work does not deal with the competition between parallel and oblique modes, it is the first numerical study that looks at the nonlinear evolution of the parallel proton firehose instability in the presence of anisotropic electrons using fully kinetic simulations, which retain a description from first principles for both electrons and ions.

This paper is organized as follows. Section 2 illustrates the employed simulation setup. Section 3 reports the main results of our simulations. In Section 4 simulation results are compared with the linear theory. Section 5 presents the results of a few tests of simulations at oblique angles with respect to the magnetic field. Finally, Section 6 presents a discussion of the simulation results in the solar wind context, and reports the conclusions.

2 Setup of the PIC simulations

In order to study the development and evolution of the parallel firehose instability in the solar wind, we performed 1D kinetic simulations using the ECsim code (Energy Conserving semi-implicit method, see Lapenta, 2017; Lapenta et al., 2017; Gonzalez-Herrero et al., 2018, 2019). ECsim is a semi-implicit PIC code which, unlike hybrid methods, is capable to retain kinetic electron and ion information, such as wave-particle interaction and non-Maxwellian VDFs. It has been demonstrated that the code is stable and energy conserving over a wide range of spatial and temporal resolutions (Lapenta et al., 2017; Gonzalez-Herrero et al., 2018). In particular, the latter aspect is crucial for the current work, because it allows us to resolve only the scales of interest. This means that we do not need to resolve the electron Debye length as in traditional PIC codes (Cohen et al., 1989; Brackbill & Forslund, 1982; Lapenta, 2017). As a consequence, we can follow the dynamics of both electrons and protons from a kinetic point of view for very long times, thus exploring the interplay of electrons and protons, and the long-term evolution of the firehose instability.

We model a collisionless plasma composed of electrons and protons with real mass ratio $\mu=~{}m_{p}/m_{e}=~{}1836$, immersed in a uniform background magnetic field $\bm{B}_{0}$. The magnetic field is along the $x$-direction ($\bm{B}_{0}=B_{0}\hat{e}_{x}$), and its magnitude $B_{0}=2.5\ 10^{-4}$ G is such that $\omega_{e}/\Omega_{e}=63.24$ and $\omega_{p}/\Omega_{p}=2709.98$, with $\omega_{j}=\sqrt{4\pi e^{2}n_{j}/m_{j}}$ and $\Omega_{j}=\sqrt{eB_{0}/c\,m_{j}}$ being the plasma and the cyclotron frequencies (for species $j$), respectively, $e$ is the elementary charge, $m$ is the species mass and $c$ is the speed of the light in vacuum. The initial values of the electron number density ($n=25$ cm*-3*), magnetic field magnitude and $\omega_{e}/\Omega_{e}$ are typical for the solar wind around 0.3 – 0.5 AU, (e.g. Venzmer & Bothmer, 2018; Tong et al., 2019). A simulation box with the length of $60\,d_{i}$ has been employed, where $d_{i}=c/\omega_{p}$ is the ion inertial length. The size of the box is chosen so that more than $20$ wavelengths of the most unstable mode fit into the box. A cell size $\Delta x\simeq 0.074\,d_{i}$ and a temporal step $\Delta t=0.5\,\omega_{p}^{-1}$ have been chosen. We therefore fully resolve the ion inertial length, and characteristic electron scales are resolved only marginally. This is sufficient as the electron firehose instability occurs at proton scales (Gary et al., 1998; Li & Habbal, 2000). Both species are initially described by bi-Maxwellian velocity distribution functions with no drift velocity:

[TABLE]

with $u_{j}=\sqrt{2k_{B}T_{j}/m_{j}}$ is the thermal velocity of the species $j$.

In all the simulations we use $10^{4}$ particles per cell per species, and adopt periodic boundary conditions for fields and particles. Simulations performed with different resolution and number of particles per cell yield similar results, showing a good convergence of the code. This convergence study together with the possibility of easily comparing the numerical outcomes with the linear theory of the periodic proton firehose instability, led us to infer that the code is correct and accurate, even when electron scales are resolved only marginally.

3 Results of the PIC simulations

In order to determine the role of the interplay of electrons and protons in the firehose instability, we compare four cases: the pure proton firehose instability (PFHI) where only the protons are anisotropic (case 1), the PFHI where electrons are also anisotropic but characterized by low plasma beta (case 2), the pure electron firehose instability (EFHI) where only the electrons are anisotropic (case 3), and the case characterized by the temperature anisotropy of both species (PFHI + EFHI) and high electron plasma beta (case 4), see Table 1. The values of initial anisotropies $A_{j}=T_{j,\perp}/T_{j,\parallel}$ are perhaps lower with respect to the average solar wind conditions. This choice is mainly motivated by simulation reasons, where the computational resources constitute a serious constraint. The simulations with higher values of $A_{j}$ are characterized by lower growth rates and hence they run significantly longer before achieving the saturation. All the other parameters such as electron and proton number densities, magnitude of the magnetic field, ratio between the electron plasma frequency and electron cyclotron frequency are typical for the solar wind as mentioned in Section 2. Finally we note that the values of anisotropy and plasma beta in the solar wind are affected by large fluctuations with values that approach those in our simulations (e.g. Kasper et al., 2003; Hellinger et al., 2006; Bale et al., 2009).

In Figure 1 the evolution of the proton velocity distribution function is displayed for the PFHI case 1. Proton VDFs are shown at the initial time step (Figure 1(a)), during the development of the instability (Figure 1(b)-(e)) and after the instability saturation (Figure 1(f)). One can notice that protons are scattered towards higher perpendicular velocities by resonant wave-particle interactions. The process continues until a condition of marginal stability is reached. The isotropic species, i.e. the electrons, does not play any particular role in the development of the instability. The electron temperature is nearly constant during the simulation and their VDF remains Maxwellian.

Case 2 of Table 1 is characterized by both species anisotropic but with electron plasma beta too low to excite the electron firehose instability. Similarly to case 1, only the proton firehose instability develops, and the evolution of the proton VDF is very similar to that shown in Figure 1 so it is not shown here.

Figure 2 shows the evolution of the electron VDF in the EFHI case 3. The anisotropy triggers the instability and the free energy of the particles is converted into magnetic energy. The anisotropy is reduced during the evolution of the instability through the scattering between particles and electromagnetic waves (Figure 2(b)-(e)) until a complete isotropization (Figure 2(f)). Similarly to case 1, the isotropic species, the protons, does not participate to the development of the instability. If one compares Figure 1 and Figure 2, i.e. the evolution of the PFHI with the evolution of the EFHI, it can be clearly seen that the latter develops on shorter time scales than the former, so the process of wave-particle interactions is more efficient and rapid in isotropizing the anisotropic species. Electrons get heated in the perpendicular direction mainly through non-resonant wave-particle interaction, as it should be for the parallel EFHI, which is usually not resonant. As mentioned by Messmer (2002), the non-resonant wave-particle interaction results in electrons with no preferential $v_{\parallel}$ being pitch-angle scattered. This results in an increase of $v_{\bot}$. However, the electron VDF at $t=18000\,\omega_{p}^{-1}$ has the characteristic butterfly shape and shows a slight signature of a resonance. This is certainly due to the high level of anisotropy. The EFHI becomes indeed resonant with the electrons for high anisotropy values (Paesold & Benz, 1999).

Figure 3 reports the evolution of the electron and proton VDFs in the PFHI + EFHI case 4. The particle VDFs at $t=0$ (Figure 3(a)-(b)), during the instability (Figure 3(c)-(f), (h) and (j)), and at the instability saturation (Figure 3(g), (i), (k) and (l)) are shown. Electrons start to isotropize much earlier than protons (Figure 3(c)) and absorb a part of the free energy before wave-proton interaction takes place. Once electrons reach the condition of marginal stability, protons start to be scattered by the interaction with electromagnetic waves until their isotropization. Comparing Figure 1 with Figure 3, one can note that in the PFHI + EFHI simulation, where both species are initially anisotropic, the exchange of energy between protons and the waves starts earlier than in the PFHI case where only protons are anisotropic.

To check in detail how particles interact with the generated electromagnetic waves in the PFHI + EFHI case 4, we present the VDFs for electrons and protons integrated over the parallel and perpendicular directions (Figure 4). In Figure 4(a)-(b) the initial anisotropic distributions are shown. At $t=28000\,\omega_{p}^{-1}$ the perpendicular electron distribution function remains Maxwellian whereas the parallel VDF becomes a flat-top distribution function (Figure 4(c)). Via the process of scattering on the electromagnetic waves generated due to the instability, the proton distribution function gains suprathermal tails in the direction parallel to the background magnetic field (Figure 4(d)). This represents a signature of a preferential resonant interaction. Finally, while the electrons at the end of the simulation show perfect isotropy with the distribution functions along the parallel and perpendicular directions being essentially identical (Figure 4(e)), the protons are still slightly anisotropic, with the thermal spread in the parallel direction still larger than the thermal spread in the perpendicular direction (Figure 4(f)).

Figure 5 shows the evolution of proton and electron temperature anisotropies, $A_{p}$ and $A_{e}$, respectively for cases 1, 2 and 4 reported in Table 1. Both species, when they are unstable with respect to the firehose instability, are subjected to parallel cooling and perpendicular heating. The result is a reduction of the initial anisotropy with time and hence the corresponding increase of $A_{p}$ and $A_{e}$. Evolution of the proton temperature anisotropy for the PFHI is essentially the same for cases 1 and 2. This is due to the fact that the electrons, in both the cases, are stable with respect to the parallel electron firehose. The parallel proton firehose instability affects only marginally the electrons when they are stable with respect to the electron firehose instability. Indeed, in these two cases the electron temperature anisotropy $A_{e}$ presents only small oscillations around the stable initial condition. This means that the proton firehose instability can be properly described also via hybrid simulations if electrons are not subject to the electron firehose instability.

For the PFHI + EFHI (case 4), characterized by the temperature anisotropy of both species and high electron plasma beta, one can notice that the electrons reach the isotropic state on a shorter time scale than protons. In this case, the interplay between electron and proton anisotropies modifies the onset of the proton firehose instability, which leads to an earlier isotropization of the proton distribution functions, if compared with cases 1 and 2. Moreover, PFHI + EFHI case 4 exhibits the final proton temperature anisotropy that is slightly lower compared to that obtained for the two PFHI cases 1 and 2, where only the protons are unstable with respect to the firehose instability. This result is in agreement with the works by Yoon & Sarfraz (2017) and Yoon et al. (2019), where the coupling between electrons and protons is considered necessary to explain the broad distribution of particles near isotropic conditions in the solar wind. In addition, we note that the final value of $A_{p}$ is around $0.7$, which is in agreement with the most probable level of proton anisotropy observed in the solar wind at 1 AU (Kasper et al., 2003).

Figure 6 shows the evolution of the out-of-the-plane magnetic field $B_{z}$ for the three analyzed cases. It demonstrates the development of the wave fluctuations and the increase of their amplitude with time as particle isotropize due to the firehose instability. The more rapid development of the instability in the PFHI + EFHI case with respect to the PFHI case can be easily noticed. Comparing the three panels of Figure 6, it is possible to identify characteristic time scales for the rise of the magnetic fluctuations. The instability time scale depends on the presence or absence of the interplay between proton and electron anisotropies. In case of anisotropic electrons and isotropic protons (EFHI case) the instability develops relatively fast, and the magnetic fluctuations reach the peak amplitude around $20000$ $\omega_{p}^{-1}$. When bi-Maxwellian distributions are initialized for both protons and electrons, the resulting kinetic instability takes more time to develop with respect to the EFHI case, but it appears almost twice faster than in the PFHI case. In all cases it is clear from Figure 6 that the parallel firehose instability is characterized by propagating modes ($\omega_{r}\neq 0$).

4 Comparison of simulation results with the linear theory

In this section we perform a comparison between the linear theory and our PIC simulations.

The linear stability analysis is based on the linearized Vlasov-Maxwell dispersion relation for the wave propagation parallel to the background magnetic field in a homogeneous and bi-Maxwellian electron-proton plasma. Linear theory can predict the fastest growing plasma modes and describe their dispersive characteristics such as growth rates, wave frequency and instability thresholds. However, linear theory cannot describe neither the interaction between the different plasma modes nor the energy transfer between plasma species. Only PIC and hybrid simulations contain the full description of the nonlinear physics.

For the firehose unstable modes the normalized dispersion relation is as follows (Shaaban et al., 2017):

[TABLE]

where $\tilde{k}=~{}kc/\omega_{p}$ is the normalized wave-number $k$, $\tilde{\omega}=\omega/\Omega_{p}$ and

[TABLE]

is the plasma dispersion function (Fried & Conte, 1961) of argument $\xi_{j}^{\pm}=\left(\omega\pm\Omega_{j}\right)/\left(ku_{j,\parallel}\right)$, with the imaginary part $\Im\left(\xi_{j}^{\pm}\right)>0$.

We solve the dispersion relation (2) numerically for the same four sets of plasma parameters used in our simulations and listed in Table 1 (see Section 3).

Figure 7 displays the unstable solutions of the firehose modes shown by the growth rates in the top panel and the wave frequencies in the bottom panel as functions of the normalized wavenumber. For case 1 the contribution of electrons is minimized by considering them isotropic, and therefore only the proton firehose instability is excited by the anisotropic protons with $A_{p}=0.1$, see black line. The growth rate and wave frequency of the unstable proton firehose mode obtained for case 1 serves as a reference for the rest of the results that are obtained for the other cases. The presence of anisotropic electrons can markedly alter the dispersive characteristics of the unstable proton firehose modes (Michno et al., 2014; Maneva et al., 2016; Shaaban et al., 2017). Therefore, in cases 2 and 4 (unstable proton firehose modes under the influence of the electron anisotropy for two different electron plasma betas), it is clear that anisotropic electrons with low plasma beta $\beta_{e,\parallel}=0.1$ have a negligible effect on the PFHI. For cases 1 and 2, only the PFHI is excited and the wave frequency is right-hand (RH) polarized (which we define as the frequency with positive $\tilde{\omega}_{r}$). For sufficiently large electron plasma beta, i.e. $\beta_{e,\parallel}=4.0$, the left-hand (LH) EFHI is excited at high normalized wavenumbers $\tilde{k}>1$ by anisotropic electrons with $A_{e}=0.1$. This leads also to the enhancement of the PFHI at low normalized wavenumbers $\tilde{k}<1$. The wave frequency plot in the bottom panel confirms the presence of both instabilities. The growing wave is RH-polarized at small wavenumbers and LH-polarized at high wavenumbers (the red curve changes the sign at at $\tilde{k}\sim 1$). In case 3 we obtain the unstable solutions of the EFH modes. With isotropic protons, the EFHI exhibits a larger growth rate and extends to larger wavenumbers compared to those obtained for case 4. Straightforward comparison between the results obtained for cases 3 and 4 shows clearly that anisotropic protons with $A_{p}<1.0$ have inhibiting effects on the EFHI, decreasing both growth rates and unstable wavenumbers. The corresponding wave frequency is negative along the bulk of the solution including the part of the normalized wavenumber space around the maximum of the growth rate. This confirms the LH polarization of the unstable EFH mode.

In order to compare our PIC simulations and the linear theory, we show in Figure 8 the fluctuations of the magnetic energy, normalized to the initial magnetic field $\bm{B}_{0}$, and compare their growth with relative growth rates obtained from the linear theory, for the four cases listed in Table 1. In all the simulations the system evolves from an initially unstable regime to the state of marginal stability, due to the interaction between waves and particles. We note that the very low level of magnetic fluctuations close to $t=0$, before the development of the instabilities, depends on whether or not electrons are initially isotropic. This is linked to the fact that we do not use any driving except the initial temperature anisotropies (i.e. we do not seed the instability). Indeed, in all our simulations the instability grows starting from the numerical noise. Different level of initial magnetic fluctuations does not affect the linear growth rate of the instability, as expected. The magnetic fluctuations grow with time as a signature of an exchange of energy between the magnetic field and the particles. The saturation is reached once the anisotropy in the electron or/and proton velocity distributions has almost entirely disappeared.

Time evolution of the magnetic energy fluctuations for the PFHI, is almost identical in cases 1 and 2, as indicated in Figure 8 with the black and the blue solid lines respectively. The magnetic energy fluctuations in both linear and nonlinear regimes of the instability are similar. The theoretical growth rate ($\gamma$) for these two sets of unstable conditions is shown in Figure 8 as a dashed black line corresponding to $\gamma=0.24\,\Omega_{p}$. Both cases lead to the same evolution of the magnetic energy fluctuations, as for small values of the electron plasma $\beta$ the electron firehose does not develop and the electron temperature anisotropy has no effect on the growth rate of the parallel proton firehose instability.

For the EFHI (case 3), the evolution of the magnetic energy fluctuations is shown in Figure 8 with the purple solid line. The maximum growth rate in this case is around three times higher than the PFHI growth rate. Indeed, the maximum growth rate according to the linear theory is $\gamma=0.77\,\Omega_{p}$ and is indicated in Figure 8 with the dashed purple line.

The main finding of this work is related to the evolution of the magnetic energy fluctuations for the PFHI + EFHI (case 4), displayed in Figure 8 with a red solid line. In this case, as shown above, the interplay between the electron and proton anisotropies significantly modifies the development of the instability. The growth rate of the parallel proton firehose instability is enhanced by the electron anisotropy in comparison with the pure proton instability (dashed black line in Figure 8). Indeed, this case is characterized by the linear growth rate $\gamma=0.54\,\Omega_{p}$ (dashed red line).

In order to verify that our simulations properly captured the linear stage of electron and proton parallel firehose instabilities and to have an additional comparison with the linear theory, Figure 9 shows the temporal evolution of the normalized wavenumber $k$ for cases 1, 3 and 4, obtained as fast Fourier transforms in space of the out-of-the-plane component of the magnetic field $B_{z}$. In Figure 9 one can remark the more rapid development of the electron firehose instability and its extension to higher wavenumbers compared to the other two cases. Furthermore, comparing the PFHI + EFHI case 4 and the PFHI case 1 it is possible to notice the earlier development of the firehose instability in the presence of the electron temperature anisotropy and the higher wavenumbers which characterize case 4 in comparison with case 1. These results are in qualitative agreement with those obtained via linear theory as shown in the top panel of Figure 7.

Figure 9 (c) shows also that other modes develop after the saturation of the fastest growing mode around $t=40000\ \omega_{p}^{-1}$ (cf. Figure 8). In particular, we note that a very energetic mode with $ck/\omega_{p}\sim 1$ grows around $t=65000\ \omega_{p}^{-1}$. At this time, the evolution of the oscillations of the magnetic field energy displays an increment (see red solid curve in Figure 8).

5 A test of oblique cases

Three-dimensional (3D) and two-dimensional (2D) full PIC simulations are required in order to obtain a more complete picture of the processes that take place in the coupled dynamics of parallel and oblique firehose instabilities, when the kinetic description for both species is retained. This kind of problems can hardly be modeled in multi-dimensional simulations with real mass ratio due to extremely challenging computational effort. In an attempt to test the coupling between electrons and protons in oblique cases, we performed 1D simulations at oblique angles with respect to the background magnetic field, which allow the development of oblique electron and proton firehose instabilities. Using the same parameters as described in Sections 2 and 3, we present a comparison between the PFHI case 1 and the PFHI + EFHI case 4 at $\theta=0^{\circ}$, $\theta=35^{\circ}$, $\theta=45^{\circ}$ and $\theta=60^{\circ}$, where $\theta$ is the angle between the directions of the magnetic field and the $x$-axis of the simulation box.

Figure 10 shows the temporal evolution of the magnetic energies for eight different runs, corresponding to PFHI + EFHI case 4 and PFHI case 1 for different orientations of the simulation box with respect to the initial background magnetic field. In this way we can compare the growth rates of parallel and oblique proton firehose instabilities and reveal how the presence of anisotropic electrons with $A_{e}<1$ affects the development of the proton firehose instability at each propagation angle. Each simulation starts from an unstable condition and the instability grows until the magnetic energy saturation, which corresponds to the quasi-isotropization of the proton VDF. In case 4, at oblique angles, the initial peak corresponds to the electron firehose and the later peak corresponds to the proton firehose.

Our result reported in Section 3, calculated at $\theta=0^{\circ}$, is confirmed to be valid also at oblique angles of propagation. Indeed, also for $\theta=35^{\circ}$ (black lines) and $\theta=45^{\circ}$ (red lines) the case 4, characterized by the anisotropy of both electrons and protons with $A_{j}<1$ (solid lines), presents a higher growth rate and a faster magnetic energy saturation than the case where only the protons are anisotropic (dashed lines). The difference in the growth rates becomes less evident at higher angles of propagation (blue curves in Figure 10). However, the impact of the electrons on the development of proton firehose instabilities is clearly present in the temporal evolution of proton temperature anisotropies. Indeed, the final proton temperature anisotropy is lower ($A_{p}$ is higher) in PFHI + EFHI case 4 compared to the PFHI case 1 at all the different propagation angles analyzed (see Figure 11).

We confirm that the parallel proton firehose instability has a lower threshold and a higher growth rate as compared to the oblique proton firehose instability (Hellinger & Matsumoto, 2000; Kasper et al., 2002; Hunana & Zank, 2017; Hunana et al., 2019).Although for these sets of parameters the dominant oblique electron firehose instability saturates on fast, electron scales and the electron temperature anisotropy is reduced much earlier than the proton temperature anisotropy, the interplay between the two species is still present and clearly affects the development of the proton oblique firehose instability.

Figure 11 shows the temporal evolution of proton temperature anisotropies in the PFHI + EFHI case 4 and PFHI case 1 at the same four angles $\theta$ listed above. In all the cases analyzed, the proton VDFs react to the increase of the magnetic energy, reducing their initial anisotropy until a stable quasi-isotropic state is reached. For each angle of propagation, the isotropization is faster when the electrons are also unstable with respect to the oblique firehose instability (solid lines), as compared to the case 1 where the electrons are initially isotropic (dashed lines). Moreover, the proton isotropization is faster in the parallel firehose instability than in the oblique firehose instability, due to the higher growth rate of the former. The electron temperature anisotropies are not shown in Figure 11 for the sake of clarity, since the electron isotropization in the PFHI + EFHI case 4 occurs on fast scales, while they remain isotropic in the PFHI case 1. However, it has to be pointed out that the interplay of instabilities at different angles is not captured in our 1D setup and has to be studied in more realistic 2D and 3D simulations.

6 Discussion and conclusions

In order to explain a low level of anisotropy recorded for both electrons and protons around 1 AU and beyond (Matteini et al., 2007), it has been suggested that electromagnetic micro-instabilities could play a fundamental role in limiting the difference among the diagonal elements of the pressure tensor (see e.g. the review by Marsch, 2006). Indeed, plasma electromagnetic instabilities in the solar wind are processes that regulate and limit the anisotropies in the particle velocity space, so the observed velocity distribution functions for both protons and electrons are at the margin of stability (Štverák et al., 2008), (Kasper et al., 2016). The firehose instability is thought to have a significant role in regulating electron and proton temperature anisotropies in the solar wind. Verifying this hypothesis requires fully kinetic simulations of plasma conditions in the solar wind.

We have studied the parallel firehose instability for protons and electrons via 1D fully kinetic simulations performed with ECsim, a semi-implicit PIC code. The main focus of this work is the interplay between anisotropic protons and electrons with $A_{j}<1$. Our non-linear fully kinetic simulations demonstrate that the electron temperature anisotropy, $A_{e}<1$, leads to a faster development of the proton parallel firehose instability, as suggested by the linear theory (Michno et al., 2014; Shaaban et al., 2017).

The presence of an electron temperature anisotropy has a limited impact on the level of the magnetic field fluctuations at saturation. The case characterized by the temperature anisotropy of both species and high electron plasma beta (PFHI + EFHI case 4) exhibits the final proton temperature anisotropy which is slightly lower compared to that obtained for the pure proton firehose instability (PFHI case 1), where only the protons are anisotropic. The final value of $A_{p}$ around $0.7$ is in agreement with the most probable level of proton anisotropy observed in the solar wind at 1 AU (Kasper et al., 2003). The proton firehose instability affects only marginally the electrons when they are stable with respect to the firehose instability. Indeed, in these cases the electron temperature anisotropy presents only small oscillations around the stable initial conditions.

The use of a semi-implicit PIC code has given us the opportunity to use a large simulation box that does not limit the number of unstable modes (more than 20 wave lengths of the fastest growing mode fit into the box). In addition we use a realistic mass ratio and a frequency ratio $\omega_{e}/\Omega_{e}=63.24$, that is relevant for the solar wind (Tong et al., 2019). The results of our fully kinetic simulations are in good agreement with the ones obtained with the linear theory. This is a confirmation that ECsim models the physics correctly even at relatively coarse resolutions.

An important limitation of our work is that it does not take into account the competition between parallel and oblique proton firehose instability. Our simulation captures only the parallel character of the analyzed unstable mode. We note that it is still an open question which instability limits the observed proton temperature anisotropies. On the one hand, hybrid simulations reveal that both parallel and oblique proton firehose instabilities are relevant in the solar wind context and show comparable growth rates for a wide range of parameters (Hellinger & Matsumoto, 2000). On the other hand, Wind/SWE observations (Kasper et al., 2002) indicate that there exists a phenomenological constraint on the proton temperature anisotropy which is compatible with the threshold of the parallel firehose instability. Despite taking into account only the parallel instability, our analysis is relevant for understanding the consequences of the combination of ion and electron temperature anisotropies in the linear and non-linear regimes of the instability. The study of the oblique instabilities goes beyond the scope of this work and will be a subject of further investigations. More realistic and more computationally demanding case to describe the full spectrum of firehose instabilities goes beyond the scope of this work as two- or maybe even three-dimensional simulations (see Howes, 2015) would be necessary.

This work focuses on the periodic proton firehose instability and aims at understanding how it is affected by an electron temperature anisotropy. Even if the oblique electron firehose instability has generally a lower threshold and higher growth rate than the parallel electron firehose instability, numerical simulations by Camporeale & Burgess (2008) demonstrate that quasi-parallel, propagating modes of the electron firehose instability are excited in the non-linear stage of the instability. In particular, Camporeale & Burgess (2008) show the shift of the most unstable modes toward smaller angles of propagation until the time when the oblique, non-propagating modes with the highest growth rates saturate and the electron firehose instability becomes driven by modes at any angle of propagation. Camporeale & Burgess (2008) state that over a certain period of time even exactly parallel modes (propagation angle of 0*∘*) dominate. We note that, as mentioned in Section 1, the three properties - being non-resonant (for electrons), parallel and propagating - are often used as definitions of the parallel firehose instability (Li & Habbal, 2000; Gary & Nishimura, 2003; Camporeale & Burgess, 2008), and all are well captured in our simulations.

Another limitation of our simulations is that we assumed initial bi-Maxwellian distributions for protons and electrons. The presence of suprathermal populations, which are ubiquitous in space plasmas (Marsch, 2006; Štverák et al., 2008; Berčič et al., 2019), also affects the instability thresholds and may alter the competition between parallel and oblique firehose modes (Lazar et al., 2011; Astfalk & Jenko, 2016; López et al., 2019; Shaaban et al., 2019c). Moreover, linear and non-linear investigations of the properties of parallel and oblique firehose instabilities have determined that the growth rates of these two instabilities are strongly affected by the presence of alpha particles and by their properties (Dasso et al., 2003; Hellinger & Trávníček, 2006; Ofman, 2019). More complex structures of particle velocity distribution functions, as well as more realistic plasma composition, introduce other sources of free energy for kinetic instabilities, which will be investigated in future work.

This work was supported by a PhD grant awarded by the Royal Observatory of Belgium to one of the authors (A. M.). These simulations were performed on the supercomputers SuperMUC (LRZ) and Marconi (CINECA) under PRACE and HPC-Europa3 allocations. Authors thank M. E. Innocenti for helpful discussions. A. N. Z. thanks the European Space Agengy (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme. S. M. S. acknowledges support by a FWO Postdoctoral Fellowship (Grant No. 12Z6218N).