Computational and theoretical study of the wave-particle interaction of protons and waves

We study the wave-particle interaction and the evolution of electromagnetic waves propagating through a plasma composed of electrons and protons, using two approaches. First, a quasilinear kinetic theory has been developed to study the energy transfer between waves and particles, with the subsequent acceleration and heating of protons. Second, a one-dimensional hybrid numerical simulation has been performed, with and without including an expanding-box model that emulates the spherical expansion of the solar wind, to investigate the fully nonlinear evolution of this wave-particle interaction. Numerical results of both approaches show that there is an anisotropic evolution of proton temperature.


Introduction
The problem of acceleration and heating of plasmas due to the interaction of particles with propagating waves has received special attention during the last decades, especially in the field of laboratory and space plasma physics.In the case of the solar wind (Axford andMcKenzie, 1992, 1996;Kohl et al., 1998;Marsch, 1998;Cranmer et al., 1999a,b;Esser et al., 1999;Hu and Habbal, 1999;Tu and Marsch, 1999;Cranmer, 2002), recent observations and theoretical results seem to indicate that most of the acceleration process occurs within a few solar radii from the Sun and the main mechanism is due to resonant absorption of ion-cyclotron waves (Isenberg, 2001;Cranmer, 2002;Hollweg and Isenberg, 2002).However, the detailed processes for the energy transfer between waves and different particle species are still an open question.To address these issues, we investigate the wave-particle interaction and evolution of circularly polarized electromagnetic proton-cyclotron waves propagating parallel to the background magnetic field.Further studies of the solar wind turbulence in the neighborhood of the break of the inertial range spectra as a function of wave vector, with components parallel (k ) and perpendicular (k ⊥ ) to the magnetic field, revealed that the fluctuation spectrum is anisotropic and that the power distribution sometimes is greater at quasi-perpendicular wave vectors (k ⊥ k ) (Matthaeus et al., 1990;Horbury et al., 2005;Dasso et al., 2005) than at quasi-parallel propagation (k k ⊥ ).However, at very long wavelength (smaller k), there is still enough energy available for the quasi-parallel propagating waves to dominate the oblique wave modes (where k k ⊥ ) (Matthaeus et al., 1990(Matthaeus et al., , 1996a,b;,b;Leamon et al., 2000;Smith et al., 2001Smith et al., , 2006;;Horbury et al., 2005Horbury et al., , 2008;;Bale et al., 2005Bale et al., , 2009)), and we focus our study on this part of the wave spectrum range.
The analysis presented follows two approaches.First, we have developed a quasilinear theory (Davidson and Ogden, 1975;Yoon, 1992;Yoon et al., 2003;Isenberg and Vasquez, 2007;Moya et al., 2011;Pavan et al., 2011) to understand the energy cascade that transfers wave energy from longer to shorter wave modes, with the subsequent energy transfer to the particles.This analysis corresponds to the first nonlinear order in Vlasov's formalism (Alexandrov et al., 1984;Krall and Trivelpiece, 1986).Second, we performed a one-dimensional hybrid simulation (Gary et al., 1997;Ofman et al., 2001;Araneda et al., 2002) of the system using an expanding box model (Velli et al., 1990;Liewer et al., 2001;Matteini et al., 2006;Hellinger and Trávničeck, 2008;Ofman et al., 2011) where a thin box of plasma moves away from the Sun in a moving frame at the local solar wind speed, in order to investigate the full nonlinear wave-particle interaction of the cascade process.All these effects; i.e., energy cascade, expansion, and nonlinear wave-particle interaction, are included in the study to show how the shape of the particle velocity distribution functions is controlled and regulated in kinetic plasmas.
This article is organized as follows.In Sect. 2 we show the basic equations of the quasilinear theory for the evolution of the macroscopic parameters of the distribution function, and we present numerical results for the case of the propagation of circularly polarized electromagnetic waves parallel to an ambient magnetic field, through a plasma composed of electrons and protons with thermal anisotropy.In Sect. 3 we present the equations and results for the same problem as in Sect.2, but with the use of unidimensional hybrid simulations with and without the inclusion of expansion effects.Finally, in Sect. 4 we compare the results of Sects. 2 and 3 and summarize the conclusions of this article.

Dispersion relation
We consider a plasma in an external magnetic field B 0 = B 0 ẑ composed of electrons and protons drifting with velocity V relative to a fixed frame (the "Lab" frame) along the background magnetic field.U = V /V Ap is the normalized drift speed of the protons, with V Ap = B 0 / 4πn p m p as the proton Alfvén velocity.m p and n p are the proton mass and density, respectively.We assume neutrality (i.e., zero net charge such that n e = n p ) and impose a zero-current condition along B 0 (U e = U , where U e is the drift velocity of the electrons).
The normalized dispersion relation for proton-cyclotron waves with left polarization, moving in the direction of the external magnetic field, in the case of a bi-Maxwellian distribution is (Gomberoff andValdivia, 2002, 2003;Gomberoff et al., 2004;Moya et al., 2011) where we have assumed that V Ap c.In Eq. ( 1), x y = ω k / p and y = ck/ω pp are the normalized complex frequency and wave number, respectively, with p = eB 0 /m p c and ω pp = (4π n p e 2 /m p ) 1/2 as the proton cyclotron and plasmas frequencies, respectively.e is the proton charge and c is the speed of light.Also, Z is the plasma dispersion function (Fried and Conte, 1961), ϕ y = x y − 1 − yU/yβ , β j = v th,j /V Ap where v th,j = 2K B T j /m p , with j = , ⊥ the parallel and perpendicular (with respect to the background magnetic field) thermal velocities of protons, respectively, and K B is the Boltzmann constant.Finally, we define A p = T ⊥ /T − 1 = β 2 ⊥ /β 2 − 1 as the proton thermal anisotropy.Typically, values of A p between 2 and 5 have been reported in the solar wind (Kivelson and Russell, 1995;Kohl et al., 1998;Cranmer, 2002Cranmer, , 2005;;Aschwanden, 2006;Kamide and Chian, 2007).
For now, we shall assume β ⊥ , β 1, such as in coronal holes (Gary, 2001;Aschwanden, 2006).Therefore, the argument of the Z function is much larger than one and we can use the semi-cold approximation (i.e., large argument asymptotic expansion) for protons and also consider electrons as cold (Gomberoff and Elgueta, 1991;Astudillo, 1996;Gomberoff andValdivia, 2002, 2003;Moya et al., 2011).We then write x y = x + iγ and assume |γ | |x|.Upon separation of real and imaginary parts of Eq. ( 2), we obtain the dispersion relation in the semi-cold regime.Thus, and where In Fig. 1 we show the branches of the dispersion relation Eq.
(3) for U = 0.3 in the normalized x-y space.Solid and dashed lines correspond to the two different solutions of the dispersion relation in the semi-cold regime.
Figure 2 shows the corresponding growth rates for the two solutions of the dispersion relation obtained from Eq. (3) for β 2 ⊥ = 0.05 = 5β 2 (A p = 4).The solid and dashed lines correspond to the solid and dashed curves in Fig. 1, respectively.It is observed that γ (y) is positive only in a narrow region near y = 1.5, which corresponds to an instability region.Also, for y > 1.5, the growth rate becomes negative and leads to a resonant absorption region.Depending on the magnetic field energy spectrum these two effects will compete and will change the shape of the distribution function.3) for U = 0.3.

Temporal evolution of macroscopic parameters
Assuming that the macroscopic parameters of the distribution function vary slowly (compared with the wave period) in time, the nonlinear temporal evolution equation of the proton distribution function f 0 , in the quasilinear approximation, is given by (Alexandrov et al., 1984;Krall and Trivelpiece, 1986;Yoon, 1992;Moya et al., 2011) where, in cylindrical coordinates is the first order perturbation of the distribution function and E k is the Fourier spectrum of the electric field.
Because of the bi-Maxwellian form of f 0 , we can define (Moya et al., 2011) and write a set of ordinary differential equations for all the parameters of the equilibrium distribution function, given by where τ = p t is the normalized time.Of course, at every time step, a solution of these equations will require us to solve the dispersion relation and integrate the functions K i (τ ) over k, because they depend on the parameters of the distribution function and the frequency ω k .Finally, to close the system of equations, we use the equation for the temporal evolution of the magnetic field spectral energy per unit length (Alexandrov et al., 1984;Krall and Trivelpiece, 1986;Moya et al., 2011) where with L as the reference length of the plasma.Thus, we have expressed the coupled system of Eqs. ( 11)-( 13) as an integro-differential system, where the temporal derivatives of the parameters correspond to k-space integrals.We note that in all the equations, we have both the frequency and growth rate, thus, to solve the system we need to explicitly solve x(y) and γ (x, y) from Eqs. (3) and (4).

Numerical results
In order to solve the system of quasilinear equations, we use a discrete grid in y-space with N y = 400 points for −2 < y < 2, hence, the separation between points is dy = 0.01.The origin has been avoided for obvious reasons.The time step is chosen to be dτ = 0.025.Knowing the magnetic field spectrum and the value of the parameters at a time τ we can solve the dispersion relation to obtain x and γ , as a function of y, at this particular time.Then, we can calculate the integrals defined in the system of equations , to evaluate the time derivative of each parameter.So, with that information, we use a 4th order Runge-Kutta method (García, 2000) to evolve the whole system to the next time step τ + dτ .Numerical integration is performed from τ = 0 until τ = 1000 using the Alfvén branch (solid line in Fig. 1).As initial values we use the same values of U p and β j used in Figs. 1 and  2, and we choose a Gaussian initial magnetic field spectrum The range in y-space was chosen in order to compare with the results of the hybrid simulations shown in the next section.The Gaussian profile (Eq.16) of the initial magnetic energy spectrum was chosen to concentrate most of the wave energy in low frequency waves, with negligible power for y > 1 (hence x < 1), e.g., below the proton resonance frequency as has been observed in the solar wind (Cranmer et al., 1999a,b).Even though the free energy is equally partitioned between waves and particles (Moya et al., 2011), it still has a non-negligible amount of the free energy in the γ > 0 region as in Yoon (1992).Due to the shape of γ (x, y) of the solid curve in Fig. 2, only a small portion of the initial waves will effectively interact with the particles.Furthermore, due to the shape of K 1 (τ ) in Eq. ( 11), in every time step the total time derivative of U is essentially zero, so that the proton drift speed can be considered essentially a constant.This is due to the fact that the contribution of the wave momentum to the total momentum is of the order of (V Ap /c) 2 1.Thus, for low frequency waves, in this approximation, there is no transfer of momentum between particles and waves.It is important to mention that in the quasilinear approximation, the whole temporal dependence of the macroscopic parameters of the distribution function depend strongly on the imaginary part of the frequency.In the case of the chosen parameters γ ∼ 10 −2 , thus integrating until t ∼ 10 3 is long enough to observe significant quasilinear effects, but not to consider higher order nonlinear effects, such as mode coupling, etc.
The temporal evolution of the two temperatures of the protons, with respect to their initial values, is shown in Fig. 3.We see evidence of a parallel heating (dashed curve in Fig. 3) and a perpendicular cooling (solid curve in Fig. 3) of the plasma.For large times the heating saturates and it seems that the system is slowly approaching a metastable situation where T = T ⊥ .At the end of the integrated interval, τ = 1000, the parallel temperature T is T ∼ 1.4T (0) and the perpendicular temperature T ⊥ exhibits a small cooling of no more than 5 %.Thereby, owing to the changes in both temperatures, the anisotropy also evolves as time goes on, as shown in Fig. 4. As time progresses, A p (τ ) decreases but, due to the slowness of the process, beyond τ ∼ 1000 it is not possible to draw conclusions about the complete behavior of the anisotropy for longer times.Towards the end of the numerical integration, A p decays to a final value of A p ∼ 2.5.
Finally, to close the quasilinear system of equations, the total free energy available due to the thermal anisotropy of the proton distribution function is absorbed by the electromagnetic field to preserve the total energy of the system (during the integration until τ = 1000 the total energy of the system remained constant to better than 0.1 %).In Fig. 5 we show an amplified portion of the magnetic field energy spectrum profile as a function of y near y = 1.3 (Fig. 5a), and as a function of x near x = 1.1 (Fig. 5b).The initial Gaussian spectrum changes due to the presence of regions of absorption and instabilities (see Fig. 2).Even though the absorption zone is greater than the instability zone at τ = 0, there is little energy for y > 2 (x > 1.4) so that the waves cannot transfer a noticeable amount of energy to the protons.Thus, the net effect of the propagation of waves through the electron-proton plasma is the generation of instabilities for 1 < y < 1.6 (0.9 < x < 1.3).In Fig. 5 it is shown that the particle energy loss is transfered to the electromagnetic field as an emergence of wave modes of higher energy density than the initial ones.Thus, for the initial set of parameters in this quasilinear approximation, the whole system evolves by transferring energy from the particles to the waves, slowly approaching a metastable equilibrium.

Hybrid model
Due to the large difference in mass between electrons and protons, an efficient, and very common way to model a plasma system is to consider electrons as a massless fluid.Thus, we can consider protons as kinetic particles with a certain velocity distribution function and electrons as a charged fluid with bulk velocity U e .As a fluid, the temporal evolution for U e is given by the momentum equation where m e is the electron mass and P e = n e k B T e is the pressure of the electron fluid.Here, T e is the electron temperature and we have assumed the quasi-neutrality condition n e ≈ n p ≈ n, where n is the average density of the plasma.
In addition to the evolution of particles, we are also interested in following the space and time evolution of the electromagnetic field.This temporal evolution is given by Maxwell's equations: where ρ and J are the charge and current densities, respectively.Since we are interested in low frequency waves, it can be shown that the displacement current term in Ampere's equation is of order O((V Ap /c) 2 ) and it can be neglected.Thus, we can deduce an expression for the current in terms of the curl of the magnetic field Similarly, the current density of a plasma with electrons and ions is given by the vectorial sum of proton and electron currents J = en(U p − U e ).Here, U p is the fluid bulk velocity of protons.Thus, neglecting the LHS of Eq. ( 17), solving the equation for E, and using Maxwell's equations (Eqs.19), we obtain a set of equations for the hybrid model: where x i and v i are the position and velocity of the i-th proton and, in this approximation, electron temperature T e remains constant and equal to zero.

Hybrid expanding box model
Numerical simulation of plasmas in intrinsically spherical geometry has always been complicated due to the problem of self-forces and the selection of the shape function for the macroparticles (Dawson, 1983;Birdsall and Langdon, 1985;Matsumoto and Omura, 1993) Velli et al., 1990;Grappin and Velli, 1996;Liewer et al., 2001).The radial position of the box, relative to the origin of S is given by R = R(t) = R 0 + U 0 t , where R 0 x is the position of the box at t = 0.Then, to transform coordinates to a frame moving with the box (the S frame with coordinates x , y , z ) a two step transformation process is done.First, a Galilean shift between x and x , and then a stretching in y and z.Namely, where a = 1 + U 0 t/R 0 .These transformations imply that, for an observer moving within the box, the box does not change its volume, but for an observer at rest with respect to the S frame, the box is expanding as it moves away from the origin.Then, using these transformations ( 25)-( 27) it is possible to represent the time and the spatial derivatives in the S frame, and then to express Eqs.( 21)-( 24) for the hybrid model in the expanding box frame (Liewer et al., 2001) as where the quantities are measured in the moving frame.The derivative operators are given by and P = diag(0, 1, 1), L = diag(2, 1, 1) and A(t) = diag(1, 1/a, 1/a) are diagonal matrices.Also, a can be written as a = 1 + τ , where τ is the time normalized to the proton cyclotron frequency p , and = U 0 /(R 0 p ).

Numerical results
As a first approximation to the problem, and to compare with the quasilinear results, we have performed a standard numerical hybrid simulation in one spatial dimension (the x coordinate) and three dimensions in velocity space, in the presence of a background magnetic field B 0 = B 0 x.To solve the fields equation, we use a periodic system in x, with N x = 512 cells of size dx = 0.5 and 1000 particles per cell, where we have chosen the same spatial normalization as in the quasilinear case, namely the ion inertial length λ p = V Ap / p .The particle and field equations are integrated in time using a rational Runge-Kutta method (Ofman et al., 2001(Ofman et al., , 2011)), whereas the spatial derivatives are calculated using a pseudo-spectral fast Fourier transform method.The time step used was dτ = 0.025 and the simulation was performed until τ = 1000.
In the case of standard hybrid simulations, the natural choice for a reference frame is the electron frame.Thus, to keep quasi-neutrality and parallel zero current, the initial condition for proton drift is U = 0.The velocity distribution function is initialized as a bi-Maxwellian using a random number generator and, in order to compare the quasilinear and hybrid models, we choose the same initial values for both proton temperatures β 2 ⊥ = 0.05 = 5β 2 and T e = 0 for electrons.The parallel x component of the magnetic field remains constant, equal to B 0 , throughout the simulations.For the initial parameters mentioned, we have performed simulations with an initial Gaussian perpendicular magnetic field with random phases, where, in order to compare, the amplitude and exponent of the Gaussian was chosen to have the same magnetic spectrum as in the quasilinear method (Eq.16).
In Fig. 6a, we show the time dependence of perpendicular and parallel temperatures during the simulation for the hybrid models.Comparing with the quasilinear method, the figure shows that in the case of the standard hybrid simulation, parallel heating (dashed curve in Fig. 6) is slower, and perpendicular cooling (solid curve in the figure) is similar to the quasilinear solution, due to the fully nonlinear resonant absorption present in hybrid models.Thus, as it can be seen in Fig. 6b, the final values obtained correspond to T ⊥ ∼ 0.95T ⊥ (0) and T ∼ 1.1T (0) for perpendicular and parallel temperatures, respectively.However, the combined effects produce a final, nonzero value (A p ∼ 3.4) of the thermal anisotropy.Compared to the quasilinear case, as shown in Fig. 7, it is observed that, although the time evolution of both curves does not match, the slopes of both curves are similar.Thus, qualitatively the two approaches agree.Also, like in the quasilinear approach, in our simulations the fluid drift U had no significant changes.It is important to mention that it is expected that both models differ.In the quasilinear approach, we follow just one solution (the Alfvén branch) of the dispersion relation and the nonlinear terms are approximated by the quasilinear theory.On the other hand, in hybrid models, the simulations consider all the branches of the dispersion relation, including the instabilities propagating antiparallel to the background magnetic field.Also, being completely nonlinear, hybrid simulations include several nonlinear effects as coupling (wave-wave interaction) between all the modes allowed by the dispersion relation.
On the other hand, for the case of the expanding box hybrid simulations, Fig. 6 shows how the expansion induces an increase in the perpendicular cooling of the plasma.While there are no significant differences between both models in the time dependence of parallel temperatures (dashed and dotted curves in Fig. 6a), the expansion hastens the decline in this quantity as time evolves in the case of perpendicular temperature (dot-dashed curve).Also, in Fig. 6b we show both temperatures with respect to their initial values.The figure shows that parallel temperature increases by 5 % while the perpendicular temperature decreases approximately to a value of T ⊥ ∼ 0.8T ⊥ (0).Thus, the thermal anisotropy in the expansion model relaxes faster than in the standard model without expansion, as it can be seen in Fig. 7, reaching a value of A p (τ = 1000) ∼ 2.9.Also, as the expansion is regulated by a(τ ) = 1 + τ (in our simulations we chose = 10 −4 ) from the figure it is clear that notable effects only occur at long enough times τ ∼ 200.It is important to mention that in the case of the solar corona and solar wind, that correspond to spherically expanding plasmas, ∼ 10 −5 .However, if we can see qualitative differences between the simple and the expanding box hybrid models at times of order τ ∼ 10 3 we need to amplify the parameter to be able to draw con- clusions about the differences and similarities between both methods (Liewer et al., 2001).
The acceleration in the relaxation of the temperatures, due to the spherical expansion of the plasma, also influences the behavior of electromagnetic waves.In Fig. 8 we show the time evolution of the perpendicular magnetic field energy until the end of the simulation.Owing to the interaction between protons and waves, for times below τ ∼ 200, the average perpendicular magnetic field energy increases monotonically.For larger times, the influence on the expanding box parameter a begins to be relevant, so that there is a considerable decrease in the perpendicular magnetic field energy compared to the case of the standard hybrid simulation (solid curve in the figure).
Starting from two different kinetic approaches, the first theoretical based on quasilinear theory and the second based on computational simulations of hybrid models, we have done numerical studies on the interaction of waves and protons.
For the chosen parameters, our results indicate that, due to the shape of the K j functions in our quasilinear method, the main change occured in the parallel temperature, while in the case of the hybrid code without expansion, the main effect of the interaction between particles and waves was a decrease of the perpendicular temperatures.However, the quasilinear and the standard hybrid simulations agree on the macroscopic evolution of the drift velocity parallel to the ambient magnetic field, and also agree on the evolution of the slope of the thermal anisotropy, even when both approaches differ in the absolute time profile of the temperatures and anisotropies.
We also performed numerical simulations using an expanding box hybrid model using the same initial parameters as in the quasilinear and standard hybrid models.Our results show how the expansion produces a relaxation of the system in the perpendicular plane.This model shows a decrease of the average magnetic field energy and also the typical cooling of an expanding gas.Nevertheless, due to the a = 1 + τ parameter, those effects are significant only for long enough times.
In conclusion, our results allow us to obtain and compare the basic properties of the wave-particle interaction, in simple electron-proton plasmas, using different models.Numerical results show that both quasilinear and nonlinear methods qualitatively agree in the evolution of the macroscopic plasma parameters, and this seems to suggest that resonant absorption and the energy cascade mentioned above may be relevant in the heating of solar wind plasma.

Fig. 1 .
Fig. 1.Normalized dispersion relation branches.Solid and dashed curves correspond to the two solutions of Eq. (3) for U = 0.3.

Fig. 2 .
Fig. 2. Normalized growth rates of the two solutions of the dispersion relation for β 2 ⊥ = 0.05 = 5β 2 (A p = 4).Solid and dashed lines correspond to solid and dashed curves in Fig. 1, respectively.

Fig. 3 .
Fig. 3. Temporal evolution of both proton temperatures with respect to their initial values.Solid and dashed lines correspond to perpendicular and parallel temperatures, respectively.

Fig. 5 .
Fig. 5. Normalized magnetic field energy spectrum ε for τ = 250 (solid), τ = 500 (dashed), τ = 750 (dotted) and τ = 1000 (dot dashed).It is observed that as time advances, there is an emergence of waves at higher modes than what was originally available.(a) ε as a function of the wave number y.(b) ε as a function of the frequency x.
Fig. 6.(a) Normalized proton temperatures as a function of time for both hybrid simulations (with and without expansion).Dotted and dot-dashed curves correspond to parallel and perpendicular temperatures in the expanding box model simulation (with = 10 −4 ), respectively.On the other hand, dashed and solid curves correspond to parallel and perpendicular temperatures for simple hybrid simulations ( = 0).(b) The same evolution shown in panel (a), but with respect to the initial values of both temperatures for the simple and expanding box simulations.

Fig. 7 .Fig. 8 .
Fig. 7.Proton thermal anisotropy for standard = 0 (solid) and expanding box = 10 −4 (dashed) hybrid simulations.Here it can be seen that the expansion effects are notable only from long enough times during the simulation.In order to compare, the anisotropy for the quasilinear solution is also included (dotted curve).

2012 1366 P. S. Moya et al.: Computational and theoretical study of the wave-particle interaction coordinates
. Recently, a formalism in which a plasma in spherical expansion, like the solar wind, has been modeled in a Cartesian system.In this expanding box model, a thin shell of plasma (the box) is considered moving radially away from the origin of certain frame S with www.ann-geophys.net/30/1361/2012/Ann.Geophys., 30, 1361-1369, x , y , z at a constant velocity U 0 x (