the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
On heating of solar wind protons by the parametric decay of largeamplitude Alfvén waves
Horia Comişel
Yasuhiro Nariyuki
Yasuhito Narita
Uwe Motschmann
By threedimensional hybrid simulations, proton heating is investigated starting from a monochromatic largeamplitude Alfvén wave with lefthanded circular polarization launched along the mean magnetic field in a lowbeta plasma. We find that the perpendicular scattering is efficient in three dimensions and the protons are heated by the obliquely propagating waves. The thermal core proton population is heated in three dimensions as well in the longitudinal and parallel directions by the fieldaligned and obliquely propagating sound waves out of the parametric decay. The astrophysical context is discussed.
Early in situ measurements at 1 AU from the VELA satellite (Bame et al., 1975) reveal that the velocity distribution function of solar wind protons is broader in the direction perpendicular to the mean magnetic field (hereafter the z direction) than in the parallel direction. This velocity anisotropy indicates a higher perpendicular temperature than the parallel one. Marsch et al. (1982) have found by using Helios 1 and Helios 2 data that such an anisotropic plasma heating occurs in high speed solar wind streams from 0.3 to 1.0 AU. This problem of anisotropic heating of ions in solar wind and solar corona is vast and has been discussed for a long time in space plasma physics (Ofman, 2010). Theoretical models (Tu and Marsch, 1995) based on cyclotron resonant or nonresonant processes have been proposed to explain the anisotropic heating of solar corona and solar wind.
Marsch and Tu (2001) have shown for the first time the observational evidence for the occurrence of the pitchangle scattering of solar wind protons, driven by resonance with ion cyclotron waves propagating away from the Sun. The perpendicular broadening of the sunward part of the measured distributions has been explained through the pitchangle scattering of solar wind protons resonantly interacting with the outward parallelpropagating Alfvén waves. In a later paper, Marsch and Bourouaine (2011) have shown that the antisunward part of the proton distribution functions can be similarly shaped by the proton diffusion by the oblique fast magnetosonic and Alfvén waves propagating away from the Sun. According to numerical simulation studies (Araneda et al., 2007), the fieldaligned part describing the tail or the proton beam of the velocity distribution functions can originate in the parametric decay of the Alfvén waves, a process predicted by theories (Sagdeev and Galeev, 1969) and supported by observations (Spangler et al., 1997).
Parametric instabilities play an important role in the dissipation of the largeamplitude Alfvén waves with parallel or quasiparallel propagation with respect to the mean magnetic field and in plasma heating by means of the ion Landau damping mechanism. Parametric instabilities, including decay, modulational, and beat instabilities, have been extensively analyzed by theoretical studies (Viñas and Goldstein, 1991a, b) or numerical magnetohydrodynamics (MHD) (Del Zanna et al., 2001; Ghosh and Goldstein, 1994; Ghosh et al., 1993, 1994) and particleincell or hybrid (Gao et al., 2013; Matteini et al., 2010a, b; Nariyuki et al., 2012, 2014; Terasawa et al., 1986; Verscharen et al., 2012) simulations. In the MHD picture, the plasma heating by the Alfvén wave can occur through generation and steepening of ion acoustic waves. A shock wave is formed as a result of the wave steepening at a late (and nonlinear) saturation stage of the parametric decay. In the kinetic picture, hybrid simulations prove that the heating mechanism is completed by kinetic effects and a beam can be created in the ion distribution function due to the nonlinear trapping of protons (Araneda et al., 2008; Matteini et al., 2010a, b). The velocity beam formation is however restricted by the conditions of lowbeta plasmas (Matteini et al., 2010a).
Here we address the question “Is the stochastic ion heating stronger in a 3D parametric decay?” Our question is motivated by two preceding studies. First, Ghosh et al. (1993, 1994) and Ghosh and Goldstein (1994) discovered from the 2D numerical MHD study that a parallelpropagating Alfvén wave collapses into obliquely propagating daughter waves by the parametric decay. Second, more recently, Gao et al. (2013) confirm in the 2D hybrid simulation that obliquely propagating Alfvén waves are indeed excited by the fieldaligned parametric decay, and propose a heating mechanism of the ambient plasma in a stochastic fashion. When the daughter Alfvén wave propagates obliquely to the mean magnetic field, the types of particle trajectories can be more diverse (see the illustration in Fig. 1). The finding and the assumed mechanism above by Ghosh et al. (1993, 1994), Ghosh and Goldstein (1994), and Gao et al. (2013) are still limited to a 2D numerical setup. Obliquely propagating waves are limited to a plane spanning parallel and perpendicular to the mean magnetic field in the 2D setup, whereas the wavevectors can have a higher degree of freedom in the azimuthal directions around the mean magnetic field.
We perform a 3D hybrid plasma simulation for the parametric decay, and track the time evolution of the proton distribution functions. We find that the stochastic heating (i.e., pitchangle scattering) occurs more quickly and the ions are heated most strongly in the 3D treatment. Our finding that the particles can be more quickly heated by the 3D parametric decay can be tested by in situ measurements by the upcoming heliospheric missions such as Parker Solar Probe (Fox et al., 2016) and Solar Orbiter (Müller et al., 2013).
We perform hybrid simulations with the AIKEF hybrid code (Müller et al., 2011), conducted in a threedimensional configuration: the size of the simulation box in each direction is L=288 d_{i}, the grid size is 1d_{i}, and 1000 superparticles are used for each computational cell. Here, ${d}_{i}={V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}$ is the ion inertial length, while V_{A} and Ω_{p} are the Alfvén velocity and ion frequency, respectively.
The parametric decay modeled in the actual study is a threewave process starting from a largeamplitude monochromatic Alfvén pump wave propagating parallel to the mean magnetic field B_{0}, a spectrum of electrostatic ion acoustic waves also at parallel propagation, and a spectrum of Alfvén daughter waves at antiparallel propagation. The amplitude of the pump wave with lefthanded circular polarization is normalized to the value of the ambient magnetic field and has a value of 0.2, while its wavenumber and the resonant frequency are k_{0}≈0.218 V_{A}∕Ω_{p} and ω_{0}≈0.19 Ω_{p}, respectively. The initial fluctuating magnetic field (B_{⟂}) and bulk velocity (u_{⟂}) satisfy the relation ${\mathit{u}}_{\u27c2}={k}_{\mathrm{0}}/{\mathit{\omega}}_{\mathrm{0}}{\mathit{B}}_{\u27c2}$. The resonant frequency ω_{0} is determined from the dispersion relation ${k}_{\mathrm{0}}^{\mathrm{2}}={\mathit{\omega}}_{\mathrm{0}}^{\mathrm{2}}/(\mathrm{1}{\mathit{\omega}}_{\mathrm{0}})$ for the lefthanded waves (Terasawa et al., 1986). The seed amplitudes for the daughter sound waves are implicitly included by the simulation noise. A low value of β=0.01 is used for the plasma beta parameter for each species of particles which is relevant for the solar corona and inner heliosphere studies.
The protons are treated in the hybrid scheme as particles, while the electrons are considered as a massless fluid. The values of the β parameter and of the pump wavenumber k_{0} are selected such that the decay instability has growth rates larger than other parametric instabilities (e.g., beat instability and modulational instability expected for lefthanded polarized waves) and can be safely evaluated by MHD theories. By using the analytic study discussed by Terasawa et al. (1986) in the twofluid description of plasma, the growth rate of the decay instability has a maximum value of γ_{th}=0.0358 corresponding to a compressional wave excited at wavenumber $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}=\mathrm{0.385}$. More weaker, the beat and modulational instabilities are estimated at wavenumbers $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}\approx \mathrm{0.218}$ and $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}\approx \mathrm{0.075}$, respectively.
The simulated magnetic field, density, and bulk velocity fluctuations are first averaged in the real space over one of the perpendicular directions (x direction). After the averaging we obtain a 2D representation of wavevectors with a parallel (${k}_{\parallel}\equiv {k}_{z}$) and a perpendicular (${k}_{\u27c2}\equiv {k}_{y}$) component. The power spectrum is constructed in the $({k}_{\parallel},{k}_{\u27c2})$ domain by Fourier analysis of the 2D spatially averaged fluctuations. In our setup, the Alfvén pump wave and the fieldaligned compressional – and Alfvén – daughter waves have the Fourier modes (${m}_{\parallel},{m}_{\u27c2}$) with values of (10,0), (18,0), and (−8,0), respectively (${m}_{\parallel (\u27c2)}={k}_{\parallel (\u27c2)}L/\mathrm{2}\mathit{\pi}$), according to the wave–wave coupling rules. The negative sign expresses the backward (antiparallel) propagation of the Alfvén daughter wave. In our astrophysical scenario, the Alfvén pump wave is propagating away from the Sun.
Figure 2a shows the time evolution of the antisunwardpropagating Alfvén pump (10,0), the compressional daughter (18,0), and the sunwardpropagating Alfvén daughter (−8,0) modes excited by the fieldaligned parametric decay. The linear growing of the daughter modes terminates close before time tΩ_{p}=300. At later times, the increase is very slow, while the pump wave remains stronger until the end of the simulation time. Besides the fieldaligned decay, moderateoblique propagating daughter waves are excited. Here we show, plotted by dashed and dotted lines, the oblique modes with perpendicular Fourier numbers ${m}_{\u27c2}=\mathrm{8}$ and ${m}_{\u27c2}=\mathrm{12}$. The growth rates are computed as slopes of the density modes, $\mathit{\gamma}=\mathrm{1}/\left\mathit{\delta}\mathit{\rho}\right\mathrm{d}/\mathrm{d}t\left\mathit{\delta}\mathit{\rho}\right$. The fieldaligned growth rate γ_{∥} has a value of 0.036 close to the above estimated γ_{th}. The growth rates for the moderateoblique daughter waves have close values to the parallel daughter mode. Viñas and Goldstein (1991b) show that there is a trend to eliminate/reduce the differences between the growth rates of the oblique and fieldaligned decay instabilities while plasma beta values are decreasing to lower values. Following Viñas and Goldstein (1991a, b) analytical treatment, our calculations indeed reveal that the oblique growth rates (at propagation angles of 10 and 20^{∘}) have close values to the fieldaligned growth rate for β=0.01, thus proving the predicted tendency of merging at lowbeta values.
The rootmeansquared (rms) density fluctuations δρ and the crosshelicities σ_{c} are represented as a function of time in Fig. 2 (middle panels). The normalized crosshelicity is defined as ${\mathit{\sigma}}_{c}=(<\mathit{\delta}\mathit{B}\cdot \mathit{\delta}\mathit{u}>)/(<\mathit{\delta}\mathit{B}{>}^{\mathrm{2}}+<\mathit{\delta}\mathit{u}{>}^{\mathrm{2}})$, where $<\mathit{\delta}\mathit{B}>$ are the averaged magnetic field fluctuations normalized to the mean field B_{0}, while $<\mathit{\delta}\mathit{u}>$ are the averaged bulk velocity fluctuations normalized to the Alfvén velocity v_{A}, and has in our study a value of +1 for the antisunwardpropagating Alfvén pump wave.
When the decay instability is set on, the density fluctuations start to increase and the crosshelicity starts to decrease. At a time just before $t=\mathrm{300}{\mathrm{\Omega}}_{p}^{\mathrm{1}}$, the compressional fluctuations reach a maximum value associated with the saturation of the instability while the sharp decreasing of the crosshelicity terminates within a narrow plateau. A weaker attenuation of the crosshelicity continues down to the value of −0.5 at the latest time of simulation. In earlier studies concerning the parametric decay (Del Zanna et al., 2001), the crosshelicity of the waves tends to change from positive to negative crosshelicity values for lowbeta simulations. As one can see in Fig. 2, at time tΩ_{p}≈300 when σ_{c}=0, the pump wave is still dominant over the parallelpropagating daughter mode, thus suggesting that a broadband spectrum of obliquely propagating waves is developed.
The power spectrum of magnetic field and density fluctuations is given in Fig. 2b at time tΩ_{p}=200 during the linear growth of the daughter modes. The spectrum of the decomposed magnetic field fluctuations shows lefthanded antisunwardpropagating waves (or righthanded sunwardpropagating waves) drawn by black solid line and lefthanded sunwardpropagating waves (or righthanded antisunwardpropagating waves) given by gray solid line. The spectrum of the incompressible sense of magnetic energy ($\mathit{\delta}{B}_{x}^{\mathrm{2}}+\mathit{\delta}{B}_{y}^{\mathrm{2}}$) is dominated by the pump wave at $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}=\mathrm{0.218}$. The first broadband gray peak and the third black peak (counting from the lefthand side of the plot) are localized around the lower and upper wavenumbers of the sideband daughter modes (${m}_{\parallel}=\mathrm{10}\pm \mathrm{18}$ corresponding to kV_{A}∕Ω_{p}=0.61 and $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}=\mathrm{0.17}$, respectively). The compressional daughter mode is the strongest peak in the density spectrum at the wavenumber close to the predicted value from our analytical study ($k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}=\mathrm{0.385}$). Its second harmonics can also be identified as the second following peak. At smaller wavenumbers close to the pump wavenumber ($k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}=\mathrm{0.218}$), additional broadband compressional modes are accompanying the ion acoustic daughter wave. Most of the peaks observed in the magnetic field spectrum correspond to the wave–wave couplings of the pump wave with the fundamentals and the harmonics of the daughter compressional wave driven by the decay instability. Exceptions are the peaks localized at $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}\approx \mathrm{0.4}$ which can be related to the beat instability and correspond to the coupling of the pump wave with the compressional mode excited at the same wavenumber. According to MHD studies (Bekhor and Drake, 2003), the beat instability driving a fast magnetosonic mode is efficient at larger ion beta plasmas and moderateoblique propagation angles.
While Fig. 2 based on analyzing the wave spectrum brings evidence of the decay instability, Fig. 3 reports the particle heating process. The upper panels of Fig. 3 present the particle distribution functions in the phase space $z{v}_{\parallel}$ at two different stages for the evolution of decay instability. The time evolution of the velocity distribution functions is usually helpful to emphasize the role of the kinetic regime in the saturation of the instability via particle trapping and wave particle interactions. The left panel of Fig. 3 refers to the early stage just after the saturation of the instability at a time of tΩ_{p}=300, while the right panel corresponds to a later time of tΩ_{p}=600. The proton phase space $z{v}_{\parallel}$ shown in the upper panels of Fig. 3 is similar to former studies. Matteini et al. (2010a) explain the spatial modulation and the modulation in enhancement of the parallel electric field (Fig. 5 in their paper) due to the broader spectrum of ion acoustic waves excited by the largeamplitude Alfvén mother wave. The spatial modulation in Fig. 3 is weaker according to the smalleramplitude pump wave used in our simulation. During the saturation of the instability we do not observe the presence of phasespace vortices leading to the formation of a proton beam. This is a consequence of the low values for the electron beta (β=0.01) and ion beta (β=0.01). At very low electron temperatures, the contribution of the electron pressure term to the electric field ($\mathrm{\nabla}{P}_{\mathrm{e}}=\mathrm{\nabla}n{k}_{{\mathrm{BT}}_{\mathrm{e}}}$) is small and the particle density fluctuations are less efficient in coupling to the electric field fluctuations. In a first stage, the protons are accelerated by the parallel electric field produced by the density fluctuations. In the later stage (tΩ_{p}=600), the particles are smoothly heated and the resonant protons accelerated by the ion acoustic waves are mixed with the thermal core of the distribution.
The lower panels of Fig. 3 report the protonreduced velocity distribution functions constructed in the (${v}_{\parallel},{v}_{\u27c2}$) plane obtained at the initial time and at the same simulation epochs as in the upper panels. The reduced velocity distribution functions ${f}_{r}({v}_{\parallel},\pm {v}_{\u27c2})$ are computed by counting the number of particles d$N={v}_{\u27c2}\mathrm{d}{v}_{\u27c2}\mathrm{d}{v}_{\parallel}\int f({v}_{\u27c2},{v}_{\parallel},\mathrm{\Phi})\mathrm{d}\mathrm{\Phi}$. Here, ${v}_{\u27c2}=\mathrm{sgn}(\mathit{\pi}\mathrm{\Phi})\sqrt{{v}_{x}^{\mathrm{2}}+{v}_{y}^{\mathrm{2}}}$ is the velocity component perpendicular to the mean magnetic field, ${v}_{\parallel}\equiv {v}_{z}$ is the parallel velocity, and the integral over the azimuthal angle $\mathrm{\Phi}=\mathrm{arctan}{v}_{x}/{v}_{y}$ is done within the interval [0,2π]. The contour lines are given for fractions of 60 %, 40 %, 30 %, 20 %, 10 %, 5 %, and 1 % of the maximum phase space density from the inner to outer parts of the distribution functions.
Due to the transversal wave field imposed as the initial condition, the velocity distribution functions at time tΩ_{p}=0 are shifted towards the initial bulk velocities. The deformation of the distribution functions with respect to the Maxwellian shape due to the presence of a wave field of forces can drive an apparent temperature anisotropy; see, e.g., Verscharen and Marsch (2011). The wave effect on the velocity distribution function can be described by a model distribution function which is equivalent to a Maxwellian distribution shifted by the mean fluid velocity of the particles associated with the local magnetic field (Nariyuki, 2011; Verscharen and Marsch, 2011). In many numerical simulations, e.g., Markovskii et al. (2009), such a modified Maxwellian distribution is used at the initial setup to properly account for the motion of particles in the wave field. The wave field effect on the distribution function by applying the pump wave is however much less important for the distribution functions evaluated at later times because the ratio of the bulk and thermal energies describing the velocity shift of the Maxwellian distribution becomes dominant by considering the wave decay and oscillation of the daughter waves. The symmetrical sets of contour levels with respect to the ${v}_{\u27c2}=\mathrm{0}$ axis are slowly merging with the time evolution of the velocity distributions.
The dashed lines at time tΩ_{p}=600 representing segments of circles describe the diffusion plateaus for the resonant protons with positive velocities scattered by backwardpropagating daughter waves (Isenberg and Lee, 1996),
where v_{0∥} is the initial value of the parallel proton velocity v_{∥} satisfying the following resonance condition:
The above equations are derived within the quasilinear theory of wave particle interaction in magnetized turbulent plasma (Kennel and Engelmann, 1966). The segment of the circle represented by the rightmost dashed line in Fig. 3 (bottomright panel) defines the cyclotron diffusion plateau obtained by numerically solving the integral in Eq. (1) including Eq. (2) and the cold plasma dispersion relation modeled by
The resulting level plateau has the center localized at a value of ${V}_{\mathrm{center}}\approx {V}_{\mathrm{A}}/\mathrm{2}$ along the negative axis of the parallel proton velocities. The inner dashed segments are obtained by slightly shifting V_{center} towards the ${v}_{\parallel}=\mathrm{0}$ axis. Similar diffusion plateaus can be derived for the lefthand side of the plot (with sunward velocity component) corresponding to the resonant scattering of protons by antisunwardpropagating waves.
It is straightforward to notice by comparing the distribution functions at the intermediate (tΩ_{p}=300) and final (tΩ_{p}=600) times that the contour levels are moderately enlarging in both the parallel and perpendicular directions following the diffusion plateaus of energy conservation driven by the pitchangle scattering of protons. Observational evidences of diffusion plateaus formed by solar wind protons have been reported starting with the paper of Marsch and Tu (2001) followed by later observational (He et al., 2015; Heuer and Marsch, 2007; Marsch and Bourouaine, 2011; Tu and Marsch, 2002) or particleincell simulation (Gary and Saito, 2003) studies. Marsch and Bourouaine (2011) have found that the oblique propagation of waves is the key factor which enables solar wind protons to scatter along the plateau levels of the velocity distribution functions.
Here we will check whether fieldaligned propagating waves can resonate and whether scatter protons by pitchangle diffusion mechanisms or obliquely propagating waves are necessary to explain the plateau levels observed in Fig. 3. The resonance velocity defined in Eq. (2) is plotted by a solid line along the wavenumber axis in Fig. 4. The right vertical axis shows the values of the resonance velocity in terms of Alfvén velocity V_{A}. The frequency–wavenumber (ω−k) spectrum corresponding to the sunwardpropagating Alfvén daughter waves is determined at time tΩ_{p}=600 and represented in a graycoded scale in the same figure. By dashed line is represented the cold plasma dispersion relation. Besides the Alfvén daughter wave observed at wavenumber $k{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}\approx $ 0.18 as the strongest mode driven by the decay instability, a trend of additional normal modes is excited and forms a relatively broadband (turbulent) spectrum at larger wavenumbers and frequencies. The maximum value of the parallel velocity ${v}_{\parallel \mathrm{max}}\approx \mathrm{0.4}{V}_{\mathrm{A}}$ reported in the velocity distribution functions for the outmost level contour corresponds to a resonance velocity at the parallel wavenumber ${k}_{\parallel}{V}_{\mathrm{A}}/{\mathrm{\Omega}}_{p}\approx \mathrm{1}$. Due to the lack of waves in this wavenumber region, the observed parallelpropagating daughter waves cannot explain the formation of plateau levels in the velocity distribution functions. The remaining candidates capable of scattering the particles along the segments of circles at moderate parallel velocities are the obliquely propagating modes. At oblique inclination angles, the dispersion relation branch becomes less tilted with respect to the wavenumber axis while the resonance velocity is shifted towards smaller (absolute) wavenumbers; therefore the resonance condition is expected to be fulfilled. However, an extensive analysis of the pitchangle scattering of protons by obliquely propagating waves in the framework of quasilinear theory is beyond the objectives of the actual study, while the determination of the frequency–wavenumber spectrum of magnetic field fluctuations at oblique angles is hard to achieve due to the requirement of sufficient spatial resolution.
The important role of the obliquely propagating daughter waves in perpendicularly heating the particles can be alternatively emphasized by the comparison with the results obtained from additional simulations carried out by decreasing the dimensionality from the 3D down to the 2D and 1D configurations, while all the other physical and numerical parameters are basically maintained the same. In the 1D box spatial variations are allowed only in directions parallel and antiparallel to the background magnetic field, whereas in the 2D simulation spatial variations are allowed in both the parallel/antiparallel direction and one perpendicular direction. Figure 5 shows the time evolution of the parallel and perpendicular proton temperatures obtained in the actual and additional 2D and 1D setups. The temperatures are determined by computing the thermal velocities (as the secondorder velocity moment) obtained by subtracting the bulk velocity from the full particle velocities according to the definition of kinetic temperature. The time evolution of the parallel temperature of protons is similar for all the simulation runs. In contrast, the perpendicular temperature for the 3D setup starts to increase and becomes larger by a factor of 4 with respect to the temperatures obtained in the downgraded configurations. The three simulations carried out with the same physics in one dimension, two dimensions, and three dimensions demonstrate that the 3D simulation yields the strongest proton heating.
With respect to this result we have done a quantification of the differences between the results of the 2D and 3D setups in what the amplitude and slope of the oblique modes mean. At the linear stage of the instability growth, the Alfvén oblique mode (−8,8) shown in Fig. 2 has close amplitude and slope as the fieldaligned daughter mode (−8,0). The correspondent mode in the 2D setup is about 3 times weaker in amplitude than the fieldaligned daughter mode computed in this configuration. The perpendicular projection of the obliquely propagating density mode (18,8) is more unstable and grows faster than its counterpart from the 2D setup. On the other hand, the 3D case opens new channels of the parametric instabilities in the sense that there is a larger population of obliquely propagating waves due to a degree of freedom in the azimuthal direction around the mean magnetic field. As the oblique propagating waves play an important role in the heating process of lowbeta plasmas, the enhanced oblique daughter waves lead to a more efficient ion heating in the perpendicular direction.
Viñas and Goldstein (1991a, b) discovered new channels of parametric instabilities when the dimension of the analyzed MHD system is increased from 1D to 2D, thus allowing the daughter and sideband modes to obliquely propagate with respect to the fieldaligned Alfvén pump wave. Some of these new instabilities, e.g., the filamentation instability, have been identified in later MHD simulations (Ghosh et al., 1993, 1994) by comparison with the analytical predictions of Viñas and Goldstein (1991b). To emphasize new instabilities particularly enhanced in 3D with respect to 2D setups, an analytical study is needed by considering both the oblique and azimuthal propagation angles of daughter modes. However, analytic treatments concerning oblique instabilities developed by parallelpropagating Alfvén pump waves have not continued since the Viñas and Goldstein (1991a, b) papers. An outcome of these former theoretical studies confirmed in the present numerical study is that the obliquely propagating daughter waves are stronger and play a more important role in the dynamics of the heating process in lowbeta plasmas than the 1D fieldaligned pump case.
The used spatial resolution for the field quantities (magnetic field, electric field, and velocity moments) is close to the ion inertial length and the proton gyroradius (ρ_{i}∼0.1d_{i}) or smaller spatial gradients cannot be resolved. The magnetic field within a numerical cell is overall homogeneous, with the linear interpolations at the particle position between mesh points or due to the wave magnetic field. Thus, the perpendicular projection of the proton motion is nearly a circle and this circular gyration is resolved by about 100 time steps. Gradients become important over about 10 gyroradii and not just over one gyration. We are warned that numerical heating could have some contribution in our simulations. Among various candidate mechanisms causing numerical heating one may specify the numerical noise given by the statistical representation of the distribution functions, the rounding error or cutoff error when evaluating the differential operator, the absorption of the numerically arising electric (possibly the electrostatic field) by the ions, and the random scattering due to the numerically fluctuating magnetic field (here the magnetic diffusion may be applicable). The numerical free energy occurring in the system can be converted into wave energy. This wave energy can be absorbed by particles and heating of the plasma. The heating effects described above can be compensated by using a suitable resistivity parameter, a smoothing procedure for the magnetic field, and numerical tests including various parameters. We have tested simulation runs with or without using a pump wave by varying the number of particles per cell, time steps δt, and grid sizes to find out sufficient energy accuracy (within 5 % for 500 elapsed ion gyroperiods). Thus we conclude that the numerical heating does not play a significant role compared to the physical heating.
We performed for the first time a 3D hybrid simulation study on the plasma heating problem associated with the parametric decay. The analysis of wave–wave coupling driven by decay instability and the time evolution of the main wave modes and crosshelicity bring evidence that obliquely propagating waves are excited early by the fieldaligned Alfvén pump wave.
We draw the following conclusions.

The pitchangle scattering is efficient in three dimensions and the plasma becomes heated more quickly by the obliquely propagating daughter waves. The stochastic heating can be verified by the upcoming solar wind measurements by Parker Solar Probe and Solar Orbiter. A temperature rise by a factor of about 4 is obtained in our 3D hybrid simulation study by a time of 300 to 600 ion gyroperiods. When applying the mapping of the elapsed time to the radial distance from the Sun advected by the solar wind (Comişel et al., 2015), we obtain a radial distance of about 0.1 to 0.2 AU from the Sun at a time of 300 to 600 gyroperiods.

Thermal core particle population is effectively heated in three dimensions as well in the longitudinal (to the wavevector) and parallel (to the mean magnetic field) directions by the fieldaligned and the obliquely propagating sound waves out of the parametric decay, confirming the lessons from the earlier studies (Gao et al., 2013; Laveder et al., 2002).
Needless to say, our conclusions are limited to a beta parameter of 0.01. Threedimensional hybrid simulations provide more realistic predictions for the waveheating problem in nonlinear space plasma dynamics. We propose to study the following items to extend the 3D simulations. Viñas and Goldstein (1991a, b), Ghosh et al. (1993, 1994), and Ghosh and Goldstein (1994) discovered the importance of the oblique waves in the parametric decay in dependence of the beta parameter. One may expect that plasma beta parameter β could be the key factor in controlling the perpendicular ion heating driven by the parametric decay instability in space plasmas.
Data from our hybrid simulations are stored at the Institut fuer Theoretische Physik – Technische Universitaet Braunschweig. Data can be obtained by writing to the following email addresses: h.comisel@tubraunschweig.de or comisel@spacescience.ro.
HC carried out the simulations, calculation, and writing. Yasuhir. Nariyuki came out with analysis and opinion. Yasuhit. Narita contributed with writing and conclusive remarks. UM handled the discussion and finalization.
The authors declare that they have no conflict of interest.
This work is financially supported by a grant of the Deutsche
Forschungsgemainschaft (DFG grant MO539/201). HC acknowledges the
hospitality at the University of Toyama for hosting the research visit. We
acknowledge the John von Neumann Institute for Computing (NIC) for providing
computing time on supercomputer JURECA at Juelich Supercomputer Centre (JSC).
Edited by: Yoshizumi Miyoshi
Reviewed by: three anonymous referees
Araneda, J. A., Marsch, E., and Viñas, A. F.: Collisionless damping of parametrically unstable Alfvén waves, J. Geophys. Res., 112, A04104, https://doi.org/10.1029/2006JA011999, 2007. a
Araneda, J. A., Marsch, E., and Viñas, A. F.: Proton Core Heating and Beam Formation via Parametrically Unstable AlfvénCyclotron Wave, Phys. Rev. Lett., 100, 125003, https://doi.org/10.1103/PhysRevLett.100.125003, 2008. a
Bame, S. J., Asbridge, J. R., Feldman, W. C., Montgomery, M. D., and Gary, S. P.: Evidence for local ion heating in solar wind high speed streams, Geophys. Res. Lett., 2, 373–375, https://doi.org/10.1029/GL002i009p00373, 1975. a
Bekhor, S. H. and Drake, R. P.: Plasma heating via parametric beating of Alfvén waves, with heliospheric applications, Phys. Plasmas, 10, 4800, https://doi.org/10.1063/1.1619975, 2003. a
Comişel, H., Motschmann, U., Buechner, J., Narita, Y., and Nariyuki, Y.: Ionscale turbulence in the inner heliosphere: radial dependence, Astrophys. J., 812, 175 pp., https://doi.org/10.1088/0004637X/812/2/175, 2015. a
Del Zanna, L., Velli, M., and Londrillo, P.: Parametric decay of circularly polarized Alfvén waves: Multidimensional simulations in periodic and open domains, Astron. Astrophys., 367, 705–718, https://doi.org/10.1051/00046361:20000455, 2001. a, b
Fox, N. J., Velli, M. C., Bale, S. D., Decker, R., Driesman, A., Howard, R. A., Kasper, J. C., Kinnison, J., Kusterer, M., Lario, M., Lockwood, M. K., McComas, D. J., Raouafi, N. E., and Szabo, A.: The solar probe plus mission: Humanity's first visit to our star, Space Sci. Rev., 204, 7–48, https://doi.org/10.1007/s1121401502116, 2016. a
Gao, X., Lu, Q., Li, X., Shan, L., and Wang, S.: Parametric instability of a monochromatic Alfven wave: Perpendicular decay in low beta plasma, Phys. Plasmas, 20, 072902, https://doi.org/10.1063/1.4816703, 2013. a, b, c, d
Gary, S. P. and Saito, S.: Paricleincell simualtions of Alfvéncyclotron wave scattering: Proton velocity distributions, J. Geophys. Res., 108, 1194, https://doi.org/10.1029/2002JA009824, 2003. a
Ghosh, S. and Goldstein, M. L.: Nonlinear evolution of a largeamplitude circularly polarized Alfvén wave: low beta, J. Geophys. Res., 99, 13351–13362, https://doi.org/10.1029/94JA00474, 1994. a, b, c, d
Ghosh, S., Viñas, A. F., and Goldstein, M. L.: Parametric instabilities of a largeamplitude circularly polarized Alfvén wave: Linear growth in twodimensional geometries, J. Geophys. Res., 98, 15561–15570, https://doi.org/10.1029/93JA01534, 1993. a, b, c, d
Ghosh, S., Viñas, A. F., and Goldstein, M. L.: Nonlinear evolution of a largeamplitude circularly polarized Alfvén wave: high beta, J. Geophys. Res., 99, 19289–19300, https://doi.org/10.1029/94JA00474, 1994. a, b, c, d
He, J., Wang, L., Ru, C., Marsch, E., and Zong, Q.: Evidence of Landau and cyclotron resonance between protons and kinetic waves in solar wind turbulence, Astrophys. J. Lett., 800, L31, https://doi.org/10.1088/201418205/800/2/L31, 2015. a
Heuer, M. and Marsch, E.: Diffusion plateaus in the velocity distributions of fast solar wind protons, J. Geophys. Res., 112, A03102, https://doi.org/10.1029/2006JA011979, 2007. a
Hoshino, H. and Goldstein, M. L.: Time evolution from linear to nonlinear stages in magnetohydrodynamic parametric instabilities, Phys. Fluids BPlasma, 1, 1405, https://doi.org/10.1063/1.858971, 1989.
Isenberg, P. and Lee, M.: A dispersive analysis of biospherical pickup ion distributions, J. Geophys. Res., 101, 11055–11066, https://doi.org/10.1029/96JA00293, 1996. a
Kennel, C. and Engelmann, F.: Velocity space diffusion from weak plasma turbulence in a magnetic field, Phys. Fluids, 9, 2377–2388, 1966. a
Laveder, D., Passot, T., and Sulem, P. L.: Transverse dynamics of dispersive Alfvén waves. II. Driving of a reduced magnetohydrodymanic flow, Phys. Plasmas, 9, 305–314, https://doi.org/10.1063/1.1417511, 2002. a
Markovskii, S. A., Vasquez, B. J., and Hollweg, V.: Proton heating by nonlinear fieldaligned Alfvén waves in solar coronal holes, Astrophys. J., 695, 1413–1420, https://doi.org/10.1088/0004637X/695/2/1413, 2009. a
Marsch, E. and Tu, C.Y.: Evidence for pitch angle diffusion of solar wind protons in resonance with cyclotron waves, J. Geophys. Res., 106, 8357, https://doi.org/10.1029/2000JA000414, 2001. a, b
Marsch, E., Mühlhäuser, K.H., Schwenn, R., Rosenbauer, H., Pilipp, W., and Neubauer, F. M.: Solar wind protons: Threedimensional velocity distributions and derived plasma parameters measured between 0.3 and 1 AU, J. Geophys. Res., 87, 52–72, https://doi.org/10.1029/JA087iA01p00052, 1982. a
Marsch, E. and Bourouaine, S.: Velocityspace diffusion of solar wind protons in oblique waves and weak turbulence, Ann. Geophys., 29, 2089–2099, https://doi.org/10.5194/angeo2920892011, 2011. a, b, c
Matteini, L., Landi, S., Del Zanna, L., Velli, M., and Hellinger, P.: Kinetics of parametric instabilities of Alfvén waves: Evolution of ion distribution functions, J. Geophys. Res., 115, A09106, https://doi.org/10.1029/2009JA014987, 2010a. a, b, c, d
Matteini, L., Landi, S., Del Zanna, L., Velli, M., and Hellinger, P.: Parametric decay of linearly polarized shear Alfvén waves in oblique propagation: One and twodimensional hybrid simulations, Geophys. Res. Lett., 37, L20101, https://doi.org/10.1029/2010GL044806, 2010b. a, b
Müller, D., Marsden, R. G., Cyr, O. C. St., and Gilbert, H. R.: Solar Orbiter Exploring the sunheliosphere connection, Sol. Phys., 285, 25–70, https://doi.org/10.1007/s1120701200857, 2013. a
Müller, J., Simon, S., Motschmann, U., Schüle, J., Glassmeier, K.H., and Pringle, G. J.: A.I.K.E.F.: Adaptive hybrid model for space plasma simulations, Comp. Phys. Comm., 182, 946–966, https://doi.org/10.1016/j.cpc.2010.12.033, 2011. a
Nariyuki, Y.: On entropymaximized velocity distributions in circularly polarized finite amplitude Alfvén waves, Phys. Plasmas, 18, 052112, https://doi.org/10.1063/1.3590857, 2011. a
Nariyuki, Y., Hada, T., and Tsubouchi, K.: Nonlinear dissipation of circularly polarized Alfvén waves due to the beaminduced obliquely propagating waves, Phys. Plasmas, 19, 082317, https://doi.org/10.1063/1.4748296, 2012. a
Nariyuki, Y., Hada, T., and Tsubouchi, K.: Collisionless Damping of Circularly Polarized Nonlinear Alfvén Waves in Solar Wind Plasmas with and without Beam Protons, Astrophys. J., 793, 138, https://doi.org/10.1088/0004637X/793/2/138, 2014. a
Ofman, L.: Wave modeling of solar wind, Living Rev. Sol. Phys., 7, 846–855, https://doi.org/10.12942/lrsp20104, 2010. a
Sagdeev, R. Z. and Galeev, A. A.: Nonlinear plasma theory, edited by: O'Neil T. and Book, D., W. A. Benjamin, New York, p. 7, 1969. a
Spangler, S. R., Leckband, J. A., and Cairns, I. H.: Observations of the parametric decay instability of nonlinear magnetohydrodynamic waves, Phys. Plasmas, 4, 846–855, https://doi.org/10.1029/2001JA000150, 1997. a
Terasawa, T., Hoshino, M., Sakai, J.I., and Hada, T.: Decay instability of finiteamplitude circularly polarized Alfven Waves: A numerical simulation of stimulated Brillouin scattering, J. Geophys. Res., 91, 4171–4187, https://doi.org/10.1029/JA091iA04p04171, 1986. a, b, c
Tu, C.Y. and Marsch, E.: Magnetohydrodynamic structures waves and turbulence in the solar wind – observations and theories, Space S. Rev., 73, 1–210, https://doi.org/10.1007/BF00748891, 1995. a
Tu, C.Y. and Marsch, E.: Anisotropy regulation and plateau formation through pitch angle diffusion of solar winf protons in resonance with cyclotron waves, J. Geophys. Res., 107, 1249, https://doi.org/10.1029/2001JA0000150, 2002. a
Verscharen, D. and Marsch, E.: Apparent temperature anisotropies due to wave activity in the solar wind, Ann. Geophys., 29, 909–917, https://doi.org/10.5194/angeo299092011, 2011. a, b
Verscharen, D., Marsch, E., Motschmann, U., and Müller, J.: Parametric decay of oblique Alfvén waves in twodimensional hybrid simulations, Phys. Rev. E, 86, 027401, https://doi.org/10.1103/PhysRevE.86.027401, 2012. a
Viñas, A. F. and Goldstein, M. L.: Parametric instabilities of circularly polarized largeamplitude dispersive Alfvén waves: excitation of parallelpropagating electromagnetic daughter waves, J. Plasma Phys., 46, 107–127, https://doi.org/10.1017/S0022377800015981, 1991a. a, b, c, d, e
Viñas, A. F. and Goldstein, M. L.: Parametric instabilities of circularly polarized largeamplitude dispersive Alfvén waves: excitation of obliquelypropagating daughter and sideband waves, J. Plasma Phys., 46, 129–152, https://doi.org/10.1017/S0022377800015993, 1991b. a, b, c, d, e, f, g