Ion acceleration and heating by kinetic Alfvén waves associated with magnetic reconnection

Magnetic reconnection, which plays a key role in the conversion of magnetic energy to thermal and kinetic energy of particles, is fundamentally important in magnetized solar, geophysical, and laboratory plasmas. The onset of magnetic reconnection has been observed at multiple scales, such as solar flares in solar corona, substorms in the magnetosphere, Earth's dipole field, interplanetary surroundings,^2,3 and tokamak plasmas.³⁸ Moreover, perturbations associated with reconnection can span an even larger range of scales due to nonlinear processes. When fluctuations occur on kinetic scales they can efficiently release the free energy available in the magnetic configuration into particle acceleration and heating. Understanding the processes responsible for acceleration and heating is of key importance in the study of reconnection.

Many experimental and numerical studies have been performed to examine the mechanisms of the acceleration of particles and the conversion of energy during reconnection. Based on the experiments on toroidal plasma merging, Ono et al.⁴ have shown that ions can be accelerated toroidally up to the order of the Alfvén velocity through contraction of the reconnected field lines. Meanwhile, ions are directly heated by conversion of up to 80% (⁠$±20%$⁠) of the dissipated magnetic energy into ion thermal energy. Through a two dimensional (2-D) Magnetic Reconnection Experiment (MRX), Yoo et al.⁵ have investigated ion acceleration and heating within the ion diffusion region in collisionless plasmas. They have suggested that the ions are accelerated near the separatrices to up to $0.5vA$ by the in-plane electric field before they are heated downstream inside the ion diffusion region, where heating may result from the local ion remagnetization, local frictional drag due to the high density, or ion kinetic effects. A numerical study on the energy transfer in the ion decoupling region has also been performed.⁶ It has been shown that in the ion diffusion region the magnetic energy converted to ion thermal energy is 3–5 times larger than that converted to the bulk kinetic energy of the ions. Using 2-D particle-in-cell (PIC) simulations and test particle simulations, Drake et al.⁷ have examined the ion heating downstream of the reconnection X-line. The trajectories of the test particles indicate that they gain an effective thermal speed as well as an outflow speed upon entering into the Alfvénic exhaust. Moreover, Cheng et al.⁸ have studied ion dynamics during collisionless driven magnetic reconnection using a 2-D PIC model. They have found that the ions become unmagnetized in the reconnection current layer, where the magnetic field reverses and ion motions are dominated by an orbit meandering effect, and around the separatrix region where ions flow from upstream to downstream. In these two regions, ions gain both kinetic and thermal energy mainly from the out-of-plane inductive electric field although ion dynamics are not the same.

Investigations have been presented to study the ion heating through resonant cyclotron interactions with high frequency Alfvén waves,⁹ which could interact with ions via the cyclotron resonance condition $ω−k∥v∥=±Ωi$⁠, according to the linear theory, in which ω is the wave frequency, $k∥$ is the parallel wave number, $v∥$ is the ion velocity component parallel to the ambient magnetic field, Ω_i is the ion cyclotron frequency, and $+/−$ signs represent the left/right hand polarizations. They have also found a strong anisotropy of ions with $T⊥/T∥≈30∼40$⁠. The resonant interactions with ion-cyclotron waves can also lead to the ion acceleration and heating. Gary et al.¹⁰ have studied the wave-ion scattering effects between the Alfvénic cyclotron waves and ions in the solar wind based on hybrid simulations, finding that ions are anisotropically heated with $T⊥>T∥$ while weakly accelerated. An experimental study on ion heating by ion-cyclotron waves driven by velocity shear has been carried out,¹¹ in which perpendicular ion heating and acceleration are also reported. In addition, stochastic ion heating and acceleration by low frequency drift-Alfvén waves have also been observed when the ion displacement due to polarization drift becomes comparable to the wavelength of the mode.¹² 

With large amplitude obliquely propagating Alfvén waves, significant perpendicular stochastic heating has been found when the wave frequency is a fraction of the cyclotron frequency.^13,14 Furthermore, with a spectrum including multiple modes, the stochastic heating can occur even if the amplitudes are significantly small.¹⁵ They have also proposed that such a heating mechanism plays a key role in the ion heating of solar corona and acceleration of solar wind. In addition, for the low frequency waves (⁠$ω<Ωi$⁠) with small amplitudes, perpendicular ion heating may also result from the higher order electric field terms.¹⁶ 

Kinetic Alfvén waves (KAWs), which is the extension of shear Alfvén wave branch in the regime when their perpendicular wavelength is comparable with the ion Larmor radius, also lead to ion heating. The ion heating and dissipation of the low frequency KAW with an amplitude larger than a threshold have been investigated.^14,17 It is found that heating is primarily in the perpendicular direction under a small β. Chandran et al.¹⁸ have examined ion heating by turbulent Alfvén waves and KAWs using test particle simulations. They have found that the heating is anisotropic with $δv⊥2≫δv∥2$ when $β≪1$⁠. In addition, they have also noticed that Landau damping and transit-time damping of KAWs would also lead to strong parallel heating when $β≥1$⁠. Later, Chaston, Bonnell, and Salem¹⁹ have estimated the heating rate both for the parallel and perpendicular directions through observational analysis. However, they have demonstrated that the entire heating is almost in the perpendicular direction. Moreover, the parallel electric field associated with KAWs would also be expected to lead to the particle energization and parallel heating. Nevertheless, it is not common to see the direct heating by $E∥$ if the amplitude of $E∥$ is small.¹⁴ 

KAWs are believed to be of great importance in reconnection since they on one hand, can break the Alfvén speed limit, and on the other, can carry sufficient parallel Poynting flux to drive transport through the diffusion region.^20–24 Our previous work has found that the KAWs are a common feature during magnetic reconnection based on 3-D hybrid simulations under various guide field strengths.¹ However, ion acceleration and heating by the resultant KAWs still remain unclear. In this paper, we examine the ion acceleration, heating, and trapping by the KAWs based on the 3-D hybrid current sheet model, which was used to study the KAWs associated with reconnection.¹ The paper is outlined as follows. In Sec. II, we describe the simulation model. In Sec. III, the simulation results are shown and discussed. Finally, a summary is given in Sec. IV.

In the hybrid model, the ions are regarded as fully kinetic particles while the electrons are treated as a massless fluid.^25–27 The model used in this paper is identical to the model that we used for the generation of KAWs in 3-D reconnection,¹ in which a Cartesian coordinate system is applied.

In the initial current sheet, the magnetic field components are B_x = 0, $By=By0$⁠, and

in which x is the current sheet normal direction, y is the current direction, z is the direction of the anti-parallel magnetic component, $By0$ is the initial guide field, $Bz0$ is the anti-parallel component in the asymptotic region, and δ is the half-width of the initial current sheets. In the simulation, there are two current sheets located at $x=Lx/4=16di$ and $x=3Lx/4=48di$⁠, respectively, where L_x describes the domain length in x.

The ions are initialized as a drifting-Maxwellian distribution. The thermal and magnetic pressures are balanced through the current sheets

Here, the left side of Equation (2) represents the local total pressure while the right side is the background (asymptotic) pressure, and normalization factor $α=4πe2/mic2$ in the simulation units. Thus, the initial ion number density can be derived as $N(x)=N0[1+1/β0(1−By(x)2/B02−Bz(x)2/B02)]$⁠, where β₀ is the total plasma β in the asymptotic region, N₀ is the ion number density in the asymptotic region, and $B0=(Bz02+By02)1/2$ is the asymptotic magnetic field.

In this hybrid code, the ions are updated by the ion equation of motion

where $v$ is the ion particle velocity, $E$ is the electric field, $B$ is the magnetic field, $u$ is the velocity of the ion bulk flow, u_e is the electron bulk flow velocity, and $ν=νJ+νc$ is the collisional frequency introduced corresponding to an ad-hoc resistivity. There are two parts of ν, with ν_J being a current-dependent resistivity due to the spontaneous anomalous resistivity²⁸ and ν_c being a localized resistivity in association with reconnection triggered by a local enhancement.¹ The localized ν_c can be expressed as

in which $x1/4=x−Lx/4, x3/4=x−3Lx/4, z1/2=z−Lz/2$ with L_z being the domain length in z, $y1/2=y−Ly/2$ with L_y being the domain length in y, λ₀ is the scale length of the resistivity in x and z, ξ is the scale length in y assumed for simulation of the 3-D effects of a finite length X-line,¹ and ν₀ is a constant. Corresponding to the initial current sheets, there are two peaks of ν_c at $(x,y,z)=(16di,64di,128di)$ and $(x,y,z)=(48di,64di,128di)$⁠, one in each current sheet center.

The electric field is updated by the electron momentum equation

where P_e is the thermal pressure of the electron fluid. The electron fluid is assumed to be isothermal, with T_e = const in this simulation. The velocity of the electron flow is obtained from Ampere's law

The magnetic field is updated using Faraday's law

The simulation is performed in 3-D with the domain size of $Lx×Ly×Lz=64di×128di×256di$ and the grid size of $Δx×Δy×Δz=0.25di×2.0di×2.0di$⁠, where d_i is the ion inertial length in the unperturbed current sheet ambient. Periodic boundary conditions are applied to the x and y directions, and a free boundary condition is applied to the z direction. In this study, the spatial length is normalized to d_i, time to $Ωi0−1$ with $Ωi0$ being the asymptotic ion gyro-frequency, the magnetic field to the asymptotic field strength B₀, and the ion number density to N₀. The velocity is normalized to the asymptotic Alfvén velocity $vA0$⁠, the electric field E to $E0=vA0B0$⁠, and the Poynting flux S to $S0=vA0B02/α$⁠.

The case shown in this paper corresponds to case 4 of Liang et al.,¹ which studied the KAWs generated in reconnection, with the initial guide field $By0=0.5, βi0=βe0=0.1$⁠, electron to ion temperature $Te0/Ti0=1$⁠, time step $Δt=0.05,ν0=1.0, ξ=5.0$⁠, and $δ=1.0$⁠. Initially, 100 particles are loaded in each grid cell in the asymptotic region to make sure the resulting ion velocity distributions as well as the ion heating are statistically meaningful. Note that β controls some of the dispersion properties of the waves (such as the dispersion relations for the fast and slow MHD modes). Moreover, the electron-to-ion temperature ratio may also play a role in the wave stability (by controlling the strength of parallel electric field and effectively the efficiency of Landau resonance). Generally, $Te/Ti$ is small in the magnetosphere, but not necessarily in other plasma environments. It is worth mentioning that the results of the two currents are found to be similar. Thus, only half of the simulation domain $x×y×z=[0,32]×[0,128]×[0,256]$ is discussed in this paper, for which the initial current sheet is at x = 16. Note that in this study, we only focus on the case with $β≫me/mi$⁠, i.e., in the KAW regime, as in magnetic reconnection in the magnetosphere. If electron kinetic effects are taken into consideration, the electron inertial effects as well as the regime of inertial Alfvén waves (IAWs) may also be a candidate. It should be noted that inertial effects can also become important in the limit of small $k∥$⁠.²⁰ Such an effect could potentially occur in a location with strong magnetic shear, but is typically associated with significant dissipation. Within the low frequency regime in reconnection, $k∥=0$ may be present in the tearing instability that occurs in the diffusion region. The waves seen in our simulation are found in the KAW regime dominated by ion gyroradius effects.

Our previous study has demonstrated that low frequency KAWs are a common feature in magnetic reconnection with various X-line lengths and guide field strengths.¹ In this paper, we examine the ion acceleration and heating by the KAWs generated during 3-D magnetic reconnection, in which a finite X-line length $2ξ=10di$ as well as a guide field $By0=0.5$ are applied. As indicated in our previous study, with a finite X-line length, i.e., $2ξ<30di$⁠, the resulting reconnection and wave structures are 3-D in nature, and the waves in the reconnection bulge are identified as KAWs propagate outward from the reconnection region along the magnetic field lines, with a propagation speed of 0.91, consistent with the phase speed (0.89) obtained from the linear dispersion relation of KAWs. The spectrum analysis shows that the perpendicular wave number $k⊥$ is dominant and much greater than the parallel wave number $k∥$⁠. It has also been found that the resultant KAWs are associated with considerable parallel electric field $E∥$ and parallel Poynting flux $S∥$⁠. The polarization relations are also consistent with those of KAWs.

An overall view of the 3-D magnetic field lines in reconnection as well as perturbation structures corresponding to the KAWs is shown in the left plot of Fig. 1, and a zoomed-in view around the iso-structure of $S∥$ indicated by the black rectangle is shown on the right. The X-line is generated at x = 16, z = 128 and is centered at y = 64. There are four types of colored field lines in Fig. 1 Lines 1–4, with Lines 1 (red) and 2 (green) being reconnected while Lines 3 (violet) and 4 (orange) not. The arrow on each field line represents the direction of the magnetic field. Lines 1 and 2 distinguish from each other by the dissipation region (where they reconnect, with B_z = 0). For Lines 1, the dissipation region is nearby or they are just through the X-line. Thus, the X-line can be traced through z = 128 from Lines 1. In contrast, Lines 2 reconnect through the quasi-steady reconnection layer above and below the X-line. Although Lines 3 are not yet reconnected, there are perturbations in the normal direction (x) around z = 128 due to the ion bulk inflow. Thus, separatrices are formed along the boundaries between Lines 1 and 3. Different from Lines 1–3, Lines 4 still remain unperturbed as the background magnetic field.

The iso-surfaces of the parallel Poynting flux with $S∥=−0.15$ are also plotted in Fig. 1, where $S∥=(δE×δB)·B¯/|B¯|, δE$ and $δB$ are the perturbed fields relative to the initial equilibrium field configurations, and $B¯$ is the spatially averaged field. It can be seen that the perturbations are spatially structured along the magnetic field lines, with $k⊥≫k∥$ [see also Fig. 3(a)]. As shown in our previous work,¹ these wave-like perturbations that are generated from the X-line and propagate along the magnetic field lines are KAWs, carrying parallel Poynting flux. The negative $S∥$ is seen in Fig. 1 in the region with x < 16 and z > 128, or x > 16 and z < 128, because there the KAWs propagate outward against the local magnetic field. Correspondingly, there are also two positive perturbations (not shown) at x < 16 and z < 128, or x > 16 and z > 128, where the KAWs propagate outward along the direction of the local magnetic field.

To reveal the ion dynamics in reconnection, the bottom three plots in Fig. 2 show the ion velocity distributions in the v_y–v_x, v_z–v_x, and v_z–v_y space inside the quasi-steady reconnection layer, where the magnetic field is dominated by Lines 2. The ions examined are inside a box with $15≤x≤17, 58≤y≤70$⁠, and $145≤z≤155$⁠, as indicated by the blue rectangle on the xz contours of B_y, which is shown in the top plot of Fig. 2. As shown in the previous study,¹ the structure of B_y around the reconnection bulge region corresponds to the Alfvénic parallel current, with an elongated pattern along the field lines at the wave front of the KAWs (not clearly seen in the xz plane in this case with an initial $By0=0.5$⁠). Overall, there is a bulk flow speed in −y and $+z$⁠, nearly symmetric about v_x = 0, while little bulk speed is seen in x. The ions are heated in both y and z directions, but predominantly in y. This is due to the local magnetic configuration. Around the current sheet center (x = 16), B_z is quite small while B_y is large. Thus, y is approximately parallel to the magnetic field. This kind of velocity distribution is consistent with previous studies of reconnection.^7,8,29,30 Reasonable mechanisms for the ions in the steady reconnection region to gain the bulk speed are through the acceleration by the reconnection electric field^7,8 and/or by slow shocks.²⁹ 

For the outward propagating KAWs generated during reconnection, the wave front is ahead of the leading bulge of field configuration (approximately z > 165 in Fig. 2), where the magnetic field lines are almost aligned in the yz plane since local $δBx$ is very small. Consequently, we examine the propagation of the KAWs primarily in the yz plane, at x = 13, $y≤64$⁠, and $z≥128$⁠. As shown in Fig. 3(a), five time instants, with t = 90, 120, 140, 160, and 190, are selected to track the wave patterns. It is worth mentioning that the perturbations of B_y elongated along the magnetic field lines (⁠$k⊥≫k∥$⁠) are KAWs. At each time, four regions, R₁–R₄, are selected to record the ion phase space properties and to examine the local acceleration and heating of ions by the KAWs, as indicated by the four blue rectangles along the path of the traveling waves. The four regions R₁–R₄ are centered at y = 48, y = 36, y = 24, and y = 12, respectively. Note that there is a shift in z with time of these regions following the wave pattern. The black arrows in Fig. 3(a) indicate the local magnetic field direction. From Fig. 3(a), top to bottom, the propagation of the waves can be summarized as follows. As reconnection takes place, the KAW perturbations start to be generated from the X-line and propagate outward (to z > 128 in Fig. 3), opposite to (since x < 16 for Fig. 3) the local magnetic field at an early time t < 90. At t = 90, the waves have reached region R₁ while the other three regions remain unperturbed. From t = 90 to 140, the waves continue propagating, passing regions R₂ and reaching R₃ in turn. Later at t = 190, the waves have almost filled up the four regions along the propagation path. Note that beside the propagation in the parallel direction, there is also a small displacement with time in the perpendicular direction due to the flow convection.

Figures 3(b) and 3(c) present the velocity distributions of ions in region R₂ at the five times in the $v⊥$–$v∥$ and $v⊥$–v_x space, where $v⊥$ and $v∥$ are the perpendicular and parallel particle velocities in the yz plane and v_x is the other perpendicular component. Note that $v⊥$ in the yz plane and v_x can be approximately regarded as two independent perpendicular components since $δBx$ is negligible in region R₂, as mentioned above. At the early time t = 90 when the waves have not reached region R₂ yet, the velocity distributions in both planes are Maxwellian, and the ions behave as unperturbed, cold background particles that are isotropic with a zero bulk velocity. At t = 120, some of the ions are accelerated as the wave perturbations have reached this region, forming a new ion beam with a finite parallel bulk velocity. At t = 140, the accelerated ions have begun to thermalize, in both the parallel and perpendicular directions. Note that there is also a small drift velocity $vd≈0.14$ in the perpendicular direction in both the core population and the accelerated beam due to the convection. Later at t = 160, more and more ions have been accelerated into the parallel beam, and the accelerated beam ions continue to be heated. Meanwhile, the core population and the accelerated beam start to interact with each other. Finally, at t = 190, the accelerated beam becomes more isotropic. Moreover, the perpendicular drift speed has also increased to $vd≈0.33$ as the reconnection bulges and thus a stronger convection electric field propagates through region R₂. The increase in v_d is relevant to the relative positions between region R₂ and the reconnection bulge. Seen from Fig. 3(b), the relative field-aligned velocity v_r between the accelerated beam and the core beam is found to be $vr≈0.85$–0.95, which is $∼1.13$–$1.27v¯A$⁠, consistent with previous studies by test particle simulations³¹ and observations of the solar wind.³² Besides, no bulk velocity is found in the x direction through the five time samples.

It is also seen from Fig. 3(b) that the parallel bulk velocity of the accelerated ion beam is about $u∥≈−0.9$⁠, which is nearly the same as the phase velocity of the wave front $Vpkaw=−0.91$⁠.¹ The averaged local Alfvén speed is $|v¯A|≈0.75$ and thus $|u∥|≈1.2|v¯A|$⁠. (Note that theoretically the wave phase speed of KAWs is larger than v_A when $βe>me/mi$⁠). Previous studies have shown that the ions can be accelerated up to the local Alfvén velocity in the outflow region of reconnection. However, here we have shown that the presence of KAWs can break this limit and accelerate ions to a super-Alfvénic velocity. Furthermore, the ions may be trapped in the wave potential as they are accelerated because $v∥$ agrees well with the wave phase speed. Previous studies have found that electrons can be trapped in KAWs and the distribution function becomes flattened.³³ In those simulations, the maximum trapping velocity is about 0.5, which is much closer to the core beam with a bulk $u∥=0$ in velocity space. In contrast, in our simulation, the trapping velocity associated with the accelerated ions is larger, and there are two ion populations rather than one flattened beam.

Our simulation shows that the ion heating by KAWs in reconnection is through different mechanisms in two different stages. In the earlier stage, when the ions are accelerated to form the beam, the ions are heated in the directions parallel and perpendicular to the magnetic field by Landau resonance and perpendicular stochasticity, respectively. In a later stage, ion beam-plasma instability^34,35 further contributes to the ion heating. In the following, we present results associated with these heating mechanisms.

In order to understand the acceleration of ions, shown as the accelerated beam in Fig. 3, the motion characteristics of two typical ions P1 and P2 belonging to the accelerated beam are examined. Figure 4(a) depicts the evolution of $v∥$ of ions P1 and P2 (blue and red solids, respectively) as well as the evolution of the mean parallel component $v¯∥$ (black dash-dot line) averaged over 800 accelerated particles from t = 70 to 160, while Figs. 4(b) and 4(c) show the particle energy E_ion in the laboratory frame of P1 and P2, respectively, as a function of the parallel velocity $−v∥′=v∥w−v∥$ in the wave frame within the same time interval, where $v∥w$ is the parallel velocity of the KAWs. At early times, e.g., t = 70, the ions have smaller $|v∥|$ and their kinetic energy is oscillating with smaller amplitudes. In the meantime, no visible accelerated beam is formed, since the ions have not been influenced and accelerated by the waves yet. Within $70<t<100$⁠, the ions are continually accelerated (⁠$v∥$ is more negative), gaining energy from the wave and resulting in the increases of both $|v∥|$ and E_ion. These two ions join the accelerated ion beam within $100<t<110$⁠. Later for t > 110, the parallel velocities of ions P1 and P2 remain nearly constant, with small oscillations around the propagation speed of the waves, indicating that the ions are trapped by the waves. In addition, the ion energy exhibits an oscillating feature around $−v∥′≈0$ after about t = 110, showing the characteristics of oscillations in the Landau resonance. The energy oscillation of particles is a combination of both the parallel motion and the perpendicular gyro-motion around the guiding center.

Before discussion of on heating, we first examine the trapped particle motion in the KAWs in t > 110. The left panel of Fig. 5 illustrates the contours of B_y (one component of the dominant magnetic field perturbation for KAWs) in the yz plane at x = 13, t = 120, 130, 140, and 150, in the region associated with the projected trajectories of the two typical ions P1 (blue) and P2 (red). The dots in each trajectory projection indicate the instantaneous positions of the ions at each corresponding time. At t = 120, it is seen that the wave front has already encountered the two ions and their dynamics will be influenced by the wave during the subsequent interval. By comparing the four time instants, it is also found that both ions move together with the propagation of certain wave phases, which provides evidence of particle trapping by the wave potential. Note that at t = 140, the two ions are just located within region R₂, contributing to the accelerated ion beam examined in Fig. 3. From t = 110 to 160, the average speed of the two ions is $∼0.9$⁠, which is consistent with the $u∥$ of the accelerated beam, and also consistent with the wave phase speed $|Vpkaw|=0.91$⁠.

Comparing the trajectories of the two particles in the left panel of Fig. 5, ion P1 seems to have a wavy movement in the anti-parallel field direction, while ion P2 seems to move more along a straight path. To understand this difference, the 3-D trajectory of ion P1 is shown in the right panel of Fig. 5, in which the small blue balls associated with the colored trajectory correspond to the blue dots in the contour plots, and the projections of trajectory in the x = 14, y = 60, and z = 150 planes are shown by gray lines. It is seen that the velocity of ion P1 also has a v_x component, which contributes to its gyro-motion. Consequently, the trajectory of ion P1 is a combination of the direct motion with the wave potential in the anti-parallel field direction and the gyro-motion in the perpendicular direction. Note that the motion in the anti-parallel field direction is dominant, which together with the gyro-motion results in a wavy projection in the yz plane. By comparison, the gyro-motion of ion P2 is much weaker, and thus its trajectory is nearly a straight line along the magnetic field.

The motion of the two ions discussed above are further illustrated in Fig. 6, in which the evolution of the velocity components of ion P1 (top row) and ion P2 (bottom row) are shown for t = 110 to 160, within the same time interval as Fig. 5. For each ion, the velocity evolution is examined in the $v⊥$–$v∥$⁠, v_x–$v∥$⁠, and $v⊥$–v_x space in the laboratory (i.e., simulation) frame and the $v⊥′$–$v∥′$ space in the wave frame, where $v⊥′=v⊥−v⊥w, v∥′=v∥−v∥w$⁠, and $v⊥w$ is the perpendicular velocity of the KAWs. First of all, it is seen that the parallel speeds of the two ions are $|v∥P1|≈|v∥P2|≈0.9$⁠, which is again consistent with the parallel speed of the waves. Second, the amplitudes of $v⊥$ and v_x of the same ion are roughly equal, indicating that v_x is consistent with the perpendicular gyro-motion as mentioned above. However, ion P2 has a much smaller perpendicular speed than ion P1, resulting in a much smaller gyroradius and thus a straighter trajectory projection in the yz plane. Moreover, there are small increases of the mean $v⊥$ in later times for both P1 and P2, as also seen in Fig. 3. Third, it is more explicit to see that there is a good consistency between the wave propagation and particle movements in the wave frame, as shown in the last column of Fig. 6, which, again, suggests that the ions are trapped by the waves.

For ion P1, there is a drift speed in the perpendicular direction (see Fig. 6), leading to the asymmetry of the gyro-motion and contributing to the oscillations of ion perpendicular energy in both the laboratory frame and the wave frame, with the frequency relative to its gyro-frequency. In the parallel direction, however, the period of energy oscillations is found to be about 30–$60Ωi0−1$ based on the guiding center motion of particles since the parallel scale length of the waves is much longer than the gyro-radius. For trapped particles, it is possible to estimate the nonlinear bounce frequency, ω_B, by considering the force exerted by the parallel electric field of the waves on the particle,

where m_i is the proton mass, e is the proton charge, and χ is a field-aligned coordinate taken about the minimum of the wave potential. Note that for the KAWs, the dominant wave vector is $k⊥$ while $k∥$ is very small. The parallel wave length (⁠$λ¯∥∼50$–60d_i within t = 100–200) is very large compared with the displacement of the trapped particles in the wave potential so that $k∥χ≪1$⁠. So it is reasonable to use the small angle approximation here for the trapped particles. The bounce frequency normalized to the proton cyclotron frequency can then be written as

From Faraday's law, $k∥δE⊥∼ω δB⊥$⁠, we have

By following the trapped particles in the simulation, the wave fields experienced by the trapped particles are estimated as $δE∥/δE⊥∼0.002/0.1=0.02, δB⊥/B0∼0.15$⁠, and the wave frequency $ω/Ωi0∼1/8$ − 1/3, which results in $ωB/Ωi0∼0.0194$–0.0316. Correspondingly, the bounce period is estimated to be about 32–$52Ωi0−1$⁠. In addition, the magnetic mirror force, $−μ∇∥B∼−i(Ti/B0)k∥δB∥$⁠, may also contribute to the trapping due to the KAWs, with a corresponding bounce frequency

where $ρi=(Ti/mi)1/2/Ωi$⁠. From the simulation, we have $ρi∼0.7, k∥∼0.13, δB∥/B0∼0.1$⁠, and $Ωi/Ωi0∼1$ for the trapped ions, and thus the corresponding bounce frequency is $ω/Ωi0∼0.0288$⁠, which corresponds to a bounce period of $∼35Ωi0−1$⁠. The estimation indicates that both the parallel electric field and the magnetic mirror force can lead to a bounce frequency consistent with that measured from the particle orbits, but the trapping due to the magnetic mirror force is more likely dominant because the mirror force is larger than the force of the parallel electric field. As the trapped particles move along the magnetic field, there is an effective potential energy, W_eff, associated with the waves, with

which is estimated to be $∼0.03$ based on the simulation data, again consistent with the parallel energy change of ion P1 in the wave frame (not shown).

Figure 7 presents the ion velocity distributions in the $v⊥$–$v∥$ space at t = 140 for all the four regions R₁ to R₄. At this moment, the leading wave front has reached region R₃ and on its way to R₄. Thus, the ions in region R₄ are still in an undisturbed Maxwellian distribution. In region R₃, only a small portion of the ions are accelerated, and the accelerated ion beam is still developing. This situation is quite similar to region R₂ at t = 120. In contrast, a large population of ions have been accelerated in region R₁, with a parallel bulk velocity of the accelerated ions $|u∥|≈0.7$⁠, which is smaller than that of the accelerated beams in the other regions at this time. Heating of the accelerated beams is also observed. Meanwhile, the accelerated beam and the core population tend to mix with each other in region R₁. This trend further develops with time, so that later at t = 190, the parallel velocity of the accelerated beam in region R₁ (to be illustrated in Fig. 12) further decreases to $|u∥|≈0.56$ and the two populations are almost mixed together.

From Fig. 7, it is also seen that the perpendicular drift velocity v_d is larger in region R₁ than in the other regions at the moment. Meanwhile, Fig. 3(b) shows that v_d of both beams is larger at a later time than an earlier time. In the former situation, region R₁ is closer to the reconnection bulge than the other regions. As for the latter, the locations of the plotted distributions are also closer to the bulge at the later time. Such a feature indicates that ions closer to the bulge have larger v_d. Consequently, it is worthwhile to discuss the effect of the relative distance of a certain velocity distribution relative to the reconnection bulge. Here we take region R₁ at t = 190 as an instance. Figure 8 presents the contour plots of $S∥, T∥$⁠, and $T⊥$ in both the xz (together with the projected magnetic field lines) and the yz planes, where $T∥$ and $T⊥$ are the ion parallel and perpendicular temperatures, respectively. The horizontal red line in Fig. 8 indicates the position of the trailing edge of the reconnection bulge, and the black rectangle shows region R₁. At t = 140 before the time shown in Fig. 8, the edge of the reconnection bulge was at $z≈164$⁠, and the range of region R₁ was $42<y<54,156<z<184$⁠. Thus, a small population of the ions examined in this region is also located in the quasi-steady reconnection region below the bulge. At t = 190, the edge of the bulge is at $z≈182$ and region R₁ is $42<y<54, 168<z<196$⁠. At this moment, the bulge is almost at the center of region R₁, and a significant portion of the ions examined are in the quasi-steady reconnection region. The change in the relative positions between the ions and the reconnection bulge indicates that region R₁ gradually merges into the quasi-steady reconnection region with time. As a result, the velocity distribution in region R₁ at t = 190 (see Fig. 12) looks more similar to that in the quasi-steady reconnection region (see Fig. 2), both with a larger perpendicular bulk velocity. Moreover, due to the generation of the accelerated ion beam in the parallel direction, there is a strong increase in the parallel temperature, with $T∥>T⊥$⁠, inside the wave pattern. The parallel velocity separation between the accelerated beam and the core population causes the enhancement in $T∥$⁠, which may not mean an ion heating, but the ion distribution at t = 190 (see Fig. 12) shows that the two beams are evolving into a single, heated population, with a significant parallel heating. Wang and Lin³⁵ have carried out a 2-D hybrid simulation on ion beam-plasma interactions and also found that the bulk speed of the beam and core ions merged with each other due to the instability.

Previous studies have shown that parallel heating by KAWs may result from the direct heating by the parallel electric field³⁶ or from the nonlinear wave-particle interactions such as ion Landau damping when $β∼1$⁠.³¹ Cheng et al.⁸ have found that the parallel electric field can accelerate the ions with meandering orbit near the reconnection current layer (near the electron diffusion region) to gain a considerable $v∥$⁠. In our study, KAWs are found to be present throughout the leading bulge region of reconnection, where $E∥∼0.005$⁠. Our simulation shows that in this region, nonlinear wave-particle interactions are the main cause of the acceleration and heating of ions. As discussed above, ions can be trapped by the waves after a short period of energy buildup. An evidence of Landau resonance is also found. It is thus important to compare the heating rate of the ions with the damping rate of the KAWs. We calculate $T∥$ and $T⊥$ separately for the accelerated beam and the core population to avoid the contribution of the separation of their mean parallel velocity to the total dynamic temperature T. For each population, the ion temperature in each region at each time can be estimated by

in which $v¯∥$ (⁠$v¯⊥$⁠) is the mean parallel (perpendicular) velocity in the yz plane for each ion beam. The estimation on the heating rate of the accelerated beams is shown in Fig. 9 for the four regions R₁ to R₄. The asterisks, open diamonds, and open squares represent the total ion temperature T, $T⊥$⁠, and $T∥$⁠, respectively, as a function of time. The red (blue) dashed line in each plot is used to estimate the perpendicular (parallel) heating rate. Note that for different regions, the estimation is performed within different time intervals because of the time differences in the generation of the ion beams. In addition, the heating mechanism becomes complicated due to the beam-plasma interaction (i.e., interaction between the accelerated ion beam and the core population) after a duration of approximately 50–60 $Ωi0−1$ as well as the influence of reconnection as mentioned above. Consequently, we estimate the heating rate within 50–60 $Ωi0−1$ upon the generation of the accelerated ion beam. With the consideration of regions R₁–R₄, the perpendicular heating rate $γ⊥$ of the ions can be estimated as $γ⊥/Ωi0≈0.0047$–0.0066, and the parallel heating rate $γ∥$ is roughly estimated as $γ∥/Ωi0≈0.0068$–0.0093. Chaston, Bonnell, and Salem¹⁹ have estimated both the parallel and perpendicular heating rate for the ions scattered in the KAWs within fast flows moving earthward. They have shown that parallel heating rate due to Landau damping is about $10−6$–$10−1.8$⁠. Our estimation based on the hybrid simulation is consistent with their results. For the perpendicular direction, they have indicated that the heating is due to the stochastic process and the heating rate is about $10−1.34$–$101.26$⁠, which is much larger than our estimation. Note that their estimation is based on a larger $k⊥ρi=0.6$–1.6, compared with $k⊥ρi≈0.5$ in our simulation. In our study, it is found that $γ⊥$ and $γ∥$ are about the same order, with $γ∥≥γ⊥$ although $T⊥>T∥$⁠. They both contribute to a heating rate of the total temperature T of $γ/Ωi0≈0.005$–0.007. In our previous work¹ the damping rate of the KAWs was found to be $∼−0.0035$–$−0.0062Ωi0$⁠. Consequently, it is evaluated that the ion heating rate is roughly consistent with the wave damping rate. The matching between the heating and damping rate also provides evidence of ion heating by the KAWs. In addition, it should also be noted that there is an obvious heating of the accelerated beam while almost no heating is found for the core population. Consequently, the heating of the whole ion population is dominated by the accelerated ion beam. It is also worth mentioning that the ion beam-plasma interaction may also cause modifications of the KAWs and the associated ion heating.

Previous studies on ion heating by KAWs have also shown that ions can be heated in the perpendicular direction.^{13–15,17,18} Johnson and Cheng¹⁴ have found on the basis of stochastic heating theory that ions are perpendicularly heated when the amplitude of $E∥$ is small. They have also shown the possibility of the perpendicular sub-harmonic stochastic heating by the waves. In this situation, the magnetic moment μ is not conserved and particle orbits can become chaotic in a few wave periods. Such a mechanism has also been reported in obliquely propagating Alfvén waves with sufficiently large amplitude or with lower amplitude if the wave spectrum is of multiple modes.^13,15 Using test particle simulations, Voitenko and Goossens¹⁷ have proposed a mechanism for ion heating by low frequency KAWs through cross-field energization, with the criteria of the cross-field heating expressed as $k⊥ρivA/vti[(1+k⊥2ρi2)/1+(1+Te/Ti)k⊥2ρi2−u/vA]Bk0/B0>1$⁠, in which v_ti is the ion thermal speed, and $Bk0$ and B₀ are the peak wave amplitude and the background magnetic field, respectively, and u is the bulk velocity. In our simulation, $vA/vti≈1.9, k⊥ρi≈0.5, Te/Ti≈0.5$⁠, and $Bk0≈0.3$⁠. Although the bulk flow of the main population is still much smaller than v_A, with $u/vA≪1$⁠, the threshold of the cross-field ion heating is not satisfied. Chandran et al.¹⁸ have also examined the ion heating by low frequency turbulence AWs or KAWs, with test particle simulations. It is found that the conservation of the magnetic moment μ is violated when the amplitudes of the fluctuating velocity or magnetic field are sufficiently large on the ion gyroradius scale (⁠$λ⊥∼ρi$⁠), leading to chaotic ion orbits and perpendicular heating. The criteria for the violation of μ conservation is $δv/v⊥≥0.19$⁠, in which δv is the amplitude of the velocity perturbation (due to the waves) and $v⊥=T⊥$⁠. In our simulation, it is estimated that $δv/v⊥≈0.66$–0.98 on the scale of $k⊥ρi<1$⁠, consistent with breaking of the first adiabatic invariant.

To examine the perpendicular stochastic heating in our simulation, we present in Fig. 10 the 3-D trajectories [Fig. 10(a) for t = 70–160], time evolution of $μ=Wion⊥/B$ [Fig. 10(b) for t = 76–156], with $Wion⊥$ being the ion perpendicular thermal energy, and time evolution of particle velocities [Fig. 10(c) for t = 110–160] in the $v⊥$–v_x space of three particles P3, P4, and P5. The black dots in Fig. 10(a) indicate the instantaneous positions of the three particles at t = 110. The black dashed line in Fig. 10(b) represents the average μ of about 800 particles. In Fig. 10(c), the velocity evolutions are plotted within the same time intervals and in the same format as Fig. 6. In the early stage (t = 70–100) as the waves have reached these particles, since the ions are being accelerated from a smaller $v∥$ [see Figs. 4(a) and 10(a) starting from $x≈16$ and $y≈80$], $ω−k∥v∥$ is large and the cyclotron or sub-cyclotron resonance conditions may be satisfied for a significant fraction of particles. During this period of time, a significant separation between particle trajectories is seen for the three ions, while their trajectories are very close at the beginning. Meanwhile, there is also an increase in μ generally, at $t≈88$⁠, showing evidence of stochastic perpendicular heating. In the later stage at t > 110, these ions are trapped in the waves, with $v∥≈v∥w$ and thus $ω−k∥v∥≈0$⁠, and the cyclotron resonance is no longer responsible for heating. Instead, the ion Landau resonance dominates during this period of time as discussed above [see Figs. 4(b) and 4(c)], leading to parallel heating of the ions. Over the time, locally (e.g., in R₂) there are constantly more ions in the core velocity distribution being accelerated and trapped by the waves when the KAWs are continuously generated from the reconnection site.¹ As a result, the accelerated beam population is persistently present locally behind the leading front of KAWs [see Fig. 3(b)], and the fraction of the accelerated ions in the velocity distribution increases with time. Hence, ion heating and acceleration take place at the same time. As also shown in Fig. 9, the parallel and perpendicular temperatures increase within the same time intervals.

After the ions are trapped (during $110<t<160$ for these three ions), their guiding centers are constrained by the field values near the potential minimum. In this situation, the total magnetic field at a guiding center roughly remains invariant. As shown in Figs. 5, 6, and 10, there is little variation in the gyroradius as well as the total perpendicular velocities of the ions after they are fully trapped. Correspondingly, no significant variations in μ are found for these ions.

In addition to the stochastic heating associated with the increase of μ and separation of ion orbits that are initially very close, more evidence that links the perpendicular heating to the field perturbations of KAWs is shown in Fig. 11, in which variations of B_x, B_y, B_z, $B⊥$⁠, E_x, E_y, E_z, and $T⊥$ in region R₁ as a function of time from t = 70–220 are presented. It is found that there is an increase in $T⊥$ from t = 80–140 (between the two vertical dashed lines in Fig. 11), as also seen in Fig. 9 for region R₁. Such increase is due to the perpendicular heating of the ions as the KAWs propagate through. As shown in the corresponding electromagnetic perturbations, there is an increase in B_y, B_z, and E_x, among which the increase in B_y and B_z will contribute to an increase in $B⊥$ in the yz plane. Note that for the KAWs here, k_x is the perpendicular wave number. As a result, E_x is the perpendicular electric perturbation and $B⊥$ is the perpendicular magnetic perturbation of the KAWs. The coherence in E_x, $B⊥$⁠, and $T⊥$ provides the evidence of perpendicular heating by KAWs. In contrast, there is no coherence among other components, such as B_x, E_y, and E_z.

The temperature $T⊥$⁠, however, undergoes a drop at t = 140 after being heated, followed by another increase at $t≥160$⁠. Such a drop can also be found for all the four regions as shown in Fig. 9. This behavior may be due to the interaction between the accelerated ion beams and the core plasma, which can be seen from the evolution of the velocity distribution shown in Fig. 12 for region R₁ at t = 90, 110, 130, 150, 170, and 190. As discussed above, the effect of beam-plasma interaction becomes more obvious with the increase of the accelerated beam fraction. The beam-plasma interaction first causes an excitation of instability at the expense of the particle energy subsequently followed by plasma heating due to further wave-particle interactions. At t = 170, the two beams nearly merge into a single population and $T⊥$ increases again. Then, $T⊥$ oscillates around a constant value as the velocity distributions merge. As a result, the overall parallel temperature $T∥>T⊥$ due to reconnection, with $T∥$ enhanced approximately by a factor of $∼2.36$ from the initial value, and $T⊥$ enhanced by a factor of $∼1.31$⁠. Wang and Lin³⁵ have indicated that Alfvénic waves may generate if the relative stream velocity of the beams is larger than $4vA$ for a concentration of $nb/n0∼10%$⁠, where n_b is the number density of the beam particles and n₀ is the number density of the background particles. In our simulation, however, there is little chance for such Alfvénic waves to generate because the relative velocity is much smaller (about 1.13–$1.27vA$ as mentioned above) although $nb/n0$ may be larger (see the velocity distribution). Nevertheless, Winske and Omidi³⁷ have shown that ion-ion cyclotron instability may be excited for a concentration of $nb/n0∼50%$ at a much smaller relative velocity $∼vA$⁠. The beam-plasma interaction in our simulation is consistent with this regime.

In this study, ion acceleration and heating are investigated on account of the KAWs generated during magnetic reconnection with a 3-D hybrid model. It is found that the ions can be accelerated by the KAWs and heated by wave-particle interactions through comprehensive nonlinear physics. The heating is due to different mechanisms at different stages of the wave-particle interactions, including Landau resonance, perpendicular stochastic orbits, as well as the ion beam instability. The main results are summarized as follows.

1. The ion acceleration and evolution of the accelerated beam can be depicted as follows. In the beginning, as the KAWs propagate away from the reconnection region and occupy the leading bulge region of reconnection, they pick up and accelerate ions, leading to the development of an accelerated beam distribution in addition to the core Maxwellian population. The ion beam is fully accelerated up to a super-Alfvénic velocity, with the average parallel speed $v¯∥$ consistent with the propagation speed of the KAWs, and the accelerated population is trapped in the wave field. Due to the convection electric field of reconnection, both the core population and the accelerated beam have a small perpendicular drift velocity. The perpendicular drift speed increases as the reconnection bulge approaches with the reconnection layer.

2. Ion heating is found in both the parallel and the perpendicular directions. During the stage of ion acceleration, evidence of perpendicular stochastic heating is found, with an increase in the gyro-averaged magnetic moment μ and a separation of ion trajectories that are initialized very close. As the ions are trapped by the waves, parallel heating due to Landau resonance is present. Consequently, both parallel and perpendicular heating is seen upon the generation of the accelerated beam. The perpendicular heating is closely correlated with the strength of the perpendicular electric and magnetic fields of the KAWs. Although both the accelerated and core beams are heated, the dominant heating arises from the accelerated beam population.

3. As time goes on, more ions in the core velocity distribution are being accelerated and trapped by the waves when the KAWs are continuously propagating from the reconnection site. Effective interaction between the accelerated ion beam and the core plasma takes place in later times (in about 50–60 $Ωi0−1$ after the generation of the ion beam) with the increase of the fraction of the accelerated beam. As a result, the ion beam and the core population merge together in the velocity space. The ion beam-plasma instability also contributes to the enhancement of the ion temperature. As a result, the overall parallel temperature $T∥>T⊥$ due to reconnection.

4. The heating rate of the accelerated beam is consistent with the damping rate of the KAWs, showing that the nonlinear wave-particle interactions contribute to the ion heating.

Our simulation demonstrates the ion dynamics associated with the KAWs generated in magnetic reconnection. The acceleration and heating of electrons as well as their effects on the KAWs are not included in the hybrid simulation. More complete physics of the KAWs must also include the electron kinetic physics, e.g., the electron heating and interaction with the KAWs.

This work was partly supported by the National Natural Science Foundation of China (No. 11675038) and the Funds of Dalian Young Talents with Grant No. 2015R001 and the DoE Grant No. DE-SC0010486 and NSF Grant No. AGS 1405225 to Auburn University. Computer resources were provided by the Alabama Supercomputer Center and the cluster for plasma major of Peking University.