PLASMA TURBULENCE AND KINETIC INSTABILITIES AT ION SCALES IN THE EXPANDING SOLAR WIND

Turbulence in magnetized weakly collisional space and astrophysical plasmas is a ubiquitous nonlinear phenomenon that allows energy transfer from large to small scales and, eventually, to plasma particles. Properties of plasma turbulence and its dynamics remain an open, challenging problem (Petrosyan et al. 2010; Matthaeus & Velli 2011). The solar wind constitutes a natural laboratory for plasma turbulence (Alexandrova et al. 2013; Bruno & Carbone 2013) since it offers the opportunity of its detailed diagnostics. Turbulence at large scales can be described by the magnetohydrodynamic (MHD) approximation, accounting for the dominant nonlinear coupling and for the presence of the ambient magnetic field that introduces a preferred direction (Boldyrev et al. 2011). Around particle characteristic scales, the plasma description has to be extended beyond MHD, and at these scales, a transfer of the cascading energy to particles is expected. The solar wind turbulence indeed likely energizes particles: radial profiles of proton temperatures indicate an important heating that is often comparable to the estimated turbulent energy cascade rate (MacBride et al. 2008; Cranmer et al. 2009; Hellinger et al. 2013). This energization proceeds through collisionless processes that may have feedback on turbulence. In the solar wind, the problem is further complicated by a radial expansion that induces an additional damping; turbulent fluctuations decrease due to the expansion as well as due to the turbulent decay. The expansion thus slows down the turbulent cascade (see Grappin et al. 1993; Dong et al. 2014). Furthermore, the characteristic particle scales change with radial distance affecting possible particle energization mechanisms.

Understanding of the complex nonlinear properties of plasma turbulence on particle scales is facilitated via a numerical approach (Franci et al. 2015a; Servidio et al. 2015). Direct kinetic simulations of turbulence show that particles are indeed on average heated by the cascade (Parashar et al. 2009; Markovskii & Vasquez 2011; Wu et al. 2013; Franci et al. 2015a); moreover, turbulence leads locally to complex anisotropic and nongyrotropic distribution functions (Valentini et al. 2014; Servidio et al. 2015). Furthermore, expansion naturally generates particle temperature anisotropies (Matteini et al. 2012). The anisotropic and nongyrotropic features may be a source of free energy for kinetic instabilities. In situ observations indicate the existence of apparent bounds on the proton temperature anisotropies that are consistent with theoretical kinetic linear predictions (Hellinger et al. 2006; Hellinger & Trávníček 2014). These linear predictions have, however, many limited assumptions (Matteini et al. 2012; Isenberg et al. 2013); in particular, they assume a homogeneous plasma that is at odds with the presence of turbulent fluctuations. On the other hand, the observed bounds on the proton temperature anisotropy (and other plasma parameters) and enhanced magnetic fluctuations near these bounds (Wicks et al. 2013; Lacombe et al. 2014) indicate that these kinetic instabilities are active even in the presence of turbulence.

In this Letter, we directly test the relationship between proton kinetic instabilities and plasma turbulence in the solar wind using a hybrid expanding box model that allows us to study self-consistently physical processes at ion scales. In the hybrid expanding box model, a constant solar wind radial velocity v_sw is assumed. The radial distance R is then $R={R}_{0}(1+t/{t}_{{\rm{e}}0})$, where R₀ is the initial position and ${t}_{{\rm{e}}0}={R}_{0}/{v}_{\mathrm{sw}}$ is the initial value of the characteristic expansion time ${t}_{{\rm{e}}}=R/{v}_{\mathrm{sw}}={t}_{{\rm{e}}0}(1+t/{t}_{{\rm{e}}0}).$ Transverse scales (with respect to the radial direction) of a small portion of plasma, comoving with the solar wind velocity, increase ∝ R. The expanding box uses these comoving coordinates, approximating the spherical coordinates by the Cartesian ones (Liewer et al. 2001; Hellinger & Trávníček 2005). The model uses the hybrid approximation where electrons are considered as a massless, charge-neutralizing fluid and ions are described by a particle-in-cell model (Matthews 1994). Here, we use the two-dimensional (2D) version of the code, fields and moments are defined on a 2D x–y grid 2048 × 2048, and periodic boundary conditions are assumed. The spatial resolution is Δx = Δy = 0.25d_p0, where ${d}_{{\rm{p}}0}={v}_{{\rm{A}}0}/{{\rm{\Omega }}}_{{\rm{p}}0}$ is the initial proton inertial length (v_A0: the initial Alfvén velocity, Ω_p0: the initial proton gyrofrequency). There are 1024 macroparticles per cell for protons that are advanced with a time step ${\rm{\Delta }}t=0.05/{{\rm{\Omega }}}_{{\rm{p}}0}$, while the magnetic field is advanced with a smaller time step ${\rm{\Delta }}{t}_{B}={\rm{\Delta }}t/10.$ The initial ambient magnetic field is directed along the radial z-direction, perpendicular to the simulation plane ${{\boldsymbol{B}}}_{0}=(0,0,{B}_{0})$, and we impose a continuous expansion in the x- and y-directions. Due to the expansion, the ambient density and the magnitude of the ambient magnetic field decrease as $\bar{n}\propto \bar{B}\propto {R}^{-2}$ (the proton inertial length d_p increases ∝ R; the ratio between the transverse sizes and d_p remains constant; the proton gyrofrequency Ω_p decreases as ∝R⁻²). A small resistivity η is used to avoid accumulation of cascading energy at grid scales; initially, we set $\eta ={10}^{-3}{\mu }_{0}{v}_{{\rm{A}}0}^{2}/{{\rm{\Omega }}}_{{\rm{p}}0}$ (μ₀ being the magnetic permittivity of vacuum) and η is assumed to be $\propto \bar{n}.$ The simulation is initialized with an isotropic 2D spectrum of modes with random phases, linear Alfvén polarization ($\delta {\boldsymbol{B}}\perp {{\boldsymbol{B}}}_{0}$), and vanishing correlation between magnetic and velocity fluctuation. These modes are in the range 0.02 ≤ kd_p ≤ 0.2 and have a flat one-dimensional (1D) power spectrum with rms fluctuations = 0.24 B₀. For noninteracting zero-frequency Alfvén waves, the linear approximation predicts $\delta {B}_{\perp }\propto {R}^{-1}$ (Dong et al. 2014). Protons initially have the parallel proton beta ${\beta }_{{\rm{p}}\parallel }=0.8$ and the parallel temperature anisotropy ${A}_{{\rm{p}}}={T}_{{\rm{p}}\perp }/{T}_{{\rm{p}}\parallel }=0.5$ as typical proton parameters in the solar wind in the vicinity of 1 AU (Hellinger et al. 2006; Marsch et al. 2006). Electrons are assumed to be isotropic and isothermal with β_e = 0.5 at t = 0.

The initial random fluctuations rapidly relax and a turbulent cascade develops. Figure 1 shows the evolution of the 1D power spectral density (PSD) ${P}_{{\boldsymbol{B}}\perp }={P}_{{\boldsymbol{B}}\perp }(k)$ of the magnetic field ${{\boldsymbol{B}}}_{\perp }$ perpendicular to ${{\boldsymbol{B}}}_{0}.$ On large scales, the initial flat spectrum evolves to a power law. This large-scale power law remains clearly visible until $t\sim 0.7{t}_{{\rm{e}}0}$, although its slope slowly varies in time, passing from about −3/2 to −5/3 (these estimated slopes are, however, quite sensitive to the chosen range of wave vectors). The variation of large-scale slopes (kd_p ≲ 1) is likely connected with the decay of the large-scale fluctuations due to the cascade and the expansion as the inertial range is likely quite narrow. This problem is beyond the scope of the present Letter and will be a subject of future work (note that a similar steepening is also observed in MHD expanding box simulations; see Dong et al. 2014); this Letter is mainly focused on ion scales.

Zoom In Zoom Out Reset image size

Around kd_p ∼ 1 there is a smooth transition in ${P}_{{\boldsymbol{B}}\perp }$ separating the large-scale power-law slope and a steeper slope at sub-ion scales (Franci et al. 2015b). The PSD amplitudes decay in time partly due to the expansion and partly due to the turbulent damping. Note that there are some indications that the position of the transition shifts to smaller kd_p with time/radial distance (compare the blue, green, orange, and red curves in Figure 1); a similar trend is observed for the proton gyroradius since it increases only slightly faster than d_p. At later times, the fluctuating magnetic energy is enhanced at ion scales around kd_p ∼ 0.4 ÷ 1 (compare the red and black curves in Figure 1); this indicates that some electromagnetic fluctuations are generated at later times of the simulation.

Figure 2 summarizes the evolution of the simulated system, which goes through three phases. During the first phase, the system relaxes from the initial conditions and turbulence develops; the level of magnetic fluctuations increases at the expense of proton velocity fluctuations. The fluctuating magnetic field $\delta {B}_{\perp }/\bar{B}$ reaches the maximum at about $t\sim 0.008{t}_{{\rm{e}}0}.$ During this phase, a parallel current j_z is generated; $\langle {j}_{z}^{2}\rangle$ normalized to ${\bar{B}}^{2}/{d}_{p}^{2}$ reaches a maximum at t ∼ 0.035t_e0 indicating the presence of a well-developed turbulent cascade (Mininni & Pouquet 2009; Valentini et al. 2014). After that, the system is dominated by a decaying turbulence; the fluctuating magnetic field initially decreases faster than $\bar{B}$ until about 0.3t_e0.

Zoom In Zoom Out Reset image size

During the second phase, protons are heated. For negligible heat fluxes, collisions, and fluctuations, one expects the double adiabatic behavior or CGL (Chew et al. 1956; Matteini et al. 2012): the parallel and perpendicular temperatures (with respect to the magnetic field) are expected to follow ${T}_{{\rm{p}}\perp }\propto \bar{B}$ and ${T}_{{\rm{p}}\parallel }=\mathrm{const}.,$ respectively. T_p⊥ decreases slower than $\bar{B}$ during the whole simulation; protons are heated in the perpendicular direction while in the parallel direction the heating lasts till about t ∼ 0.25t_e0, whereas afterward protons are cooled. The parallel and perpendicular heating rates could be estimated as (see Verscharen et al. 2015)

$$\begin{eqnarray}&&{Q}_{\parallel }=\displaystyle \frac{{\bar{n}}^{3}{k}_{B}}{{\bar{B}}^{2}}\displaystyle \frac{d}{{dt}}\;\left(\displaystyle \frac{{T}_{\parallel }{\bar{B}}^{2}}{{\bar{n}}^{2}}\right)\ \mathrm{and}\ {Q}_{\perp }=\bar{n}{k}_{B}\bar{B}\displaystyle \frac{d}{{dt}}\;\left(\displaystyle \frac{{T}_{\perp }}{\bar{B}}\right).\end{eqnarray} \tag{ 1 }$$

A more detailed analysis indicates that between t = 0.1t_e0 and t = 0.7t_e0 the parallel heating rate ${Q}_{\parallel }$ smoothly varies from about 0.2 Q_e and $-0.2{Q}_{{\rm{e}}},$ whereas ${Q}_{\perp }$ is about constant at ∼0.2 Q_e; here, ${Q}_{{\rm{e}}}=\bar{n}{k}_{B}T/{t}_{{\rm{e}}}.$ In total, protons are heated until t ∼ 0.7t_e0, and the heating reappears near the end of the simulation t ≳ 0.95t_e0. Note that the perpendicular heating rate is a nonnegligible fraction of that observed in the solar wind, where ${Q}_{\perp }\approx 0.6{Q}_{{\rm{e}}}$ (Hellinger et al. 2013); however, the proton heating in 2D hybrid simulations is typically quite sensitive to the used electron equation of state (Parashar et al. 2014) and also to the resistivity and the number of particles per cell used (Franci et al. 2015a). The turbulent heating is, however, not sufficient to overcome the expansion-driven perpendicular cooling as in the solar wind (Matteini et al. 2007). During the third phase, t ≳ 0.7t_e0; there is an enhancement of the parallel cooling and perpendicular heating that cannot be ascribed to the effect of the turbulent activity. For a large parallel proton temperature anisotropy, a firehose instability is expected. The presence of such an instability is supported by the fact that the fluctuating magnetic field increases (with respect to the linear prediction), suggesting a generation of fluctuating magnetic energy at the expense of protons. To analyze the role of different processes in the system, we estimate their characteristic times (Matthaeus et al. 2014). The bottom panel of Figure 2 compares the turbulent nonlinear eddy turnover time ${t}_{\mathrm{nl}}={k}^{-3/2}{({P}_{{\boldsymbol{B}}\perp }(k)/{\mu }_{0}{m}_{{\rm{p}}})}^{-1/2}$ at kd_p = 1 (see Matthaeus et al. 2014; the expansion time t_e, and the linear time t_l of the oblique firehose (Hellinger & Matsumoto 2000, 2001) estimated as ${t}_{{\rm{l}}}=1/{\gamma }_{{\rm{m}}},$ where γ_m is the maximum growth rate calculated from the average plasma properties in the box assuming bi-Maxwellian proton velocity distribution functions (Hellinger et al. 2006). The expansion time t_e is much longer than t_nl at kd_p = 1 (as well as at the injection scales). The expanding system becomes theoretically unstable with respect to the oblique firehose around t ∼ 0.47t_e0 but clear signatures of a fast proton isotropization and of a generation of enhanced magnetic fluctuations appear later t ≳ 0.7t_e0. This is about the time when the linear time becomes comparable to the nonlinear time at ion scales. After that, ${t}_{{\rm{l}}}{{\rm{\Omega }}}_{{\rm{p}}}$ slightly increases as a result of a saturation of the firehose instability, whereas ${t}_{\mathrm{nl}}{{\rm{\Omega }}}_{{\rm{p}}}$ at kd_p = 1 is about constant (note that Ω_p decreases as R⁻²). This may indicate that the instability has to be fast enough to compete with turbulence; however, the 2D system has strong geometrical constraints. Also the stability is governed by the local plasma properties. Figure 3 shows the evolution of the system in the plane $({\beta }_{{\rm{p}}\parallel },{A}_{{\rm{p}}}).$ During the evolution, a large spread of local values in the 2D space $({\beta }_{{\rm{p}}\parallel },{A}_{{\rm{p}}})$ develops. Between t ≃ 0.1t_e0 and t ≃ 0.65t_e0, the average quantities evolve in time following $\langle {A}_{{\rm{p}}}\rangle \propto \langle {\beta }_{{\rm{p}}\parallel }{\rangle }^{-0.86}.$ This anticorrelation is qualitatively similar to in situ Helios observations between 0.3 and 1 AU (Matteini et al. 2007). During the third stage, when the strong parallel temperature anisotropy is reduced, both local and average values of ${\beta }_{{\rm{p}}\parallel }$ and A_p appear to be bounded by the linear marginal stability conditions of the oblique firehose (Hellinger & Trávníček 2008), although relatively large theoretical growth rates ${\gamma }_{{\rm{m}}}\sim 0.1{{\rm{\Omega }}}_{{\rm{p}}}$ are expected.

The evolution in the real space is shown in Figure 4, which shows the magnitude of the perpendicular fluctuating magnetic field δB_⊥and the proton temperature anisotropy A_p at different times (see also the animation that combines the evolution in Figures 3 and 4). The modes with initially random phases rapidly form vortices and current sheets. Despite the overall turbulent heating, a strong temperature anisotropy ${T}_{{\rm{p}}\perp }\lt {T}_{{\rm{p}}\parallel }$ develops owing to the expansion and a firehoselike activity develops in the form of localized waves/filaments with enhanced $\delta {B}_{\perp }.$ These fluctuations appear in regions between vortices where ${B}_{\perp }/\bar{B}$ is enhanced, i.e., in places where the angle between the simulation plane and the local magnetic field θ_B ($\approx \mathrm{arccos}{B}_{\perp }/\bar{B}$) is less oblique (reaching $\sim {60}^{{\rm{o}}}$ and below). This is in agreement with the theoretical expectations, while the oblique firehose is unstable for moderately oblique wavevectors with respect to the magnetic field near the threshold; farther away from the threshold, the unstable modes become more oblique (see Figure 3), and these oblique angles are locally available between vortices where we observe the enhanced level of magnetic fluctuations. These geometrical factors may be responsible for the late appearance of the instability (but constraints on the instability timescales imposed by the turbulent nonlinearities are likely also important). The localized wave packets are Alvénic, have wavelengths of the order of 10 d_p, and propagate with a phase velocity of about 0.1 v_A, in agreement with the expectation for the nonlinear phase of the oblique firehose. Furthermore, the parallel temperature anisotropy is strongly reduced in their vicinity. These Alvénic wave packets are responsible for the enhanced level of the magnetic PSD at ion scales seen in Figure 1.

Zoom In Zoom Out Reset image size

For a linear instability, it is expected that the magnetic fluctuations increase exponentially in time during its initial phase (except when the growth time is comparable to the expansion time; see Tenerani & Velli 2013). Figure 2, however, shows that the overall magnetic fluctuations $\delta {B}_{\perp }/\bar{B}$ (with respect to the linear prediction), and $\delta {B}_{\parallel }/\bar{B}$ increase rather slowly (secularly) in time for t ≳ 0.7t_e. This behavior is expected for a long time evolution in a forced system after saturation (see Matteini et al. 2006; Rosin et al. 2011; Kunz et al. 2014). An additional analysis indicates that the expected exponential growth is indeed seen in the simulation, but only locally both in space and time. This exponential growth is obscured by the turbulent fluctuations; furthermore, it is blurred out due to the averaging over the simulation box in the global view of Figure 2.

On a microscopic level, the firehose activity leads to an efficient scattering from parallel to perpendicular direction of protons in the velocity space. Figure 5 shows the evolution of the proton velocity distribution function $f=f({v}_{\parallel },{v}_{\perp })$ averaged over the simulation box. While turbulence leads locally to complex proton distribution functions (see Valentini et al. 2014; Servidio et al. 2015), the average proton distribution function during the first two phases remains relatively close to a bi-Maxwellian shape (Figure 5, top panel). During the third phase, there appear to be clear signatures of the cyclotron diffusion (for protons with ${v}_{\parallel }\gtrsim {v}_{{\rm{A}}}$) as expected for the oblique firehose instability (Hellinger & Trávníček 2008).

Zoom In Zoom Out Reset image size

Using 2D hybrid simulations, we investigated the evolution of turbulence in a slowly expanding plasma. The numerical model shows that the turbulent heating is not sufficient to overcome the expansion-driven cooling and that the oblique firehose becomes active for a sufficiently large parallel proton temperature anisotropy and for sufficiently oblique angles of propagation.

While the modeled expansion is about 10 times faster than in the solar wind, the ratio between the expansion and the nonlinear eddy turnover timescales is quite realistic: ${t}_{{\rm{e}}}/{t}_{\mathrm{nl}}\approx 1000$ at kd_p = 1 for t ≳ 0.7t_e, which is about 4 times smaller than that of the solar wind with similar plasma parameters at 1 AU (Matthaeus et al. 2014). Note also that a similar evolution is observed for many different plasma and expansion parameters.

In the present case, both turbulence and the 2D geometry constraints strongly affect the firehose instability, and there are indications that firehose has an influence on turbulence (the mixed third-order structure functions are enhanced due to the firehose activity suggesting a stronger cascade rate; Verdini et al. 2015). The problem of the interaction between turbulence and kinetic instabilities requires further work. Three-dimensional simulations are needed to investigate the interplay between turbulence and instabilities as the most unstable modes usually are parallel or moderately oblique with respect to the ambient magnetic field. In the present case, the parallel firehose (Gary et al. 1998; Matteini et al. 2006) would be the dominant instability, but the 2D constraints strongly inhibit it. On the other hand, numerical simulations indicate that the oblique firehose plays an important role in constraining the proton temperature anisotropy in the expanding solar wind even in the case when the parallel firehose is dominant (Hellinger & Trávníček 2008). Nevertheless, the present work clearly demonstrates for the first time that kinetic instabilities may coexist with strong plasma turbulence and bound the plasma parameter space.

The authors wish to acknowledge valuable discussions with Marco Velli. The research leading to these results has received funding from the European Commission's 7th Framework Programme under grant agreement #284515 (project-shock.eu). P.H. and P.M.T. acknowledge GACR grant 15-10057S and the projects RVO:67985815 and RVO:68378289. L.M. acknowledges UK STFC grant ST/K001051/1. A.V. acknowledges the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office (IAP P7/08 CHARM).