The Role of Parametric Instabilities in Turbulence Generation and Proton Heating: Hybrid Simulations of Parallel-propagating Alfvén Waves

C. A. González; A. Tenerani; M. Velli; P. Hellinger

doi:10.3847/1538-4357/abbccd

1. Introduction

Collisionless or weakly collisional turbulent plasmas are typically found in space and astrophysical environments. As is the case of the heliosphere and the solar wind, the outflow of plasma that is continually emitted by the Sun and permeates our solar system. Large-amplitude fluctuations in the plasma velocity and magnetic field, known as Alfvénic fluctuations, are commonly observed in the solar wind. Such fluctuations are almost incompressible, and they display the typical velocity-magnetic field correlation that characterizes Alfvén waves propagating away from the Sun (Coleman 1967; Belcher and Davis 1971). In spite of such a high degree of correlation, Alfvénic fluctuations in the solar wind are characterized by a well-developed power-law spectrum that dominates the low-frequency range of the solar wind fluctuation energy (Bavassano et al. 1982; Horbury et al. 2005; Bruno & Carbone 2013). It is thought that these Alfvénic fluctuations might be generated near the Sun and that they may contribute to coronal heating and solar wind acceleration (Velli 1993; Erdélyi & Fedun 2007; Verdini et al. 2009), problems that are still under debate in the community.

In situ observations support the idea that dissipation of turbulent fluctuations might contribute significantly to plasma heating (Smith et al. 2001; Bruno & Carbone 2013; Hellinger et al. 2013). However, internal energy generation involves different channels, such as resonant (Hollweg & Isenberg 2002; Chen et al. 2019) and nonresonant (stochastic) wave–particle interactions (Chen et al. 2001; Chandran et al. 2010; Cranmer 2014), magnetic reconnection within coherent structures or scattering by current sheets (Dmitruk et al. 2004; Drake et al. 2009; Parashar et al. 2009; Matthaeus & Velli 2011; Servidio et al. 2011; Karimabadi et al. 2013; Zhdankin et al. 2013; Isliker et al. 2017; Pisokas et al. 2018), among others, and the very nature of the dissipation process(es) is still puzzling. In this regard, the evolution of proton temperature in the solar wind shows a strong departure double adiabatic expansion, and preferential particle heating in the perpendicular direction to the local magnetic field is typically observed. In addition, the proton distribution function displays many nonthermal features such as a secondary proton population with a drift velocity of the order of the local Alfvén speed (Marsch 2006). Interestingly, kinetic simulations have shown that a field-aligned proton beam may form self-consistently through the decay of an initial large-amplitude Alfvénic fluctuation (Araneda et al. 2008; Matteini et al. 2010).

Alfvén waves of arbitrary amplitude with constant total pressure are known to provide an exact solution to the compressible magnetohydrodynamic (MHD) system in a homogeneous plasma, in that nonlinearities are turned off and there no couplings with compressible modes. However, such a dynamical system is linearly unstable to parametric instabilities and large-amplitude Alfvén waves are known to decay into compressible and secondary Alfvénic modes through three- or four-wave resonances that lead to a variety of parametric instabilities, depending on the plasma beta and dispersive effects, such is the case of parametric decay (Galeev & Oraevskii 1973; Derby 1978), modulational, and beat instabilities (Sakai & Sonnerup 1983; Wong & Goldstein 1986; Jayanti & Hollweg 1993; Nariyuki & Hada 2007). Parametric instabilities of Alfvén waves (or of a spectrum of Alfvén waves) have been widely studied over the years through theoretical approaches (Goldstein 1978; Jayanti & Hollweg 1993; Malara & Velli 1996), and numerical simulations adopting both MHD (Ghosh et al. 1994; Malara et al. 2000; Del Zanna et al. 2001; Tenerani et al. 2017) and kinetic models (Terasawa et al. 1986; Matteini et al. 2010; Nariyuki et al. 2012; Verscharen et al. 2012), although most often in one-dimensional setups. In particular, the traditional parametric decay instability has attracted much attention over the years in the context of both turbulence and plasma heating. This type of decay is most efficient at low values of the plasma beta and it essentially involves the decay of a pump Alfvén wave into a lower-frequency reflected Alfvén wave and a forward sound wave. For this reason, parametric decay remains an appealing process because it provides a natural mechanism for the production of reflected modes, which is essential for the triggering of a turbulent cascade. Indeed, recently it has been proposed as a viable mechanism to initiate the turbulent cascade in the solar wind acceleration region (Chandran 2018; Réville et al. 2018), while global MHD simulations of the solar wind have also shown that the parametric decay instability can contribute substantially to solar wind heating and acceleration, thanks to the generation of compressible modes that, in the absence of kinetic effects, naturally steepen into shocks (see, e.g., Shoda et al. 2019). The traditional parametric decay has been also invoked as a possible source for the generation of inward modes and solar wind turbulence in the inner heliosphere, where an increasing content of reflected waves (cross-helicity) and an evolving turbulent spectrum are observed with increasing heliocentric distance (Bavassano et al. 2000). However, expansion effects are known to inhibit its development, essentially because the parametric decay process is strongly suppressed as the plasma beta increases at larger heliocentric distances (Tenerani & Velli 2013, 2020; Del Zanna et al. 2015). Temperature anisotropies can destabilize the parametric decay at values of the plasma beta approaching unity and above, but the anisotropy in the solar wind is not large enough to affect the instability significantly (Tenerani et al. 2017).

Despite much work on parametric instabilities, less attention has been devoted to kinetic effects in multidimensions. The multidimensional nature of parametric instabilities of a parallel-propagating Alfvén wave was first investigated via the two-fluid linear theory by Kuo et al. (1988) and later in the work by Viñas & Goldstein (1991, 1992) where it was shown that the oblique propagation of the daughter waves allows for additional parametric instabilities depending on the angle of the density perturbation with respect to the mean magnetic field. Previous numerical studies showed that while oblique modes naturally emerge when the pump wave itself is in oblique propagation or in two-dimensional turbulence (as observed, for example, in Matteini et al. 2010; Primavera et al. 2019), perpendicular and quasi-perpendicular modes can grow as the result of a different decay process of an Alfvén wave in parallel propagation, known as the filamentation and the magnetoacoustic instability, respectively. Such highly oblique modes have been reported previously in numerical simulations (Gao et al. 2013; Comişel et al. 2018, 2019).

In this paper we revisit the stability of Alfvén waves in parallel propagation using 1D, 2D, and 3D hybrid simulations to explore the combined effect of multidimensionality and kinetic proton physics at different values of the plasma beta and pump wave amplitude. We consider left-handed circularly polarized large-amplitude Alfvén waves, and we investigate how the decay process and its saturation and nonlinear stages depend on the pump wave amplitude, plasma beta, and dimensionality. We find that the overall decay process involves a superposition of modes due to parametric and filamentation instability that survives at values of the plasma beta well above unity. The instability triggers both a turbulent cascade in the perpendicular direction and a wave energy conversion process that ultimately leads to the formation of a field-aligned proton beam at the Alfvén speed, regardless of the plasma beta, and that appears to be associated with a strong particle heating. While particle heating is affected by the pump wave amplitude, it is surprisingly observed whenever the decay occurs, regardless of the dimensionality and of the plasma beta. This result thus suggests that the particle heating we observe in our simulations is predominantly a one-dimensional process driven by the decay of the wave.

This paper is organized as follows: In Section 2 we present the quasi-neutral hybrid model and the numerical setup that we have employed in this study. Section 3 is divided in three parts. In Section 3.1 we describe the global dynamics of the instability for different plasma beta, wave amplitude, and dimensionality of the problem. The spectral properties of the electromagnetic field are presented in Section 3.2 where we focus on the turbulent properties at different scales. Finally, in Section 3.3 we discuss the effect of the proton/electron beta and the wave amplitude on the proton heating and acceleration, and address the problem of wave–particle interactions by proposing a possible mechanism to explain the observed features. In Section 4 we summarize our results.

2. Model and Simulation Setup

We have employed a hybrid model where electrons are treated as a massless and isothermal neutralizing fluid, while the proton dynamics is described by the Vlasov–Maxwell equations (Equation (1)). The coupling with the electromagnetic fields is given by the low-frequency and nonrelativistic Maxwell's equations, where quasi-neutrality ( ${n}_{i}={n}_{e}=n$ ) is assumed and the electric field is determined via the generalized Ohm's law:

$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{i}}{\partial t}+{\boldsymbol{v}}\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{r}}}+\displaystyle \frac{e}{{m}_{i}}\left({\boldsymbol{E}}+\displaystyle \frac{{\boldsymbol{v}}}{c}\times {\boldsymbol{B}}\right)\cdot \displaystyle \frac{\partial {f}_{i}}{\partial {\boldsymbol{v}}}=0,\end{eqnarray} \tag{ 1a }$

$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c\ {\rm{\nabla }}\times {\boldsymbol{E}},\ \ \ \ {\boldsymbol{J}}=\displaystyle \frac{c}{4\pi }{\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 1b }$

$\begin{eqnarray}&&{\boldsymbol{E}}+\displaystyle \frac{{{\boldsymbol{u}}}_{i}}{c}\times {\boldsymbol{B}}=-\displaystyle \frac{{k}_{{\rm{B}}}{T}_{e}{\rm{\nabla }}n}{{en}}+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{{en}}+\eta {\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 1c }$

with c the speed of light, e the electron charge, k_B the Boltzmann constant, and T_e is the electron temperature. The proton number density n and the proton bulk velocity ${{\boldsymbol{u}}}_{i}$ are computed from the moments of the distribution function ( $n=\int f({\boldsymbol{r}},{\boldsymbol{v}},t)d{\boldsymbol{v}}$ and $n{{\boldsymbol{u}}}_{i}=\int {\boldsymbol{v}}f({\boldsymbol{r}},{\boldsymbol{v}},t)d{\boldsymbol{v}}$ , respectively). In this work we made use of the CAMELIA code (see, e.g., Franci et al. 2018), which is a hybrid particle-in-cell code that uses the current advance method (Matthews 1994) and Boris scheme for the particle pusher, with good stability and long-term accuracy.

The numerical setup consists of a large-amplitude, large-scale Alfvén wave propagating along the mean magnetic field, ${{\boldsymbol{B}}}_{0}$ , that we take along the x-axis. Periodic boundary conditions are imposed in all directions of the computational box. Lengths are normalized to the proton inertial length ${d}_{i}=c/{\omega }_{p}$ with ${\omega }_{p}={(4\pi {{ne}}^{2}/{m}_{i})}^{1/2}$ the proton plasma frequency. Time is expressed in units of the inverse of proton gyrofrequency ${{\rm{\Omega }}}_{{ci}}^{-1}={({{eB}}_{0}/{m}_{i}c)}^{-1}$ and velocities are normalized to the Alfvén speed ${v}_{{\rm{A}}}={B}_{0}/{(4\pi {{nm}}_{i})}^{1/2}$ . The plasma beta for both ions and electrons is defined as ${\beta }_{p,e}=8\pi {{nk}}_{{\rm{B}}}{T}_{p,e}/{B}_{0}^{2}$ . Since in almost all simulations we have ${\beta }_{p}={\beta }_{e}$ , we will just use the symbol β to indicate the proton beta, unless otherwise specified. We have included the resistive term in the generalized Ohm's law to improve energy conservation by avoiding energy accumulation at the grid scales. The resistive coefficient is defined in units of $4\pi {\omega }_{p}^{-1}$ and the associated length scale is chosen to be greater than the grid size but smaller than any other scale of interest (i.e., smaller than the proton inertial length or proton gyroradius depending on the plasma beta).

We initialize the system using an isotropic homogeneous plasma with a uniform particle density and proton velocities randomly distributed with a Maxwellian distribution function at a given temperature T_p and a fixed number of particles per cell (npp). The initial pump Alfvén wave is initialized with a wavenumber ${n}_{0}=4$ and wavevector ${k}_{0}=2\pi {n}_{0}/L$ , L being the box size in units of proton inertial length (we use a square or cube box with equal sides), that satisfies the condition $\delta {\boldsymbol{u}}=-({\omega }_{0}/{k}_{0})\delta {\boldsymbol{b}}$ , with $| \delta {\boldsymbol{b}}| =\delta {b}_{0}$ the amplitude of the pump wave normalized to the mean magnetic field magnitude B₀. The wave frequency is determined from the normalized dispersion relation ${k}_{0}^{2}={\omega }_{0}^{2}/(1-{\omega }_{0})$ for left-handed circularly polarized waves. The initial Alfvén wave is given by $\delta {b}_{z}=\delta {b}_{0}\cos ({k}_{0}x)$ and $\delta {b}_{y}=-\delta {b}_{0}\sin ({k}_{o}x)$ . The box sizes adopted in all the simulations presented in this paper have $L=128{d}_{i}$ , and the pump wave is weakly dispersive with a wavenumber ${k}_{0}{v}_{{\rm{A}}}/{{\rm{\Omega }}}_{{ci}}=0.196$ . A summary of the numerical and plasma parameters adopted in this work can be found in Table 1.

Table 1. Numerical Parameters for the Simulations Presented in This Paper

Run	Δx	Δt	ppc	β_p	β_e	δb₀	η
A1-1D	0.25	0.025	10000	0.25	0.25	1.0	0.002
A2-1D	0.25	0.025	10000	0.5	0.5	1.0	0.002
A3-1D	0.25	0.025	10000	2.0	2.0	1.0	0.002
B-1D	0.25	0.025	10000	0.5	0.5	0.5	0.002
B-1D	0.25	0.025	10000	0.5	0.5	0.25	0.002
A1-2D	0.0625	0.005	1000	0.25	0.25	1.0	0.0004
A2-2D	0.0625	0.00625	1000	0.5	0.5	1.0	0.0004
A3-2D	0.0625	0.01	1000	2.0	2.0	1.0	0.0004
B1-2D	0.25	0.025	1024	0.5	0.5	1.0	0.002
B2-2D	0.25	0.025	1024	0.75	0.75	1.0	0.002
B3-2D	0.25	0.025	1024	1.0	1.0	1.0	0.002
B4-2D	0.25	0.025	1024	2.0	2.0	1.0	0.002
B5-2D	0.25	0.025	1024	4.0	4.0	1.0	0.002
C1-2D	0.25	0.025	1024	0.5	0.5	0.5	0.002
C2-2D	0.25	0.025	1024	0.5	0.5	0.25	0.002
D-2D	0.25	0.025	1024	0.0	2.0	1.0	0.002
A3-3D	0.5	0.05	512	2.0	2.0	1.0	0.004

Download table as: ASCII Typeset image

3. Results

3.1. Global Dynamics

In Figure 1, top panel, we show the time evolution of the kinetic, magnetic, and thermal energy of three representative simulations that we use as a reference to summarize the main properties of the dynamical evolution of the system. In particular, we show results for runs A3-1D, A3-2D, and A3-3D for the 1D, 2D, and 3D cases, respectively, with $\beta =2$ and $\delta {b}_{0}=1$ . The bottom panel shows the parallel and perpendicular proton temperature evolution (black and red colors, respectively) and the temperature anisotropy ${T}_{\perp }/{T}_{\parallel }$ (green color) for the same simulations. We have introduced the parallel and perpendicular temperatures defined in terms of the decomposition of the pressure tensor according to the direction of the total magnetic field as ${p}_{\parallel }\,=\,{\boldsymbol{p}}:\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}$ and ${p}_{\perp }={\boldsymbol{p}}:({\mathbb{I}}-\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}})/2.$ The pressure tensor ${\boldsymbol{p}}=\int {(u-{v}_{p})}_{i}{(u-{v}_{p})}_{j}f({\boldsymbol{r}},{\boldsymbol{v}},t)d{\boldsymbol{v}}$ is obtained from the particle velocity distribution, and $\hat{{\boldsymbol{b}}}={\boldsymbol{B}}/\parallel {\boldsymbol{B}}\parallel$ is the direction of the total magnetic field.

Three different stages can be identified during the evolution: initially, the wave propagates without significant dispersion, the kinetic and magnetic energy oscillate around a mean value, while the proton temperature remains constant. This oscillation is possibly due to fact that the initial condition is not an exact solution to the Vlasov–Maxwell equation (Sonnerup & Su 1967). After this initial stage, the pump wave decays by conveying its energy to the particles (at $t\sim 250{{\rm{\Omega }}}_{{ci}}$ ), resulting in an increase of the overall thermal energy. Finally, the saturation of the instability slows down the particle energization process, and the system achieves a steady-state condition with almost constant kinetic, magnetic, and thermal energy ( $t\sim 400{{\rm{\Omega }}}_{{ci}}$ ). Note that for the 1D simulation shown here there is no proton heating at all. That is because for $\beta =2$ the decay instability is suppressed in 1D, and therefore the wave is not disrupted. The total energy of the system is not perfectly conserved during the simulations, in fact, there is a relative error of the order of $3 \%$ and some numerical heating is present in the simulations.

The evolution of the rms of density fluctuations and of the average proton temperatures for different plasma beta and wave amplitudes are presented in Figure 2, where the decay of the pump wave is marked by the rapid increase of density fluctuations and of the overall temperature. Here, the total temperature is defined as $T=({T}_{\parallel }+2{T}_{\perp })/3$ , and ${\rm{\Delta }}T=T(t)-T(t=0)$ represents the net change of the total temperature from the initial value. Results for 1D simulations are also plotted as a reference (dashed lines).

The left panels of Figure 2 display results for different plasma beta for an initial wave amplitude $\delta {b}_{0}=1$ . As can be seen, the evolution of the decay process is consistent with the predictions from Hall-MHD linear theory in the 1D cases. The pump wave is subject to decay instability, which becomes slower as the (electron) plasma beta increases. For the amplitudes considered here, the instability is suppressed at large plasma beta ( $\beta =2$ ), where the wave is more likely to decay via modulational or beat instabilities, with smaller growth rates and hence a slower decay process. The slight increase in density fluctuations that can be seen in the upper panel of Figure 2 for the 1D case does not correspond to the disruption of the pump Alfvén wave. Interestingly, the 2D simulations display an opposite trend to the plasma beta. In the low beta regime the 2D simulations are characterized by a growth rate similar to the 1D cases, although the decay occurs later than in the 1D, in agreement with previous studies comparing parametric decay from one to three dimensions (Del Zanna et al. 2001). However, the 2D simulations display a rapid decay process even in the large beta case and, contrary to the 1D case, the growth rate tends to increase with the plasma beta. We ascribe such differences between the 2D and 1D sets of simulations to the onset of filamentation/magnetosonic decay simultaneously to the main parametric decay process, and to the resulting nonlinear dynamics. Additionally, a large increase of thermal energy is always observed when the decay occurs. Although multidimensional simulations display a slightly larger final temperature than the 1D cases, this points to the fact that most of the heating mechanism(s) may be ascribed to 1D dynamics. We defer a discussion of proton heating to Section 3.3.

The right panels of Figure 2 show the same quantities in the left panel but for fixed beta ( $\beta =0.5$ ) at different initial wave amplitudes. As can be seen by inspection, the growth rate increases with the amplitude, as is expected from linear theory. Interestingly enough, a strong proton heating is observed as the pump wave amplitude increases, with parallel and perpendicular temperatures displaying the same trends for 1D and 2D simulations.

By way of illustration, the decay process for A2-2D and A3-2D is presented in Figure 3. In the left panels we plot the amplitude of the most unstable modes and of the pump wave (in black). The right panels show the superposition of the 2D fast Fourier transform (FFT) of ρ (red contours) and B_y (black contours) at the maximum of the curves in the left panels.

As can be seen, different kinds of daughter waves are excited, which leads to a competition between different types of parametric instabilities. We find that two types of decay are at play: a parallel and a quasi-parallel one, corresponding to the traditional parametric decay instability, and a perpendicular one, corresponding to the filamentation/magnetosonic instability. The parallel decay is evident in the enhancement of density fluctuations with wavenumber $({n}_{\parallel },{n}_{\perp })=(8,0)$ and, correspondingly, of the forward-propagating Alfvén wave with $({n}_{\parallel },{n}_{\perp })=(12,0)$ , in agreement with the three-wave coupling resonance condition. The quasi-parallel sideband modes are also observed for B_y and ρ with wavenumbers $(8,2)$ and $(4,2)$ (density is not shown), respectively. These quasi-parallel modes are the most unstable ones in the simulations with $\beta =2$ and $\beta =0.5$ , and the maximum amplitude of each daughter wave corresponds to about $10 \%$ of the amplitude of the pump wave. The oblique daughter wave leads to an enhancement of ρ and B_y at wavenumber $(0,2)$ . This mode is nonpropagating and it is weakly damped, its amplitude remaining constant throughout the simulation after the onset of filamentation/magnetosonic instability.

For the sake of completeness we display in Figure 4 the contour plots of the field-aligned component of the fluctuating magnetic field $\delta {b}_{x}$ (left panels) and of the pump wave b_y (right panels) at three different time for run A2-2D. The magnetic fluctuations of $\delta {b}_{x}$ are found to be highly anticorrelated with density perturbations (not shown), a signature of the slow mode character of the growing fluctuations that persist even after the posterior distortion of the pump wave. The combination of the pump wave with the daughter waves and also the dispersion generated by small-scale fluctuations leads to the steepening of the waveform (see the middle panel of Figure 4) that finally results in the disruption of the wave and the corresponding proton heating. The quasi-perpendicular mode $(0,2)$ can be easily identified in the contour of the b_x component since the amplitude of that mode remains constant after the saturation of the instability, contrary to the parallel and quasi-parallel modes that are highly damped after the saturation of the instability.

3.2. Spectral Properties

The decay of the parent Alfvén wave into secondary modes triggers nonlinear interactions that ultimately lead to the establishment of a turbulent cascade. At saturation of the instability, an energy spectrum spanning scales down to subproton scales develops preferentially in the perpendicular direction to the mean magnetic field. In Figure 5 we show the resulting magnetic and electric field energy spectra for the 2D cases shown in the left panels of Figure 2. In the top panels we plot the reduced 1D perpendicular spectrum of the parallel ( ${E}_{{B}_{\parallel }}$ , ${E}_{{E}_{\parallel }}$ ) and perpendicular ( ${E}_{{B}_{\perp }}$ , ${E}_{{E}_{\perp }}$ ) components of the electromagnetic field fluctuations at saturation stage. We also plot the reduced spectrum of density fluctuations for the same period and also the initial density spectrum as a reference to quantify the noise floor level in the simulations. The reduced 1D perpendicular spectrum is computed as ${E}_{\delta A}({k}_{\perp })\,=\int {{dk}}_{\parallel }A({k}_{\parallel },{k}_{\perp })$ , with $A({k}_{\parallel },{k}_{\perp })$ the 2D spectral energy density. The vertical dashed line marks the location of ${k}_{\perp }{\rho }_{i}=1$ , with ${\rho }_{i}=\sqrt{{\beta }_{i}}{d}_{i}$ the proton gyroradius.

**Figure 5.** (Top) The reduced 1D energy spectra as function of ${k}_{\perp }$ for ${B}_{\perp }$ (blue), ${B}_{\parallel }$ (orange), ${E}_{\perp }$ (green), and ${E}_{\parallel }$ (red) at the saturation stage for each simulation. The reduced spectrum for the initial proton density is plotted in cyan color as a reference of the particle noise level. (Bottom) Spectral ratio between ${B}_{\parallel }$ (gray) and ρ (black) with ${B}_{\perp }$ normalized with the KAW linear prediction for each simulation.
Download figure:
Standard image High-resolution image

**Figure 5.** (Top) The reduced 1D energy spectra as function of ${k}_{\perp }$ for ${B}_{\perp }$ (blue), ${B}_{\parallel }$ (orange), ${E}_{\perp }$ (green), and ${E}_{\parallel }$ (red) at the saturation stage for each simulation. The reduced spectrum for the initial proton density is plotted in cyan color as a reference of the particle noise level. (Bottom) Spectral ratio between ${B}_{\parallel }$ (gray) and ρ (black) with ${B}_{\perp }$ normalized with the KAW linear prediction for each simulation.
Download figure:
Standard image High-resolution image

The power spectrum of the magnetic field components in the parallel direction is less developed, and the spectrum is dominated by the daughter waves and its harmonics (not shown), but in the transverse direction the magnetic field shows a broad inertial range with a Kolmogorov-like spectrum ( $\sim {k}_{\perp }^{-5/3}$ ). A spectral break is also observed when approaching to proton scales, marking the transition to another turbulent regime at subproton scales. The spectral break occurs at the larger of the proton scales depending on the plasma beta, in agreement with previous numerical simulations (Franci et al. 2016). Turbulence at kinetic scales does not seem to follow a universal behavior, unlike the low-frequency, large-scale dynamics. A strong variability of the spectral index at subproton scales has indeed been reported in the solar wind and Earth's magnetosphere with values ranging between −4 and −2 (Alexandrova et al. 2009; Sahraoui et al. 2010; Chen 2016; Bowen et al. 2020). Kinetic Alfvén waves (KAWs) are often invoked to explain turbulent fluctuations at subproton scales, since magnetic field energy spectra measured by in situ observations and numerical simulation of plasma turbulence usually find power laws with a spectral index close to −7/3, in agreement with KAW theory. In this particular set of simulations, the spectrum at subproton scales is steeper than the KAW predictions, and the spectral index is about −3.5.

The change of turbulence regime from the large to the small scales is marked by the increase of plasma compressibility, an effect that we observe in all of our sets of simulations. In order to characterize small-scale fluctuations, we consider spectral field ratios that are known to provide a useful tool to investigate the polarization properties of turbulent fluctuations (Gary & Smith 2009; Chen et al. 2013; Chen & Boldyrev 2017; Grošelj et al. 2017; Cerri et al. 2019). In the bottom panels of Figure 5 we present the spectral ratios ${R}_{1}\equiv {C}_{1}{E}_{\delta {B}_{\parallel }}/{E}_{\delta {B}_{\perp }}$ and ${R}_{2}\equiv {C}_{2}{E}_{\delta \rho }/{E}_{\delta {B}_{\perp }}$ for each simulation. The ratios are normalized to the rms of each field and to the theoretical prediction from KAW at different plasma beta ( ${C}_{1}={\beta }_{t}(1+{T}_{e}/{T}_{i})/2+{\beta }_{t}(1+{T}_{e}/{T}_{i})$ and ${C}_{2}=4/((1+{T}_{i}/{T}_{e})(2+{\beta }_{t}(1+{T}_{e}/{T}_{i})))$ with ${\beta }_{t}={\beta }_{p}\,+{\beta }_{e}$ . KAW theory predicts a value of unity for R₁ and R₂ at subproton scales. This would correspond to strong compressive magnetic fluctuations in nearly pressure balance in the kinetic range. In our set of simulations, where we considered ${T}_{e}/{T}_{i}=1$ , the values of ${R}_{\mathrm{1,2}}$ for the magnetic field are in rough agreement with the theory, but the level of density fluctuations at scales smaller than proton scales largely exceeds the linear prediction. This could be due to nonlinearities and/or the presence of different wave activity like slow modes or whistler or other plasma modes in the subproton range. It is also noted that particle noise beyond proton scales may also dominate the density spectrum, and therefore the ratio overestimates the expected values.

The electric field shows a flattening of the spectrum at subproton scales with an index of about −0.8 at kinetic scales. This feature is observed in all simulations for different plasma beta. This corresponds to the increment of propagation velocity of the turbulence fluctuations at subproton scales and is well recognized to be due to the dominance of the Hall term at proton scales, which has been previously discussed in the context of fluid and kinetic simulations (Dmitruk & Matthaeus 2006; Howes et al. 2011; Franci et al. 2015).

The parallel electric field is developed during the collapsing stage of the pump wave, and it plays a crucial role in the particle heating observed during the decay process until saturation. Interestingly, an important/dominant contribution to the parallel electric field comes from the field-aligned component of the Hall term rather than from the electron pressure term. Indeed, the electron pressure presents the same trend as the Hall term, but with a smaller amplitude (not shown). In the top panel of Figure 6, we show the rms value of the parallel component of the Hall electric field ( ${E}_{{H}_{\parallel }}={E}_{{H}_{x}}$ ) and ${\rm{\Delta }}T$ for different beta simulations. As can be seen, there is a clear correlation between the Hall parallel electric field and particle heating in this set of simulations. In the bottom panel of Figure 6 we present the parallel (black lines) and perpendicular (red lines) spectral ratio between those two terms. Solid, dashed, and dotted lines refer to the same plasma beta cases reported in the top panel, with each curve being taken at the maximum of the parallel electric field and averaged over $50{{\rm{\Omega }}}_{{ci}}$ .

**Figure 6.** (Top) Temporal evolution of the rms of the parallel Hall electric field component (green lines) and the total temperature change for different plasma beta simulations. (Bottom) Spectral ratio of the parallel component of electron pressure to the parallel Hall term for the same simulations shown in the top panel.
Download figure:
Standard image High-resolution image

At large scales, the ratio between electron pressure to Hall term at transverse scales ( ${k}_{\perp }$ ) appears to increase as the plasma beta increases. This trend does not come as a surprise since ${P}_{1{\rm{D}}}({E}_{{\rm{\nabla }}{p}_{e}})/{P}_{1{\rm{D}}}({E}_{H})\propto {\beta }_{e}^{2}/4$ and it is consistent with the formation of slow modes along the perpendicular direction. In the parallel direction ( ${k}_{\parallel }$ ), instead, the Hall term is larger than the electron pressure term at all beta values, as a consequence of the presence of strong electron currents produced by the decay of the pump wave. However, at subproton scales, the electron pressure dominates over the Hall term in both parallel and perpendicular directions. Again, particle noise may contribute to overestimate the level of density fluctuations observed at small scales.

As already mentioned, the generation of parallel electric field fluctuations is crucial for the particle heating and acceleration observed during the decay process of the pump wave and until the steady-state condition is reached at saturation stage, a problem that we address in the next section.

3.3. Proton Heating

The proton dynamics is illustrated in Figure 7 for the simulation with $\beta =0.5$ and $\delta {b}_{0}=1$ . The left panels show contours of the proton distribution function in phase space (x– ${v}_{x}$ ) at three different times from top to bottom. The corresponding reduced distribution functions of parallel and perpendicular velocities (v_x and v_y) are shown in the right panels. We show the particle information at three different stages of the evolution. At $t=320{{\rm{\Omega }}}_{c}^{-1}$ , phase-space vortexes form with the same wavelength of the density fluctuations developed by the parametric instability, with not yet significant heating at that time of the evolution. After this initial stage, a "piston-like" mechanism mediated by the parallel electric field allows for the generation of a secondary proton population propagating parallel to the mean magnetic field. The beam travels at the Alfvén speed, with signatures of particle trapping clearly visible in the phase space. The proton beam is persistent and remains stable when a steady-state condition is reached.

**Figure 7.** (Left) Reduced distribution functions $f(x,{v}_{x})$ at different stages of the evolution for the run A2-2D. (Right) Reduced distribution function $f({v}_{x})$ for the x-component of the particle velocities (solid blue lines) and y-component of the particle velocities (dashed lines). The initial reduced distribution function is plotted in red.
Download figure:
Standard image High-resolution image

The averaged particle distribution function (PDF) over $100{{\rm{\Omega }}}_{{ci}}$ after the saturation stage is presented in Figure 8 for simulations with different plasma beta. We computed the average PDF right after the end of the instability is established and when the particle temperature is statistically constant. It can be noted that the beam forms around the Alfvén speed for all the beta cases. The distribution function in the perpendicular direction, instead, is a Maxwellian, and the total change of perpendicular temperature is not affected by the plasma beta.

**Figure 8.** Reduced distribution functions for (left) $f({v}_{x})$ and (right) f(v_y) averaged over the steady-state stage for runs B1-2D (blue), B2-2D (orange), B3-2D (green), B4-2D (red), and B5-2D (purple).
Download figure:
Standard image High-resolution image

The effect of the wave amplitude on the final distribution function is presented in Figure 9. Even if the heating of particles depends on the amplitude of the pump wave, the beam formation along the mean magnetic field is persistent. The core of the distribution is mostly affected by the finite amplitude effects, and larger tails in the distribution are found with larger wave amplitude. Since the distribution function is averaged over several gyroperiods, the number of particles in the beam is affected by the averaging.

The understanding of the overall proton heating due to the unstable behavior of Alfvén waves and the corresponding energy transport toward smaller scales is fundamental to understand the implications on the plasma heating typically observed in solar and astrophysical context. The nature of parallel and perpendicular heating comes from different physical mechanisms; therefore, we first focus on the parallel heating, which is in fact due mainly to the beam generation, and subsequently we discuss the perpendicular heating and the possible mechanism.

According to previous work, the electron beta plays a crucial role on the saturation process because the electron temperature contributes to the strength of the parallel electric field via the pressure gradient term in Ohm's law and also provides the coupling with the acoustic mode (see Equation 1(b)). It was suggested that the saturation of the instability was due to particle trapping by the field-aligned electric field generated by density fluctuations, which would lead to a mean field-aligned beam whose velocity appears to depend on the plasma beta. Hybrid simulations with ${\beta }_{e}=0$ have shown that the beam formation is suppressed and the saturation process results in the steepening of the ion acoustic wave, just like in the fluid description (Matteini et al. 2010). In the same spirit, we present in Figure 10 the results for simulations with ${\beta }_{p}=0.5$ and $\delta {b}_{0}=1.0$ for two different scenarios: a case with ( ${\beta }_{e}=0.5$ ) and a case with cold electrons ( ${\beta }_{e}=0$ ).

**Figure 10.** Reduced distribution functions for (left) f(v_x) and (right) f(v_y) averaged over the steady-state stage for runs D-2D (blue) and B1-2D (orange).
Download figure:
Standard image High-resolution image

As can be seen, the two simulations do not display significant differences in the PDFs. Not only is the perpendicular temperature achieved at the end of the process the same for both simulations (right panel) but, importantly, the field-aligned beam (left panel) is persistent in a plasma with ${\beta }_{e}=0$ , even though a less populated beam is observed for cold electrons. This is because the electron pressure gradient still contributes to the trapping of particles. The comparison between these two setups therefore points to the fact that it is the Hall term that contributes the most to the field-aligned electric field and to the beam formation.

According to Vlasov linear theory, proton heating by Alfvén waves is possible via resonant damping, and particles can exchange energy with the wave only at discrete resonance interaction restricted by the condition $\omega -{k}_{\parallel }{v}_{\parallel }=n{{\rm{\Omega }}}_{{ci}}$ , with ω the frequency of the wave, ${v}_{\parallel }$ the particle velocity, and ${k}_{\parallel }$ is wavenumber parallel to the magnetic field. Linear theory implies that for n = 0, when the phase speed of the wave is of the order of the proton parallel velocity, particles can resonate with the wave and then they can gain or lose parallel velocity ( ${v}_{\parallel }\leqslant \omega /{k}_{\parallel }$ or ${v}_{\parallel }\geqslant \omega /{k}_{\parallel }$ , respectively). This resonant wave–particle interaction at n = 0 represents two different physical interactions, the Landau damping, driven by a parallel electric field, and the transit time damping (TTD; Fisk 1976; Achterberg 1981), which is the magnetic analog of Landau damping. In TTD it is the mirror force ${{\boldsymbol{F}}}_{\mathrm{mir}}=\mu {{\boldsymbol{\nabla }}}_{\parallel }B$ , with $\mu ={{mv}}_{\perp }^{2}/2B$ the conserved particle magnetic moment, to play the analogous role of the parallel electric field in Landau damping and, as such, it can also contribute to the parallel heating. Interactions with $n\ne 0$ instead lead to cyclotron resonances in which the particles resonate with the oscillating electric and magnetic field (Hollweg & Isenberg 2002). This kind of wave–particle interaction results in the violation of the conservation of μ so that particles can experience strong perpendicular energization. In general these wave–particle resonances are important because they can lead to parallel/perpendicular heating produced by pitch-angle scattering resulting in an isotropization process of the particle distribution function by particle diffusion in velocity space (Kennel & Engelmann 1966; Lynn et al. 2012, 2013).

As we have shown in Section 3.1, the proton heating does not display major differences between 1D and 2D simulations whenever the wave is subject to parametric instabilities. In fact, the perpendicular heating is roughly the same, although a somewhat larger parallel temperature is obtained in the 2D simulations (see, for instance, the left panel of Figure 2). In this sense, the main proton heating mechanism at play seems to be the same process regardless of the dimensionality of the system and of the proton/electron plasma beta. In this way we can analyze the proton heating process for a 1D case without loss of generality. Moreover, the fact that in these numerical experiments the proton heating is essentially a one-dimensional mechanism, we can rule out contributions from stochastic heating possibly due to KAWs or in general to obliquely propagating modes and reconnection, which cannot develop in our 1D system.

In Figure 11, we present the spatial profiles of the parallel and perpendicular temperature and also the Hall term components ( ${E}_{{H}_{\parallel }}={b}_{z}{j}_{y}-{b}_{y}{j}_{z}$ and ${E}_{{H}_{\perp }}=\sqrt{{E}_{{H}_{y}}^{2}+{E}_{{H}_{z}}^{2}}\,={B}_{0}\sqrt{{j}_{z}^{2}+{j}_{y}^{2}}$ ) and the gradients of the magnetic field ( ${\partial }_{x}B=-{B}^{-1/2}{E}_{{H}_{\parallel }}$ ) with the current components ${j}_{y}=-{\partial }_{x}{b}_{z}$ and ${j}_{z}={\partial }_{x}{b}_{y}$ , at different stages of the evolution for run A1D-I. Initially, the circularly polarized Alfvén wave is characterized by a constant- $B$ state and zero parallel Hall electric field because of the initial condition (Figure 11(a)). Once the decay comes into play, the disruption of the wave together with dispersive effects generates gradients of $B$ in the field-aligned direction and the wave tends to steepen in some localized regions. There it follows an enhancement not only of density fluctuations but also of the Hall electric field through the current density components j_y and j_z. Such steepened wavefronts with enhanced Hall electric field propagate at the Alfvén speed (not shown). The combined effect of particle trapping by the growing electric field fluctuations, and acceleration from the localized Hall electric field at the steepened fronts of the Alfvén wave contribute to the acceleration of particles into a field-aligned beam at the Alfvén speed, enhancing the number of resonating particles with the fluctuations themselves. We suggest that both types of n = 0 resonances (Landau and TTD damping associated with the gradients of B) might be responsible for the parallel proton heating (Figures 11(b) and (c)). Once the resonating fluctuation energy is transformed into particle heating, the system achieves a steady state, with small gradients of B, a persistent beam, and a strong parallel and perpendicular heating (Figure 11(d)).

As mentioned above, perpendicular heating might be produced by pitch-angle scattering processes due to the acceleration and redistribution of particles in the field-aligned direction of phase space. There is also a high correlation between a large perpendicular temperature and a strong perpendicular Hall electric field (Figures 11(b) and (c)). However, we notice that the increment of parallel and perpendicular temperatures occurs simultaneously and at the same rate, favoring the idea that perpendicular heating is probably due to pitch-angle scattering. While here we have borrowed concepts from linear theory, the formation of the beam, and the subsequent wave–particle interactions leading to saturation and plasma heating are really a nonlinear process where finite amplitude effects on particle orbits should be taken into account. We plan to develop a more detailed investigation of the beam formation and its possible relation to both wave steepening and perpendicular heating in future work.

It is important to also mention that the interaction of protons with current sheets and fluctuations that develop naturally as a result of the turbulent cascade may account for the additional particle heating that we observed in the 2D and 3D simulations (e.g., Dmitruk et al. 2004; Zhdankin et al. 2013; Isliker et al. 2017; Pisokas et al. 2018). However, we observe a more efficient parallel than perpendicular heating (see Figure 2, second and third panels), contrary to what is expected from the strong turbulence perspective, where the effects of current sheets results in strong perpendicular proton heating.

4. Summary and Conclusion

We have performed hybrid simulations of large-amplitude, parallel-propagating Alfvén waves subject to parametric instabilities. We have investigated how their stability, nonlinear evolution, and saturation are affected by the amplitude of the pump wave, the plasma beta, and the dimensionality of the system. Our main findings can be summarized as follows:

1.
In multidimensional systems the initial decay process can be described as a superposition of both the parametric decay instability and the filamentation/magnetosonic instability. The former leads to parallel-propagating density fluctuations and Alfvénic modes, while the latter leads to the formation of perpendicular pressure-balanced fluctuations of density and magnetic field. The filamentation instability becomes the dominant process at beta values larger than unity, contrary to one-dimensional systems where the decay appears to be strongly suppressed.
2.
The decay process naturally results in a well-developed turbulent cascade preferentially in the transverse direction to the mean magnetic field, and that shows similar properties for all the plasma beta values considered. At saturation of the instability, the magnetic energy spectrum displays a Kolmogorov-like inertial range at large scales. At subproton scales, a steepened magnetic field and a flattened electric field spectrum is observed, displaying power-law scalings that appear to be consistent with observations and simulations of plasma turbulence. We have quantified the nature of the fluctuations at subproton scales by means of spectral ratio analysis, and we found a strong magnetic compressibility at subproton scales that is consistent with KAW linear theory, though a discrepancy with KAW theory is found in the excessive level of density fluctuations (possibly due to particle noise).
3.
When the decay occurs, the saturated state is always characterized by a heated plasma displaying a persistent field-aligned beam localized at the Alfvén speed. The beam also forms when electrons are cold ( ${\beta }_{e}=0$ ), pointing to the fact that it is the field-aligned Hall electric field to play the dominant role in accelerating the beam of particles. Such an electric field is enhanced at the steepened edges of the pump wave, and we argue that it mediates wave–particle interactions via both Landau and transit time damping. Landau and transit time damping are expected to become effective once the beam is formed, so that there is an increased number of resonating particles with the mean field-aligned electric field and compressible fluctuations, respectively.
4.
The overall proton heating is predominantly a one-dimensional mechanism. We argue that wave–particle resonances in the mean field-aligned directions contribute to parallel heating. The perpendicular heating is instead attributed to a stochastic mechanism that works as an isotropization process of the distribution function via pitch-angle scattering. A somewhat larger parallel temperature is observed in multidimensions than in 1D simulations, which may be explained in terms of the additional contribution due to protons interacting with turbulent fluctuations.

In conclusion, the decay process acts as a trigger to develop a turbulent cascade and to enhance wave–particle interactions, the latter resulting in a field-aligned beam and efficient plasma heating, reproducing in this way some features that are observed in the solar wind. Moreover, we have shown that the decay process remains efficient also at large values of the plasma beta ( $\beta \gt 1$ ), which makes these results relevant not only to space plasmas, but also to astrophysical environments where the plasma beta can reach values well above unity. It will be the subject of future work to investigate further the (nonlinear) wave–particle interactions leading to the beam formation and to the observed plasma heating to corroborate the scenario proposed here. It will also be important to extend our results to the expanding solar wind, and to investigate the evolution of a spectrum of Alfvén waves within the Expanding Box model by including a population of alpha particles (Maneva et al. 2015) to assess whether the instability is favored by the expansion or not, what is the contribution of the decay process to solar wind heating and beam acceleration, and how minor ions react/modify the overall evolution.

This research was supported by NASA grant #80NSSC18K1211. We also acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper. URL: http://www.tacc.utexas.edu.

The Role of Parametric Instabilities in Turbulence Generation and Proton Heating: Hybrid Simulations of Parallel-propagating Alfvén Waves

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Model and Simulation Setup