On heating of solar wind protons by the breaking of large amplitude Alfvén waves

By means of hybrid simulations, we present a study on plasma heating by the field-aligned parametric decay of a monochromatic left-handed polarized Alfvén wave. Simultaneous multidimensional comparisons of the wave modes and proton kinetics suggest that parametric decay of Alfvén waves and pitch angle scattering of solar wind protons are interrelated. Parametric decay mechanism yields counter-propagating Alfvén waves that can shape and broaden via pitch angle scattering mechanism 5 both the sunward and antisunward sides of the proton velocity distribution functions in agreement with in situ measurements of fast stream solar wind plasmas.


Introduction
Early in situ measurements at 1AU from VELA satellite (Bame et al., 1975) reveal that the velocity distribution function of solar wind protons is broader in direction perpendicular to the mean magnetic field compared to the parallel distribution, thus indicating that the perpendicular temperature is higher than the parallel temperature.Marsch et al. (1982) have found by using Helios 1 and Helios 2 data such plasma heating occuring in high speed solar wind streams from 0.3 to 1.0 AU.The problem of acceleration and heating of ions in solar wind and solar corona is vast and has been discussed for a long time in space plasma physics, see e.g., the reviews of Marsch (2006); Ofman (2010).Theoretical models based on cyclotron resonant-or non-resonant-processes have been proposed to explain the heating of solar corona and solar wind.All of these models are related to the pitch angle scattering of protons. of solar wind protons, driven by resonance with ion cyclotron waves propagating outwards the Sun.The ions in resonance with transversal cyclotron waves conserve their total kinetic energy in a frame moving with the phase velocity of the wave while their velocity distribution functions show a plateau observed through circular contours of constant density in phase space.Thus, the perpendicular broadening of the left hand part of the measured distributions (the sunward part) has been explained through the pitch angle scattering of solar wind protons resonantly interacting with the outward propagating waves.However, this process cannot explain the antisunward part of the core distribution.Araneda et al. (2007) suggest, on basis of the conclusions discussed by Marsch and Tu (2001), that a candidate which is able to elucidate this mystery could be the parametric decay instability, a process so rarely detected by in situ measurements in space plasmas (see e.g., Spangler et al., 1997).
Parametric instabilities play an important role in the dissipation of the large amplitude 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 (see e.g., Viñas and Goldstein , 1991) or numerical magnetohydrodynamics MHD (e.g., Ghosh et al., 1993;Ghosh and Goldstein, 1994;Del Zanna et al., 2001) and particle-in-cell or hybrid (e.g., Terasawa et al., 1986;Matteini et al., 2010a, b;Verscharen et al., 2012;Gao et al., 2013;Nariyuki et al., 2014) simulations.In the MHD picture, the plasma heating by the Alfvén wave can occur through generation and steepening of magnetosonic 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 non-linear trapping of protons (see e.g., Araneda et al., 2008;Matteini et al., 2010a, b).By using hybrid simulations, Nariyuki et al. (2012) show that monochromatic Alfvén waves with circularly left-handed polarization are quickly dissipated by a beam driven by kinetic Alfvén waves.Nariyuki et al. (2014) also discuss on the collisionless damping of the field aligned Alfvén waves at low frequency in solar wind plasmas with or without the presence of a proton beam.The velocity beam formation is however restricted by the conditions of beta plasmas.On the assumption of cold electrons, the electric field is decoupled from the density fluctuations and the formation of the beam is inhibited (see e.g., Matteini et al., 2010a).By means of 2D hybrid simulations, Gao et al. (2013) recently suggest that under such conditions relevant for the solar corona, the obliquely-propagating Alfvén daughter waves driven by the parametric decay could stochastically heat the plasma.
Here we present a numerical study on plasma heating by the decay of a large amplitude monochromatic Alfvén wave propagating along the mean magnetic field in a multidimensional simulation box immersed in a uniform low beta electron proton collisionless plasma.By analyzing the velocity distribution functions and the time history of the wave modes obtained from one-, two-, and three-dimensional hybrid simulations, we report that parametric decay and pitch angle scattering mechanisms are both involved in broadening the entire proton velocity distribution in direction perpendicular to the mean magnetic field.
The left-handed circularly polarized Alfvén pump wave with forward propagation does perpendicularly broaden one side of the particle velocity distributions while the backward propagating Alfvén daughter wave enlarges the other side, respectively.

Hybrid simulation and discussion on results
The hybrid simulations are conducted in 1D, 2D and 3D setups by using the hybrid code AIKEF (Mueller et al., 2011).The one-dimensional box has a length of L 1 =720 d i while in the two-dimensional setup, L 1 =L 2 =288 d i on each direction.Here, is the ion inertial length while V A and Ω p are the Alfvén velocity and the ion gyrofrequency, respectively.A number of N c =10000 super-particles is introduced in each computational cell.By adding the third dimension (L 3 =288 d i ), the number of super-particles is diminished to N c =1000 particles in a cell.Even so, the total number of super-particles used in the 3D system (N t ≈ 2.0 • 10 10 ) is at the limit of our computational resources but it is sufficient to reduce the influence of numerical noise in agreement with additional 3D tests carried out by using smaller sizes of the box (L 1 =L 2 =L 3 =144 d i ), and a larger number of super-particles in the simulation cell (N c =2500).The grid size is the same for all the cases and equals the ion inertial length (d i ).A low value of β=0.01 is used for the plasma beta parameter for each species of particles, namely the protons treated in the hybrid scheme as particles and the electrons considered as a massless fluid.At the initial time, a monochromatic Alfvén pump wave with left-handed circularly polarization is launched parallel to the mean magnetic field along the Oz axis of the simulation box.The amplitude of the pump wave is set to 20% of the mean magnetic field magnitude, δB/B 0 =0.2.The wavenumber of the pump wave is k 0 V A /Ω p = 0.218, corresponding to a Fourier mode of m=10 (m=25 for the 1D system).
The parametric decay of the Alfvén wave is a 3-wave process involving a large amplitude pump wave, a longitudinal daughter wave with the same direction of propagation as the pump wave and a counter-propagating transverse daughter wave.Figure 1 shows the time profiles for the magnetic field energy of the pump wave (labeled by E P ) and the backward propagating Alfvén daughter waves (labeled by E D ) for the 2D setup (left panel) and the 3D setup (right panel).Overplotted by gray solid lines are the root mean square (rms) density fluctuations (ρ) with a maximum at about 300 ion gyroperiods when the linear growing of E D terminates, i.e., the linear saturation of the instability occurs.At the nonlinear saturation time when the energy of the daughter wave overtakes the energy of the mother wave, the frequency-wavenumber power spectrum of the magnetic field is shown in the lower panels of Fig. 1.The left-handed pump wave (here we adopted the convention of negative frequency for the left-handed mode) marked with E P can be seen in the lower-right quadrant (k > 0 and ω < 0) while the counter-propagating Alfvén daughter wave (E D ) is developed in the lower-left quadrant (k < 0 and ω < 0) of the frequency-wavenumber spectrum and has the Fourier mode m=-8 (m=-23 for the 1D run, not shown here).The waves driven by the decay instability dominate the entire spectrum.Qualitatively speaking, the results presented in Fig. 1 are consistent with the earlier MHD studies, proving that the parametric decay responsible for the breaking of the Alfvén pump wave occurs irrespective of the spatial dimensions.
Figure 2 shows the particle distribution functions in the phase space z − v at two different stages for the evolution of decay instability in the 2-D and 3-D setup.The upper panels of Fig. 2 refer to the linear stage of the instability at a time of tΩ p ∼ 300 and the lower panels are obtained for the latest time of the simulation run (tΩ p = 600) at the nonlinear phase of the instability.
Here, the number of particles is represented in a color coded scale with a minimum defined by the lightest nuance of gray.At the linear stage, the protons are accelerated by the parallel electric field produced by the density fluctuations.The proton phase space z−v is close to the result from an earlier study.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 large amplitude Alfvén mother wave.The spatial modulation in Fig. 2 is weaker according to the smaller amplitude pump wave used in our simulation.During the saturation of the instability we do not observe the presence of phase-space vortices leading to the formation of the proton beam due to most probably the small value of electron beta used in simulation.The proton phase space at the linear stage of the instability does not differ too much between the 2D and 3D runs.The velocity distribution, in contrast, is different for the 3D system at the nonlinear stage of the instability as one can see in the lower panels of Fig. 2. Particles are smoothly heated and the resonant protons accelerated by the ion acoustic waves are mixed with the thermal core of the distribution.
Figure 3 reports a comparison of the proton velocity distribution functions constructed in the (v ⊥ ,v ) plane obtained from the one-, two-, and three-dimensional systems at the same simulation epochs as in Figure 2. Here, the perpendicular velocity of the particle is taken as v ⊥ = v 2 x + v 2 y .The contour lines are shown for fractions of 0.6, 0.3, 0.2, and 0.1 of the maximum phase space density from the inner to the outer part of the distribution functions.One firstly should note that the shape of the perpendicular velocity distribution (as a parallel cut to the v ⊥ axis) still preserve the remnant of the rigid displacement of the whole velocity distribution function by the transversal wave field imposed as an initial condition, see e.g., Verscharen and Marsch (2011);Nariyuki (2011).However, the time evolution of velocity distribution function is different for the three analyzed systems suggesting that the initial shifted Maxwellian form would not significantly disturb the processes triggered by the decay instability at later times.Similar with Fig. 2, the velocity distribution functions obtained from the 2D and 3D setups do not clearly differ close to the linear stage of the decay instability at time tΩ p =300 (upper panels).The particles in the uni-dimensional setup experience the strongest parallel acceleration.The final distribution functions for the 3D system (panel 3f) report a larger perpendicular acceleration suggesting a stronger perpendicular heating of the protons in comparison with the 1D and 2D systems.The time history of the ion temperature (not presented here) is consistent with this result showing a significant increasing of the perpendicular temperature with respect to that obtained from the 1D and 2D runs at the latest times of our simulations.
The velocity distribution functions (VDF) in Fig. 3 (at tΩ p =600) indicate that the parallel wave scattering is more efficient in broadening the phase space density for the 1D system.The perpendicular wave scattering, in contrast, is effective in the threedimensional run.The dashed lines drawn at time tΩ p =600 in panel 3f represent arcs of circles with centers localized along the v axis at the position corresponding to the phase velocities of Alfvén pump wave (V ph = ω 0 /k 0 V A ≈-0.9) and Alfvén daughter wave (V ph = ω/kV A ≈+0.9), respectively.They describe the particle velocity in the wave frame, v = (v − V ph ) 2 + v 2 ⊥ .According to pitch angle scattering model, the particle energy, E=mv 2 /2, is conserved in this reference system.Identical dashed lines are plotted for comparison in the other five panels of Fig. 3.We see that these arcs coincide with the most obliquely parts of contour lines while the outer contours are better overlapped than the inner ones.The overlapping is larger in comparison with that observed at earlier time (tΩ p =300 in panel 3c).We expect that at later times (tΩ p >600) when the effects of the initial wave field should disappear, the matching will be improved.By comparing panel 3c and panel 3f, we note that the time evolution of the VDF is primarily controlled by the pitch angle scattering mechanism.The inner contour line corresponding to 0.6 of the maximum phase space density (panel 3c) evolves closely bounded by the dashed line satisfying the particle energy conservation in the wave field referential (panel 3f).In contrast, panel 3d corresponding to the 1D system does not have this feature while panel 3e related to the 2D setup shows a better agreement for the left part (negative parallel velocities) of VDF.
At time tΩ p =600, the contour line corresponding to fraction of 0.1 of maximum phase space density is best described by the pitch angle scattering model in panel 3f and partially in panel 3e.Again, panel 3d shows disagreement due to the larger parallel scattering developed in the 1D system.
Another interesting result from the upper panels in Fig. 1 is that the E D fluctuation energy of daughter mode decreases We studied the proton acceleration and heating driven by the parametric decay of a large amplitude Alfvén wave in the linear and nonlinear stage of the instability in a multidimensional system.By comparing the wave modes and proton velocity distribution functions, we conclude that the pitch angle scattering on the parallel propagating daughter-and mother-Alfvén waves is the key mechanism which describes the perpendicular broadening of the particle velocity distribution functions.In the 3D system, the conditions for resonant cyclotron interactions of ions with the left-handed Alfvén waves are apparently fulfilled especially at the non-linear stage evolution of the decay process.Further studies are needed for a better understanding of this circumstance by using a larger panel of initial resonant conditions for the particles and wave fields, or different values of the ion and electron temperatures.Tsurutani et al. (2005) discuss on a close relation between the perpendicular acceleration of protons, the phase-steepened Alfvén waves, and the magnetic field decreases observed preponderantly in the fast streams of interplanetary solar wind.
Parametric decay and pitch angle scattering could play also a role for a better understanding of the connection between the mentioned observations.Nevertheless, the result discussed in this paper can invoke that that parametric decay, besides its implication in the reduction of the cross-helicity in the solar wind and in the local production of turbulence, would also play a key role in pitch angle scattering of the solar wind protons and finally in the perpendicular temperature anisotropy observed in the inner heliosphere.
Competing interests.The authors declare that they have no conflict of interest.
Acknowledgements.This work is financially supported by a grant of the Deutsche Forschungsgemainschaft (DFG grant MO539/20-1).
HC acknowledges the hospitality at University of Toyama for hosting the research visit.We acknowledge John von Neumann Institute for

Figure 1 .
Figure 1.Top: Time evolution of the normalized magnetic field energy of the Alfvén pump wave (E P ) and the counter propagating Alfvén daughter wave (E D ).Overplotted by gray are given the rms density fluctuations ρ.Bottom: Power spectrum in the wavenumber frequency domain of the magnetic field determined around time tΩp ∼ 500.

Figure 2 .
Figure 2. Top: Proton phase space z − v at time tΩp = 300 close to the linear saturation of the decay instability.Bottom: Proton phase space z − v at the final time tΩp = 600 of the simulation.

Figure 3 .
Figure 3. Proton distribution functions determined in the plane (v ⊥ ,v ) at the same times as in Fig. 2.