Particle simulation of plasma heating by a large-amplitude shear Alfvén wave through its transverse modulation in collisionless plasmas

We investigate the nonlinear dynamics of a large-amplitude shear Alfvén wave (AW) with circular polarization in collisionless plasmas by using a 2D3V fully relativistic electromagnetic particle-in-cell (PIC) simulation code. We found that when the amplitude, A = δB/B0, of the shear AW is larger than , the shear AW spontaneously becomes unstable for the modified two-stream instability (MTSI), resulting in the excitation of small-scale quasi-electrostatic waves with electric fields parallel to a uniform magnetic field. We found that the electrons are heated mostly in the direction parallel to the magnetic field due to the quasi-electrostatic lower hybrid waves by the MTSI. Subsequently, the shear AW becomes unstable with strong transverse modulation, resulting in the excitation of ion acoustic waves that can heat ions. About 70% of the AW energy can be converted to plasma heating of both electrons and ions. The plasma heating time is within about 250ωci−1, which is shorter than the collision time between protons and neutral hydrogens in the upper chromosphere. The obtained results could be very important for plasma heating of coronal loops in the upper chromosphere.


Introduction
Alfvén, in 1942 [1], found that a distinctive wave mode arises in a highly conducting, magnetized and incompressible plasma, the so-called 'shear Alfvén wave (AW)' that propagates along the magnetic field direction. The existence of the wave was experimentally verified by Lundquist [2]. The importance of the AWs was soon realized and since then a lot of studies have been done for space, astrophysical and laboratory plasmas (see [3,4]).
In solar plasma physics there remains an unsolved problem of how the solar corona is heated. The shear AWs that may be excited near the photosphere through convective motions or the interaction [5] between magnetic flux tubes or the collision [6] between magnetic flux tube and magnetosonic waves are believed to be an important energy carrier from the photospheric kinetic energy. A number of works have been done to study how the shear AWs can be damped effectively in the upper chromosphere or lower corona. Among them collisional damping as well as collisionless ion cyclotron damping [7] were proposed. The AWs with a frequency less than about 0.6ω ci are hardly damped through the ion cyclotron damping. Therefore, it is very important to investigate the effective damping mechanism(s) of finite amplitude shear AWs with a frequency less than 0.6ω ci in collisionless plasmas.
Linear resonant damping of surface AWs along the loop [8] is also investigated, associated with mode conversion into the kineticAW which has a longitudinal component of electric field and a short perpendicular wavelength. The AWs propagating in a plasma with density nonuniformity transverse to the magnetic field can damp through the phase-mixing effect [9]. Previous works [10]- [12] of the phase-mixing were studied in the MHD approximation. Recently, Tsiklauri et al [13,14] investigated the phase-mixing of shear AWs in collisionless plasmas by using a 2D3V fully relativistic electromagnetic particle-in-cell (PIC) simulation code. They found that the electrons can be effectively heated through the electric field generated in the region of density inhomogeneity. Therefore, when we study plasma heating by the wave-particle interaction of shear AWs, we must take into account full electron dynamics as well as ion dynamics, beyond the MHD approximation.
The nonlinear dynamics of the shear AWs was investigated mostly in the framework of the MHD approximation; for examples, parametric or decay instability [15,16], soliton dynamics [17,18], parametric instability by using a hybrid simulation code [19,20], self-focusing dynamics [21] through the poderomotive force, triggering of a tearing mode by AWs [22] and transverse 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT modulation including the Hall effect [23]. Furthermore, in weakly ionized plasmas, transverse modulation of shear AWs was investigated [24].
In this paper, we investigate the nonlinear dynamics of a large-amplitude shear AW with circular polarization in collisionless plasmas by using a 2D3V fully relativistic electromagnetic PIC simulation code. We found that when the amplitude, A = δB/B 0 , of the shear AW is larger than A c = v th,i c ω pe kc ω pe ω ce , the shear AW spontaneously becomes unstable for the modified two-stream instability (MTSI) [25], resulting in the excitation of small-scale quasi-electrostatic waves with electric fields parallel to a uniform magnetic field. We found that the electrons are heated mostly in the direction parallel to the magnetic field due to the quasi-electrostatic waves by the MTSI. Subsequently, the shear AW becomes unstable with strong transverse modulation, resulting in the excitation of ion acoustic waves that can heat ions. About 70% of the AW energy can be converted to plasma heating of both electrons and ions. The plasma heating time is within about 250ω ci −1 , which is shorter than the collision time between protons and neutral hydrogens in the upper chromosphere. The obtained results could be very important for plasma heating of coronal loops in the upper chromosphere.
The paper is organized as follows. In section 2 we present our simulation model, in section 3 we show our simulation results and in the final section we discuss heating of chromospheric plasma based on the present results and provide a summary.

Simulation model
We used a 2D3V, fully relativistic electromagnetic PIC code, modified from the 3D3V TRISTAN code [26]. The system size is L x = 1055 and L y = 485 , where (= 1) is the grid size. The periodic boundary conditions for both x-and y-directions are imposed on particles and fields. The ratio of ion mass to the electron mass is m i /m e = 16, the simulation time step is ω ci t = 0.003125 (ω ci is the ion cyclotron frequency), respectively. The external uniform magnetic field, B 0 , is in the x-direction and the Alfvén velocity v A = 0.25c. The ratio of electron cyclotron frequency to the plasma frequency is about ω ce /ω pe = 1, the Debye length is v th,e /ω pe = 1 , plasma beta is β = 0.02, skin depth is c/ω pe = 10.0 , while the ion skin depth is c/ω pi = 40.0 . Here ω pi is the ion plasma frequency. The electron thermal velocity is v th,e = 0.1c, electron Larmor radius is ρ e = v th,e /ω ce = 1 , the ion thermal velocity is v th,i = 0.025c and ion Larmor radius is ρ i = v th,i /ω ci = 4 .
To excite the AW with left-hand circular polarization in the simulation domain, we impose the initial magnetic and electric fields as follows: where CGS Gaussian units are being used. We choose the wavenumber of the AW as kc/ω pe = 0.179, which means that the AW with three wavelengths is excited in the simulation domain. The frequency ω of the AW is determined from the dispersion relation of cold plasma approximation, taking into account both electron and ion dynamics. Then we have the frequency of AW as ω = 0.49ω ci = 3 × 10 −2 ω pe . The drift velocities associated with the above AW for both ions and electrons are calculated from the equations of motions coupled with Maxwell's equations 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT as follows: Then the above drift velocities are added as the shifted Maxwell velocity distribution function. We checked that the small-amplitude AW with frequency larger than about 0.6ω ci damps through the ion cyclotron damping, whose damping rate agrees well with the theoretical value. We also confirmed that when the amplitude of the AW with frequency larger than about 0.6ω ci increases, there occurs the amplitude oscillation due to the phase-trapping of ions, resulting in very weak ion cyclotron damping. In the next section, we present mostly the simulation results of the AW amplitude with δB/B 0 = 0.5 and a frequency of 0.5ω ci . The accuracy of the present simulation is confirmed by which the total energy in the system is conserved within 0.08%.

Simulation results
The time evolution of finite-amplitude shear AW is characterized by three phases (I)-(III) shown in figure 1: (I) the phase of excitation of MTSI with electron heating, (II) the phase of transverse modulation of the AW and (III) the phase of ion heating due to ion acoustic waves. Figure 1(a) shows time evolution of magnetic field energies associated with shear AW in the system: B y 2 2 dx dy, B z 2 2 dx dy. Figure 1(b) shows time evolution of electric field energies associated with AW in the system E y 2 2 dx dy and E z 2 2 dx dy, while as seen in the figure the field energy E x 2 2 dx dy, which is initially zero, grows and attains to a maximum value at around 30ω ci t shown by the arrow (I). The electric field energy As seen in the phase (III) of figure 1(d), ions are gradually heated during the decay phase of transverse modulation. As seen later we found that during the transverse modulation of the shear AW, ion acoustic waves are excited and they contribute to ion heating.
The energy conversion rate from the AW energy to plasma heating is about 70% and the time for the plasma heating is about 2.5 × 10 2 ω ci −1 s.

Phase (I): excitation of MTSI with electron heating
Firstly, we study the phase (I).  is a large-scale wave pattern corresponding to the shear AW and another is a small-scale structure similar to the pattern of E x in figure 3. The characteristic scale-length in the y-direction of E x and E y is about 2.5c/ω pe = 25 . To understand the characteristics of the excited waves with small spatial structures, we performed two-dimensional Fourier transformation by using the space-time data of the electric field E x for the time period between 0 and 16384 steps both for x-and y-directions to obtain the wave dispersion relation. Figures 5(a) and (b) are dispersion relations of E x for parallel and perpendicular direction to the magnetic field, respectively. From figure 5(b) we find that the waves propagating perpendicular to the magnetic field are excited with the wavenumber k y c/ω pe = 2.5, which corresponds to the wavelength observed by the small-scale structures.
To understand the origin of the above observed quasi-electrostatic waves, we must consider the ion drift velocity U ⊥ to excite the shear AW that is approximately given by U ⊥ /c = A k x c ω pe ω ce ω pe , where A = δB/B 0 is the amplitude of the shearAW. If the shearAW has a amplitude larger than the critical value A c = v th,i c ω pe kc ω pe ω ce , the ion drift velocity U ⊥ required to excite the shear AW becomes larger than the v th,i . This condition is satisfied for the present simulation with A = 0.5, ω ce = ω pe  and k x c ω pe = 0.179. The ion drift velocity U ⊥ is larger than the ion thermal velocity v th,i = 0.025c, as seen in the figure of the initial ion velocity distribution function in the y-and z-directions (see figure 10). This situation is known to be unstable for the MTSI [25].  figure 5(a) and it is about 0.155c that is larger than the electron thermal velocity v th,e = 0.1c.  As seen in figure 5(a), the excited waves have the frequency in the range of ω ci < Re[ω] < ω ce . Therefore, the excited waves are lower hybrid waves. The propagation angle of the wave is about 60 • . Therefore, we conclude that the quasi-electrostatic lower hybrid waves excited during the phase (I) are due to the MTSI caused by the ion transverse drift associated with the shear AW. As seen in the phase (I) of figure 1(c), the excited lower hybrid waves associated with the MTSI are damped, resulting in a strong increase in the electron kinetic energy K ex in the x-direction.
Here we comment on the mass ratio used. The lower hybrid waves excited by the MTSI almost propagate perpendicular to the magnetic field, if we use the real mass ratio, m i /m e = 1836. The parameter (ω ce /ω pe = 3) used by Ott et al [25] is different from our present one (ω ce /ω pe = 1). Ott et al considered that the ions are unmagnetized. This leads to perpendicular ion heating. But in our case the growth rate of MTSI becomes comparable to the AW frequency; therefore, ions are magnetized and perpendicular ion heating is suppressed during MTSI. However, if we use the real mass ratio, the growth rate of MTSI becomes larger than the AW frequency. Again the ions are unmagnetized during MTSI. Thus, we need to further study this effect.

Phase (II): transverse modulation of the AW
In this subsection we consider the second phase (II) that is indicated in figure 1(b). This phase is characterized by the strong damped oscillation of wave electric and magnetic fields, as  Figure 7 shows time evolution of magnetic field B y at (a) ω ci t = 125 and (b) ω ci t = 156.25 and electric field E x at (c) ω ci t = 125 and (d) ω ci t = 156.25 during the late phase of transverse modulation. In this phase there appear second harmonic waves associated with the electric field E x due to strong nonlinearity, as seen in figures 6(c) and (d).

Phase (III): ion heating due to ion acoustic waves
In the last phase (III) that is shown in figure 1(d), the damped transverse modulation of the shear AW changes from strong to moderate amplitude oscillation and here there occurs an increase in ion kinetic energy K ix . To understand the characteristics of the electric field E x during the phases (II) and (III), we performed a two-dimensional Fourier analysis of E x for the case of the AW   figure 8. The excited wave is the ion acoustic waves that eventually damp to heat ions. Therefore, we conclude that the decrease in the field energy E x from about ω ci t = 150 seen in figure 1(b) is due to the damping of ion acoustic waves, resulting in the increase in ion kinetic energy in the phase (III).
As seen in the time evolution of the magnetic and electric field energies of the shear AW, their period becomes longer than the initial one, and finally about 70% of the shear AW energy can be converted into plasma kinetic energy. In figure 9 the time evolution of the phase-space plot of both electron and ions is shown at three time steps, at (a) ω ci t = 0, (b) ω ci t = 31.25 and (c) ω ci t = 156.25. The left column shows the electron phase-space, while the right column shows the ion phase-space. As seen in the electron phase-space at ω ci t = 31.25 that corresponds to the phase (I), the electron heating occurs almost uniformly in the x-direction, because the excited waves associated with the MTSI have relatively broad band spectrum in the k x direction. However, as seen in the electron phase-space at ω ci t = 156.25 that corresponds to the late phase (II), the electron heating is not uniform in the x-direction; it strongly couples with the longitudinal modulation of the electric field E x . On the other hand, the ions are heated with slow modulation in the x-direction as seen in figure 9(c). Finally in figure 10, we show the time evolution of the velocity distribution function of both electrons and ions at three time steps; ω ci t = 0, 31.25 and 250. The left column shows the three components of electron velocity distribution functions, while the right column shows the three components of ion velocity distribution functions. As mentioned before, the ion drift velocity U ⊥ is larger than the ion thermal velocity v th,i = 0.025c, as seen in the figure of the initial ion velocity distribution function in the y-and z-directions (see figure 10). This is the reason why the MTSI can generate quasi-electrostatic waves, resulting in the electron heating in slight excess in the positive x-direction. Final ion velocity distribution in the x-direction also shows a slight excess of high-energy particles in the x-direction. Therefore, we expect that due to the damping of the shear AW, both ions and electrons have bulk flows along the magnetic field. This effect could be important for the solar wind acceleration in the upper chromosphere. But we need to further explore the physical conditions that lead to the solar wind acceleration. We note here how the amplitude of shear AW changes the present simulation results. To see this effect, we did simulations with different wave amplitudes of δB/B 0 = 0.2, 0.3, 0.5 and 1.5. Even in the case of δB/B 0 = 0.2, we observed the same transverse modulation shown here. The only difference is that when the wave amplitude becomes small, the onset of the excitation of the MTSI becomes slow and its growth rate decreases. We show the comparison of the onset of the excitation of the MTSI in figure 11. As seen in figure 11(2 and 3), there does not occur a strong hump in the time evolution of electric field energy E x .

Discussions and summary
New observations [27] of coronal loops in EUV wavelengths with the Transition Region and Coronal Explorer (TRACE) demonstrated that the most influential model of coronal loops originated from Rosner et al (1978) [28], assuming that a constant pressure along the loop is inapplicable to a class of EUV loops. Those loops are found to be heated near the foot-points, with a heating scale height of 12 ± 5 mm. These results support coronal heating mechanisms operating in or near the chromosphere and transition region. We apply the simulation results obtained in the previous section to the upper chromosphere. The parameter used in the present simulation, v A > v th,e , is satisfied if the electron temperature is T e < 6 × 10 4 K, where we used the magnetic field B 0 = 100 G and the electron number density n e = 5 × 10 10 cm −3 , considering that the chromospheric height from the photosphere is about 2000 km. The protonneutral hydrogen collision frequency ν pn in this region [29] is about ν pn = 100 s −1 , which is slow compared with the collisionless damping frequency ν damp of the shear AW that is about ν damp = ω ci /250 = 4 × 10 3 s −1 (for B 0 = 100 G) as shown in the previous section. Therefore, the collisionless damping of shear AW associated with its transverse modulation is very important to heat the chromospheric plasma. As seen in the simulation, the energy transfer rate from the shear AW to plasmas is very efficient and rapid. The collisionless damping process of shear AWs found here may play an important role for coronal loop heating near the foot-points.