Ion acceleration by parallel propagating nonlinear Alfvén wave packets in a radially expanding plasma

The numerical simulation of the nonlinear evolution of the parallel propagating Alfvén waves in a radially expanding plasma is performed by using a kinetic-fluid model (the Vlasov–MHD model). In our study, both the nonlinear evolution of the Alfvén waves and the radial evolution of the velocity distribution function (VDF) are treated simultaneously. On the other hand, important ion kinetic effects such as ion cyclotron damping and instabilities driven by the nonequilibrium ion velocity distributions are not included in the present model. The results indicate that the steepened Alfvén wave packets outwardly accelerate ions, which can be observed as the beam components in the interplanetary space. The energy of imposed Alfvén waves is converted into the longitudinal fluctuations by the nonlinear steepening and the nonlinear Landau damping. The wave shoaling due to the inhomogeneity of the phase velocity is also observed.


Introduction
It is well known that nonlinear Alfvén waves are ubiquitously observed in the solar wind (Bruno and Carbone, 2013).These Alfvén waves are probably generated in the vicinity of the solar photosphere, and may play an important role in heating and acceleration of the coronal and solar wind plasmas.A lot of authors have already discussed the damping of the Alfvén waves in coronal plasmas (Ofman, 2010a;Galinsky and Shevchenko, 2013).In particular, the nonlinear damping of low-frequency, quasi-parallel propagating Alfvén waves is of great interest.This is because they can carry the momentum and energy over a long distance, since the linear col-lisionless damping of these waves is too small to dissipate them (Barnes, 1966).
On the other hand, it is also known that the proton beams are frequently observed in solar wind plasmas (Marsch, 2006).Recently, it was shown by Araneda et al. (2008) and Matteini et al. (2010) that the protons trapped by the envelope modulation (or corresponding longitudinal fluctuations) due to the parametric instabilities of Alfvén waves can be observed as the beam protons in the velocity space (Araneda et al., 2008).Past studies (Araneda et al., 2008;Matteini et al., 2010) assumed the uniform background plasmas and the periodic system.However, the radial inhomogeneity of the plasma density and the magnetic field may not be negligible to discuss the nonlinear evolution of low-frequency Alfvén waves in coronal plasmas, since the radial change of the nonlinearity and the non-WKB effects such as the reflection may occur (e.g., Chandran et al., 2011).
In the present paper, we firstly perform the numerical simulation of the nonlinear evolution of the parallel propagating Alfvén waves in a radially expanding plasma by using a kinetic-fluid model (the Vlasov-MHD model) (Nariyuki and Hada, 2007;Nariyuki et al., 2011Nariyuki et al., , 2013)).Actually, this has already been discussed by a lot of authors by using the magnetohydrodynamic (MHD) models (Suzuki and Inutsuka, 2006;Tanaka et al., 2007).In our study, the radial evolution of the Alfvén waves and the velocity distribution function (VDF) are treated simultaneously.Even if the isotropic VDFs are assumed in the Vlasov-MHD model, as far as the nonlinearity is relatively weak (δb/b <∼ 0.3) (Nariyuki and Hada, 2007), the nonlinear evolution of the Alfvén waves and the trapped ions observed in the Vlasov-MHD model

Y. Nariyuki et al.: Ion acceleration by nonlinear Alfvén wave packets
agree with those in the ion hybrid simulation (Nariyuki et al., 2011).In Sect.2, the one-dimensional (1-D) Vlasov-MHD model in the spherical coordinate and the simulation setup are introduced briefly.Section 3 contains the numerical results.We discuss and summarize the results in Sects.4 and 5, respectively.

Simulation setup
We consider a radially expanding flux tube, which is measured by the heliocentric distance r.To discuss both the nonlinear evolution of the transverse fluid motion due to the Alfvén waves and the longitudinal (radial) ion kinetics, we introduce the 1-D Vlasov-MHD model in the spherical coordinate below.For simplicity, we here neglect the gravitational force term, since it is relatively small and may not play important roles in the super Alfvénic solar wind at r > 10R s , where R s is the solar radius.
The velocity distribution functions obey the Vlasov equa- where f (x, v, t) is the single particle distribution function, and In the spherical coordinate (x s = (r, θ, φ), where x = r sin θ cos φ, y = r sin θ sin φ, z = r cos θ), dv dt can be written as where i j is the unit vector of the j = r, θ, φ component, v = v s = (v r , v θ , v φ ) = ( dr dt , r dθ dt , r sin θ dφ dt ), and a = a r i r +a θ i θ + a φ i φ .Notice that v s is the function of x.By transforming (1) into the six-dimensional spherical coordinate ((x s , v s )), we find (Larson, 1969) where The relation ( 5)-( 7) can be derived from Eq. ( 2).We here define the longitudinal distribution function We assume the spherical symmetry ( ∂f ∂θ = ∂f ∂φ = 0) and isotropic plasmas.By integrating Eq. (1) on v θ and v φ , we obtain where u j is the j component of the bulk velocity, and n = gdv r , the number density.By taking the boundary condition f → 0 when v θ,φ → ±∞, where Equation ( 11) is valid in isotropic plasmas (Nariyuki et al., 2013).We here take the ion mass m i = 1, q = 1, and where p i = nT i = (v r − u r ) 2 f dv = (v r − u r ) 2 gdv r .When electrons can be treated as the massless fluid, we obtain where and p e is the electron pressure.Then, we find that ( The first term on the right-hand side of Eq. ( 9) comes from the sixth and seventh terms on the left-hand side of Eq. ( 4).Notice that Eq. ( 9) is obtained for arbitrary f .While the second term on the right-hand side of Eq. ( 9) comes from the sixth term, we here neglect it by taking θ = π/2.
where b = b θ + ib φ and u = u θ + iu φ are the complex transverse magnetic field and velocity, By integrating Eq. ( 15) on v r , we recover the continuity equations and the momentum equation of the radial component in the ideal MHD system (Suzuki and Inutsuka, 2006;Tanaka et al., 2007).The ion pressure p i is also evaluated by the integration of g.Isotropic ions are assumed in Eqs. ( 15)-( 17) to close the system.Note that the assumption of isotropic ions is reasonable as far as we concentrate on the low-frequency wave dynamics in low beta plasmas, although the VDFs of ions in the coronal and solar wind plasmas may be anisotropic.This is because the pressure anisotropy effects in the dispersion relation of Alfvén waves are negligibly small in low beta plasmas (see Sect. 4.2).
The simulation region of the present study is from 10R s to 110R s .All the physical quantities are normalized to the parameters at r = 10R s as follows: the number density n 0 = ρ 0 /m p = 2.9 × 10 9 [m −3 ] (m p is the mass of protons), the longitudinal magnetic field b r0 = 1.3 × 10 −6 [T], the Alfvén velocity V A0 = b r0 / √ (µ 0 ρ 0 ) = 5.0 × 10 5 [m s −1 ] (µ 0 is the vacuum permeability), respectively.The space and time are normalized to the solar radius R s = 7 × 10 8 [m] and T 0 = R s /V A0 = 1.4 × 10 3 [s], respectively.We use the uniform mesh with 4000 grid points in the r direction and 1000 grid points in the v r direction.The grid spacing is r = 0.025, v r = 0.01 and the time step is t = 2.5 × 10 −4 .The transverse momentum Eq. ( 16) and the induction Eq. ( 17) are solved by the second-order rational Runge-Kutta scheme (Wambecq, 1978) for time integration and the second-order central difference method for evaluating spatial derivatives.The splitting scheme (Cheng, 1976) and the PIC (Positive Interpolation for hyperbolic Conservation laws) scheme are applied for time advancement of the Vlasov Eq. ( 15) (Umeda, 2009;Umeda et al., 2012).Note that no artificial dissipation term is included in the Vlasov-MHD system.The detailed numerical scheme to solve Eqs. ( 15)-( 17) has also been described in a past study (Nariyuki et al., 2011).
Since the solar wind may become super Alfvénic flow within r < 10 (r = 10 is the inner boundary at 10R s ), we here assume the constant radial velocity u r0 = 1 ρ v r g 0 dv r = V A0 as the initial condition, where g 0 is the initial VDF.Then, the initial density ρ = 100/r 2 due to the conservation of the mass flux.Thus, the critical point at which the solar wind becomes super Alfvénic flow is r c = 10 in the present study.The radial magnetic field is b r = 100/r 2 due to the conser- vation of magnetic flux.Note that u r > V A = b r / √ (µ 0 ρ) at r > 10.Thus, the initial ion pressure p 0 is assumed to be uniform.The initial ion beta is The nonmonochromatic right-handed circularly polarized Alfveń waves and v(r = 10) = −b(r = 10) are imposed at the lower boundary, where The imposed Alfvén waves are given as the white noise (b p is independent of ω n ).We here set ω 0 = 12.566 (ω 0 /T 0 = 0.009 [Hz]), and N = 60 (0.0076 [Hz] < ω/T 0 < 0.0104 [Hz]) (Run A), N = 100 (0.0067 [Hz] < ω/T 0 < 0.0113 [Hz]) (Run B), N = 140 (0.0058 [Hz] < ω/T 0 < 0.0122 [Hz]) (Run C).At the early stage (t < 10π/ω 0 ), the linearly increasing wave amplitude a = t (ω 0 /10π ) is set to avoid the discontinuity of the wave front.During 10π/ω 0 < t < 11, we set a = 1 and b p as < |b(r = 10)| >= 0.2 (where <> indicates time average), which is consistent with previous studies (e.g., Maneva et al., 2013).After t > 11, the imposed wave amplitude is diminished.

Results
In this section, we mainly show the results of Run B, since those of runs A and C are essentially the same as those of Run B. The difference between the runs is discussed in the last part of this section.
Figure 1   the density fluctuation by the ponderomotive force (the second term of the right-hand side of Eq. 14) (Mjølhus, 1976;Mio et al., 1976).Such longitudinal fluctuations can dissipate through the nonlinear Landau damping (Nariyuki et al., 2008(Nariyuki et al., , 2010;;Markovskii et al., 2009).Figures 2 and 3 show the snapshots of the transverse magnetic field (b ) and g at (a) t = 16.5,(b) t = 19.25, and (c) t = 22.0, respectively.As time elapses, nonlinear steepening of the envelopemodulated Alfvén waves (Cohen and Kulsrud, 1975) is observed (Figs. 1 and 2).In Fig. 3, the ions accelerated by the nonlinearly steepened wave front of the Alfvén waves are observed.Some ions with the higher velocity than the local phase velocity of the Alfvén waves move away from the steepened wave front and form the stretched structure in the phase space (Fig. 3b, c).It is important that the accelerated ions can be observed as the ion beams by the spacecraft observation (Fig. 4).For instance, since V A 10V A0 /r, the relative beam speed in Fig. 4a is v r = 0.5 1.8VA (r = 36), which is consistent with the one observed in solar wind plasmas (Marsch, 2006).In the one-fluid system, the bulk kinetic energy corresponding to the relative speed between core and beam protons is evaluated as the pressure and temperature parallel to the ambient magnetic field.For instance, the temperature evaluated by using the whole proton VDFs in Fig. 4a is about 1.5 times larger than the one of the core proton.On the other hand, the core proton temperature at t = 16.5 and r = 36 (Fig. 3a) is 0.0525, while the one at t = 22.0 and r = 36 is 0.0450.This suggests that the heating of ions by the nonlinear Alfvén waves in the present study corresponds to the generation of the relative bulk kinetic energy, which can be free energies to cause the beam instabilities.Actually, the present model also does not include the other important kinetic effects.The quantitative results of heating in the present study may be modified by the kinetic effects.To incorporate the instabilities into the present model through the anisotropic pressure is one of the important future issues.As shown in Fig. 7, the intensity of the sunward-propagating waves (z m ) is much less than the anti sunward-propagating waves (z p ). Namely, the non-WKB reflection by the background inhomogeneity (Suzuki and Inutsuka, 2006;Chandran et al., 2011) is weak, and the parametric decay instability of Alfvén waves (Tanaka et al., 2007) does not occur.The absence of the non-WKB reflection indicates that the WKB (Wentzel-Kramers-Brillouin) approximation is well satisfied in the present study.The absence of the parametric decay instability of Alfvén waves is discussed in the next section.

Wave shoaling
In this subsection, we discuss the increase in nonlinearity in the WKB Alfvén wave packets in the solar wind.It is well known that the envelope-modulated Alfvén waves are steepened by the nonlinearity (Cohen and Kulsrud, 1975).The steepening converts the energy of transverse waves into the one of longitudinal waves (parallel ion kinetic energy), as shown in the previous section.This result suggests that the simple steepening (parallel cascading) and nonlinear ion Landau damping can play an important role in damping of Alfvén waves.On the other hand, the modulational instability does not occur due to the absence of the wave dispersion in the present system (Sakai and Sonnerup, 1983).Notice that the absence of the modulational instability does not affect the decrease in energy of envelope-modulated Alfvén waves (Nariyuki et al., 2008).In contrast to the previous studies (Suzuki and Inutsuka, 2006;Tanaka et al., 2007) in which the MHD system was used, the decay instability is not clearly observed in the present study.This is consistent with the previous studies (Nariyuki et al., 2008;Markovskii et al., 2009) in which the ion hybrid simulations in uniform plasmas were carried out.Actually, the growth rate of the decay instability is small due to the small amplitude of each wave number mode in the present study.Moreover, the growth rate of the decay instability is weakened by the ion Landau damping of the longitudinal waves (Nariyuki and Hada, 2007).By comparing both Figs. 1 and 2c, the contraction of the wave length (2π/k) is also observed.This is due to the conservation of the wave frequency where the phase velocity u r + V A decreases with increasing heliocentric distance.When the wave packets satisfy the WKB approximation, the nonlinearity of the waves (Barnes, 1992) where the constants C 0 and C 1 are determined by using the boundary condition.Equation ( 21) indicates that the nonlinearity of the wave packets increases with increasing heliocentric distance.The increase in the nonlinearity with the contraction of the wave length is similar to the wave shoaling of the surface gravity waves.When u r >> V A (r >> C 1 ), ω u r k and |b| b r 2 ∝ r.In this limit, the wave number k is the quasi-invariant and the waves can be described by using the expanding box (Grappin et al., 1993;Grappin and Velli, 1996;Liewer et al., 2001;Hellinger and Travnicek, 2005;Ofman et al., 2011).On the other hand, if u r << V A (r << C 1 ), k ∝ r and |b| b r 2 ∝ r 3 .Thus, the nonlinearity |b| b r at a certain r becomes larger with increasing heliocentric distance from the critical point r c .This effect of inhomogeneity promotes the nonlinear steepening and the ion acceleration by ponderomotive force.
The polarization of the imposed waves does not affect the time evolution of the system, since the wave dispersion due to the Hall effect and the finite Larmor radius effect is neglected in the present model.Notice that the nonlinear evolution of low-frequency Alfvén waves generates the higher frequency waves.The linear damping (ion cyclotron damping) and nonlinear damping (parametric instabilities) of the dispersive Alfvén waves are affected by the polarization (e.g., Liewer et al., 2001;Hellinger et al., 2005;Matteini et al., 2010).Actually, while the polarization of Alfvén waves does not affect the time evolution of the present dispersionless system, dispersive left-handed porarized waves are unstable to both decay and modulational instabilities in low beta plasmas (Longtin and Sonnerup, 1986).To incorporate the non-MHD effects (Chandran et al., 2011) such as the kinetic effects discussed in the next subsection into the present one-fluid model is an important future issue.

Ion kinetics
In the uniform and periodic system, the phase space vortex corresponding to the density fluctuations has also been observed in the nonlinear stage (Araneda et al., 2008;Matteini et al., 2010).In the present system, ions accelerated by the ponderomotive force are observed as the rolled-up VDFs in the phase space, while the vortex structures are not observed.This may be because of the asymmetrical wave form observed in Fig 2 .Note that the trapped ions in the uniform system are the pseudo beam in the phase space, while the detrapped ions observed in the present study can survive as the real beam.However, the linear stability of the beams could not be discussed in the present system, since inhomogeneity of the velocity space is decoupled from the transverse wave dynamics due to the assumption of the isotropic velocity distribution function.Notice that the assumption of the isotropic VDFs in the present model is only valid for low beta plasmas such as those in the vicinity of the sun.Although the protons observed in the solar wind have the temperature anisotropy (Marsch, 2006;Bale et al., 2009;Maruca et al., 2011), the phase velocity of Alfvén waves evaluated by the Chew-Goldberger-Low (CGL) system is close to the Alfvén velocity in low beta plasmas (Nariyuki et al., 2013).The phase velocity of anti-sunward propagating Alfvén waves in CGL system is where p and p ⊥ are the pressure parallel and perpendicular to the total magnetic field b.Equation ( 22) indicates that even if temperature (and pressure) anisotropy is large, v φ u r + V A in low beta plasmas.For instance, by using an empirical law of fast solar wind plasmas in Hellinger et al. (2011), v φ − u r 1.019VA (r = 20), 1.020VA (r = 40), and 1.018VA (r = 60), respectively.On the other hand, as discussed below, some important kinetic effects are not included in the present model.Of course, even if the effect of temperature anisotropy on the phase velocity of Alfvén waves is negligible, the various instabilities can be caused by temperature anisotropy (Matteini et al., 2013).
In the present study, Alfvén wave packets excite the longitudinal fluctuations (deformation of the VDFs).Maneva et al. (2013) have recently confirmed the similar phase space vortex in the expanding box model.Moreover, nonresonant scattering of ions neglected in the present model can also affect the local VDFs (Araneda et al., 2008;Dong and Singh, 2013).In addition, the wave-particle resonant interactions between the Alfvén wave packets and ions become important with increasing heliocentric distance, since the ion cyclotron frequency decreases (Liewer et al., 2001).The Alfvén wave steepening observed in the present study causes the parallel cascading in the frequency space and may promote the ion cyclotron damping.It is also noteworthy that the kinetic instabilities driven by the beam components could not be described in the present system.The beam instabilities are affected by the inhomogeneity (e.g., Ofman, 2010b;Hellinger and Travnicek, 2011).To fill the gap between the past studies with the uniform system and the present study, further implementation of kinetic simulations with mesoscale models such as the expanding box models (Liewer et al., 2001;Hellinger and Travnicek, 2005;Ofman et al., 2011;Maneva et al., 2013) is required.

Conclusions
In the present paper, we discuss numerically the ion acceleration by nonlinear Alfvén wave packets in the radially expanding plasma.The results indicate that the steepened Alfvén wave packets outwardly accelerate ions, which can be observed as the beam components in the interplanetary space.The wave shoaling of the WKB Alfvén waves is also observed.On the other hand, as discussed in the previous section, important ion kinetic effects such as ion cyclotron damping and instabilities driven by the non-equilibrium VDFs are not included in the present model and the quantitative results of the present study may be modified by the kinetic effects.Actually, the effective heating and acceleration of ions by the Alfvén waves are very weak in the present study, while the ion beams are formed corresponding to the envelope modulation.The non-equilibrium VDFs observed in the runs still preserve the free energy to cause the instabilities.In order to address the effective ion heating of solar wind plasmas by low-frequency Alfvén waves, the kinetic effects and the other non-MHD effects such as turbulent cascading and collisions should be incorporated into the present model.
is the particle velocity, e = (e x , e y , e z ) the electric field, b = (b x , b y , b z ) the magnetic field, m the rest mass, and q the charge.
) ρ = r 2 ρ, u r the longitudinal velocity, p e the electron pressure, p b = |b| 2 /2, and p e(b) = r 2 p e(b) .C = b r r 2= 100, since the physical quantities are normalized to the parameters at the inner boundary as shown below.We here assume the isothermal electrons (p e /ρ = const.).
50, n is an integer, and N is the number of wave modes, respectively.The phase α n is given by a uniform random number generator with a range [0, 2π].The inner boundary condition for g is fixed at the initial value.The open boundary condition is imposed at the upper boundary (r = 110).Initially, β e = 0.02.
shows the transverse magnetic field (b θ , b φ , ±|b |) at t = 11.The imposed Alfvén wave packets have the envelope modulation corresponding to the broadband frequency spectra.Notice that the envelope modulation excites

Fig. 5 .
Fig. 5.The snapshot of (a) the transverse magnetic field and (b) g at t = 22.0 in Run A.

Figures 5
Figures5 and 6show the snapshot of the transverse magnetic field and g at t = 22.0 in runs A and C, respectively.In the resultant VDFs of all the runs (Figs.4c, 5b, and 6b), the accelerated ions can be observed.The number of generated beams increases with increasing frequency range of the imposed Alfvén waves, since the wave length of the envelope modulation becomes smaller.

Fig. 6 .
Fig. 6.The snapshot of (a) the transverse magnetic field and (b) g at t = 22.0 in Run C.

Figure 7
Figure 7 shows the Elsasser variables z p = u y − b y / √ ρ and z m = u y + b y / √ ρ at t = 22.0.As shown in Fig.7, the intensity of the sunward-propagating waves (z m ) is much less than the anti sunward-propagating waves (z p ). Namely, the non-WKB reflection by the background inhomogeneity(Suzuki and Inutsuka, 2006;Chandran et al., 2011) is weak, and the parametric decay instability of Alfvén waves(Tanaka et al., 2007) does not occur.The absence of the non-WKB reflection indicates that the WKB (Wentzel-Kramers-Brillouin) approximation is well satisfied in the present study.The absence of the parametric decay instability of Alfvén waves is discussed in the next section.

Fig. 7 .
Fig. 7.The Elsasser variables z p = u y − b y / √ ρ (gray solid line) and z m = u y +b y / √ ρ (black dashed line) at t = 22.0 in (a) Run A, (b) Run B, and (c) Run C.