PERPENDICULAR ION HEATING BY LOW-FREQUENCY ALFVÉN-WAVE TURBULENCE IN THE SOLAR WIND

Beginning in the 1960s, a number of authors developed steady-state hydrodynamic models of the solar wind in which the temperature was fixed at the coronal base and the solar wind was heated by thermal conduction (e.g., Parker 1965; Hartle & Sturrock 1968; Durney 1972; Holzer & Leer 1980; Leer & Holzer 1980). For realistic values of the coronal temperature and density, these models were unable to reproduce the large flow velocities of fast-solar-wind streams and large proton temperatures observed at 1 AU, suggesting that the solar wind is heated above the coronal base by some additional mechanism. Further evidence for extended, non-conductive heating has been provided by measurements from the Ultraviolet Coronagraph Spectrometer (UVCS), which show radially increasing minor-ion temperatures in coronal holes (the open-magnetic-field-line regions from which the fast wind emanates) at heliocentric distances r between 1.5 R_☉ and 3.5 R_☉ (Kohl et al. 1998; Li et al. 1998; Esser et al. 1999; Antonucci et al. 2000). Identifying the physical mechanisms responsible for this heating and determining the heating rates of the different particle species are among the major challenges in the study of the solar wind at the present time.

One of the first mechanisms proposed to account for solar-wind heating was Alfvén-wave (AW) turbulence (Coleman 1968). The importance of AW turbulent heating is suggested by in situ measurements of ubiquitous, large-amplitude fluctuations in the velocity, magnetic field, and electric field in the interplanetary medium with power spectra that vary approximately as power laws over a broad range of length scales (Belcher & Davis 1971; Goldstein et al. 1995; Tu & Marsch 1995; Bruno & Carbone 2005; Bale et al. 2005). The amplitudes of the measured fluctuations are positively correlated with solar-wind temperature (Grappin et al. 1990) and imply a turbulent heating rate (based on phenomenological turbulence theories) that is comparable to the observationally inferred solar-wind heating rate (Smith et al. 2001; Breech et al. 2009; Cranmer et al. 2009b; Stawarz et al. 2009).

The in situ measurements on which the above studies are based are limited to the locations where spacecraft have flown—that is, to r ≳ 0.3 AU. However, a number of considerations suggest that AW turbulent heating is important closer to the Sun as well. For example, the velocity and magnetic-field fluctuations measured in situ are often correlated in the sense of AWs propagating away from the Sun in the solar-wind frame (Belcher & Davis 1971; Tu & Marsch 1995; Bavassano et al. 2000), indicating that these waves originate at or near the Sun. Throughout the chromosphere, filamentary structures (spicules) undergo transverse motions consistent with outward-propagating AWs carrying an energy flux sufficient to power the solar wind (De Pontieu et al. 2007). The corona is also pervaded by transverse velocity patterns (seen in measurements of the Doppler shift of an Fe⁺¹² emission line), which propagate along the local magnetic-field direction at speeds comparable to the local Alfvén speed $v_{A} = B/\sqrt{4 \pi \mu _{\rm p} n m_{\rm p}}$, where B is the magnetic-field strength, n is the proton number density, and μ_p is the mean molecular weight per proton (Tomczyk et al. 2007).

Radio observations are also consistent with AW turbulent heating of the extended corona and inner solar wind. For example, Faraday rotation measurements have been used to probe magnetic fluctuations at 2 R_☉ ≲ r ≲ 15 R_☉ (Hollweg 1982; Sakurai & Spangler 1994; Mancuso & Spangler 2000). If the measured magnetic fluctuations are AWs, then they carry an energy flux that is sufficient to drive the solar wind (Hollweg 1982; Andreev et al. 1997; Hollweg et al. 2010). Coles & Harmon (1989) analyzed spectral-broadening and phase-scintillation data to measure the power spectrum of electron density fluctuations at radii as small as 5 R_☉. The density fluctuations have significant amplitudes over a broad range of spatial scales, consistent with passive-scalar mixing of background density fluctuations or entropy modes by AW turbulence (Harmon & Coles 2005). Because AWs become compressive at small scales (when the wavelength is ≲ the proton inertial length; Hollweg 1999), the rms amplitude of density fluctuations at small scales can be used to place an upper limit on the turbulent heating rate, and this upper limit exceeds the heating rates needed to produce the solar wind in theoretical models (Harmon & Coles 2005; Chandran et al. 2009).

Further support for the importance of AW turbulence for solar-wind heating is provided by several theoretical models that describe the radial evolution of AW turbulence in the solar wind, taking into account the cascade of wave energy from large scales to small scales, the loss of wave energy due to dissipation at small scales, and the wave-amplitude modulation and wave reflection that are induced by large-scale variations in the solar-wind velocity, density, and magnetic-field strength. These models incorporate observational constraints on the wave amplitudes at r < 3 R_☉ and r>0.3 AU, and lead either to temperatures, densities, and wind speeds that are consistent with observations (Cranmer et al. 2007; Verdini et al. 2010), or to heating rates that are sufficient to reproduce the observed temperature profiles (Cranmer & van Ballegooijen 2005; Verdini & Velli 2007; Chandran & Hollweg 2009). Wave reflection plays an important role in these models because the Sun launches only outward-propagating waves, and AWs propagating in the same direction in the plasma rest frame do not interact with one another (Iroshnikov 1963; Kraichnan 1965).

Taken together, the observational evidence and theoretical results described above suggest that AW turbulence plays an important, and perhaps dominant, role in the heating of the solar wind and coronal holes. However, there is another set of observations that has proven difficult to reconcile with heating by AW turbulence. In coronal holes at 1.6 R_☉ < r ≲ 3 R_☉, proton temperatures exceed electron temperatures, and minor-ion temperatures greatly exceed electron temperatures, indicating that ions are efficiently heated in the low-β conditions found near the Sun (Cranmer 2009a), where β = 8πp/B² is the ratio of the plasma pressure to the magnetic pressure. In addition, UVCS observations show that minor ions in coronal holes are heated in such a way that thermal motions perpendicular to the background magnetic field B₀ are much more rapid than thermal motions along B₀ (i.e., T_⊥ ≫ T_∥) (Kohl et al. 1998; Li et al. 1998; Antonucci et al. 2000). A similar temperature anisotropy is measured in situ at r>0.3 AU for protons in fast-solar-wind streams with β ≪ 1, despite the fact that (double) adiabatic expansion acts to decrease T_⊥/T_∥ (Marsch et al. 1982, 2004; Hellinger et al. 2006). Thus, in low-β conditions, ions are heated efficiently, and ion heating is mostly "perpendicular to the magnetic field."

These observations present a challenge for solar-wind-heating models based on AW turbulence for the following reasons. At the comparatively large scales at which most of the AW energy is concentrated (i.e., at the "outer scale" of the turbulence), the wave frequencies are much less than the proton cyclotron frequency Ω_p. For example, most of the AW energy detected remotely in the solar atmosphere is at periods ranging from several minutes to a few hours (Chashei et al. 2000; De Pontieu et al. 2007; Tomczyk et al. 2007). At r>0.3 AU, the bulk of the AW energy is at periods of hours and longer (Tu & Marsch 1995; Bruno & Carbone 2005). By comparison, Ω⁻¹_p ∼ 10⁻⁴ s in coronal holes at r = 2 R_☉ and Ω⁻¹_p ∼ 1 s at r = 1 AU. Moreover, in AW turbulence, the energy cascade is anisotropic, transporting energy primarily to small scales measured perpendicular to B₀ rather than small scales along B₀ (Shebalin et al. 1983; Goldreich & Sridhar 1995; Ng & Bhattacharjee 1996; Galtier et al. 2000; Cho et al. 2002). In wavenumber space, energy cascades primarily to larger k_⊥, and only weakly to larger k_∥, where k_⊥ and k_∥ are the components of the wavevector k perpendicular and parallel to B₀. Because the frequency of an AW is given by ω = k_∥v_A, this anisotropic energy cascade is inefficient at transporting wave energy to higher frequencies. At perpendicular scales (λ_⊥) of order the proton gyroradius ρ_p, the AW cascade transitions to a kinetic-Alfvén-wave (KAW) cascade (Bale et al. 2005; Howes et al. 2008b; Schekochihin et al. 2009; Sahraoui et al. 2009), and at λ_⊥ ≲ ρ_p the fluctuations dissipate. The rms amplitude of the magnetic-field fluctuations at λ_⊥ ≲ ρ_p is ≪B, and thus it seems plausible to assume that the KAW fluctuations damp at the same rate as linear KAWs at the same k. For the β values found in coronal holes and the solar wind, linear KAWs undergo significant electron Landau damping (Quataert 1998; Leamon et al. 1999). However, when β ≪ 1, ion thermal speeds are ≪v_A. Ions are thus unable to satisfy the resonance condition ω − k_∥v_∥ = 0 for Landau damping or transit-time damping, where v_∥ is the velocity component parallel to B₀, because ω/k_∥ ⩾ v_A for KAWs (Quataert 1998; Hollweg 1999). In addition, the KAWs produced by the anisotropic cascade do not undergo ion cyclotron damping, because ω ≪ Ω_p for these waves (and also |ω − k_∥v_∥| ≪ Ω_p) (Cranmer & van Ballegooijen 2003; Howes et al. 2008a). As a consequence, linear damping of KAWs by ions is negligible when β ≪ 1. Thus, if low-frequency AW/KAW turbulence damps according to linear Vlasov theory, then it is unable to explain the strong perpendicular ion heating that is inferred from observations of coronal holes and low-β fast-wind streams.

A number of studies have gone beyond linear Vlasov theory to investigate the possibility of perpendicular ion heating by low-frequency AW/KAW turbulence. Johnson & Cheng (2001), Chen et al. (2001), White et al. (2002), Voitenko & Goossens (2004), and Bourouaine et al. (2008) investigated the dissipation of KAWs and AWs with ω < Ω_p, finding that such waves cause perpendicular ion heating if the wave amplitude exceeds a minimum threshold. Dmitruk et al. (2004) and Lehe et al. (2009) simulated test particles propagating in the electric and magnetic fields resulting from direct numerical simulations of MHD turbulence at 0.1 ≲ β ≲ 10. They both found perpendicular ion heating under some conditions, but Lehe et al. (2009) argued that the perpendicular heating seen in both studies is due to cyclotron resonance and does not apply to the solar wind because it is an artifact of limited numerical resolution. Parashar et al. (2009) found perpendicular ion heating in two-dimensional hybrid simulations of a turbulent plasma, in which ions are treated as particles and electrons are treated as a fluid. Servidio et al. (2010) analyzed magnetic reconnection in simulations of MHD turbulence, and Drake et al. (2009) found that reconnection leads to perpendicular ion heating.⁴ In addition, Markovskii & Hollweg (2002) and Markovskii et al. (2006) investigated high-frequency secondary instabilities that are generated by KAWs near the gyroradius scale, and argued that such instabilities may be able to explain the observed perpendicular ion heating.

In this paper, we continue this general line of inquiry and address an important open problem: determining the perpendicular ion heating rate in anisotropic, low-frequency (ω < Ω_p), AW/KAW turbulence as a function of the amplitude of the turbulent fluctuations at the gyroradius scale. In Section 2, we develop a phenomenological theory of stochastic ion heating, obtaining an approximate analytic expression for the heating rates of different ion species in plasmas with β ≲ 1. We also present simulations of test particles propagating in a spectrum of AWs and KAWs to test our phenomenological theory and to determine the two dimensionless constants that appear in our expression for the heating rate. In Section 3, we discuss the implications of our results for perpendicular ion heating in coronal holes and the fast solar wind.

We consider ion heating by AW/KAW fluctuations with transverse length scales λ_⊥ (measured perpendicular to B₀) of order the ion gyroradius ρ = v_⊥/Ω (i.e., k_⊥ρ ∼ 1), where Ω = qB₀/mc is the ion cyclotron frequency, and q and m are the ion charge and mass. We assume that ρ ≳ ρ_p, where

$$\begin{equation} \rho _{\rm p } = \frac{v_{\perp {\rm p}}}{\Omega _{\rm p}} \end{equation} \tag{ 1 }$$

is the rms proton gyroradius in the background magnetic field,

$$\begin{equation} v_{\perp {\rm p}} = \sqrt{\frac{2k_{B} T_{\rm p}}{m_{\rm p}}} \end{equation} \tag{ 2 }$$

is the rms perpendicular velocity of protons, T_p is the (perpendicular) proton temperature, and m_p is the proton mass. If ρ ≫ ρ_p, then the gyro-scale fluctuations are AWs. If ρ ∼ ρ_p, then the gyroscale fluctuations are KAWs. For simplicity, we neglect the effects of fluctuations with either λ_⊥ ≫ ρ or λ_⊥ ≪ ρ. We also assume that

$$\begin{equation} \beta \lesssim 1. \end{equation} \tag{ 3 }$$

We define δv_ρ and δB_ρ to be the rms amplitudes of the fluctuating velocity and magnetic-field vectors at k_⊥ρ ∼ 1. Similarly, δE_ρ and δΦ_ρ are the rms amplitudes of the fluctuating electric field and electrostatic potential at k_⊥ρ ∼ 1. We assume that δv_ρ, δB_ρ, δE_ρ, and δΦ_ρ are related to one another in the same way that the magnitudes of the fluctuating velocity, magnetic field, electric field, and electrostatic potential are related in a linear KAW. Thus, since k_⊥ρ_p ≲ 1,

$$\begin{equation} \delta E_\rho \simeq \frac{\delta v_\rho B_0}{c}, \end{equation} \tag{ 4 }$$

δΦ_ρ ∼ ρδE_ρ, and

$$\begin{equation} q \,\delta \Phi _\rho \sim m v_\perp \,\delta v_\rho. \end{equation} \tag{ 5 }$$

The fractional change in an ion's perpendicular kinetic energy mv²_⊥/2 induced by gyroscale fluctuations during a single gyroperiod is then approximately

$$\begin{equation} \frac{2 q\, \delta \Phi _\rho }{m v_\perp ^2} \sim 2 \varepsilon, \end{equation} \tag{ 6 }$$

where

$$\begin{equation} \varepsilon = \frac{\delta v_\rho }{v_\perp }. \end{equation} \tag{ 7 }$$

When ε ≪ 1, an ion's kinetic energy is nearly constant during a single gyroperiod. The inequality ε ≪ 1 also implies that δB_ρ ≪ B₀, because of Equation (3). Thus, if ε ≪ 1, an ion's orbit in the plane perpendicular to B₀ closely approximates a closed circle. In this case, the ion possesses an adiabatic invariant of the form J = ∮pdq˜ that is conserved to a high degree of accuracy, where q˜ is the angular coordinate corresponding to the particle's nearly periodic cyclotron gyration and p is the canonically conjugate momentum (Kruskal 1962). In the limit of small ε, J is approximately equal to the magnetic moment μ = mv²_⊥/2B. The near conservation of J implies that perpendicular ion heating is extremely weak. In Appendix A, we present a calculation for electrostatic waves with k_⊥ρ ∼ 1 and ε ≪ 1 that illustrates how the leading-order terms in the time derivative of v²_⊥ are unable to cause secular growth in T_⊥.

On the other hand, as ε increases from 0 to 1, the fractional change in an ion's perpendicular kinetic energy during a single gyroperiod grows to a value of order unity. We treat the spatial variations in the electrostatic potential Φ at k_⊥ρ ∼ 1 as random or disordered, as is the case in turbulence or a spectrum of many randomly phased waves. Thus, when ε exceeds some threshold (whose value we investigate below), an ion's orbit in the plane perpendicular to B₀ becomes chaotic. In this case, the ion's orbit does not satisfy the criteria for the approximate conservation of J (Kruskal 1962), and perpendicular ion heating becomes possible (Johnson & Cheng 2001; Chen et al. 2001; White et al. 2002).

To estimate the rate at which ions are heated, we begin by considering the Hamiltonian of a particle of charge q and mass m:

$$\begin{equation} H = q\Phi +\frac{1}{2m} \left({\bf p} - \frac{q}{c} {\bf A} \right)^2, \end{equation} \tag{ 8 }$$

where A is the vector potential, and p is the canonical momentum. Hamilton's equations imply that

$$\begin{equation} \frac{d H}{d t} = q\frac{\partial \Phi }{\partial t} - \frac{q {\bf v}}{c}\cdot \frac{\partial {\bf A}}{\partial t}, \end{equation} \tag{ 9 }$$

where v = m⁻¹(p − qA/c) is the particle's velocity. The electric field is given by E = −∇Φ − c⁻¹∂A/∂t. The second term in Equation (9) is qv · E_s, where E_s = −c⁻¹∂A/∂t is the part of the electric field that has a nonzero curl. In AWs and KAWs with ω < Ω_p and k_⊥ρ_p ≲ 1, E_s is negligible compared to the total electric field in low-β plasmas (see Equation (46) of Hollweg 1999), which are our primary focus, and so from here on we neglect the second term in Equation (9).

When ∂Φ/∂t>0, a particle can gain potential energy, kinetic energy, or both. For example, if an ion interacts with an electrostatic wave with wavelength ≫ρ and frequency ≪ Ω, then the ion's guiding-center drifts with velocity cE × B₀/B²₀. The particle's kinetic energy undergoes small-amplitude oscillations due to its gyromotion. However, because its guiding center moves perpendicular to ∇Φ, there is no significant secular change in its kinetic energy. The ion's magnetic moment μ is almost exactly conserved, and the change in its total energy is almost exactly equal to the change in its potential energy.

On the other hand, if a particle enters a region in which ∂Φ/∂t>0 and then leaves this region, moving up and down the potential gradient, then it can gain kinetic energy as illustrated in Figure 1. The "wire-mesh" surface in the upper panel of this figure represents Φ(x₁, x₂, x₃, t) at x₃ = 0 at some initial time, and the lower panel shows Φ(x₁, x₂, x₃, t) at x₃ = 0 at a later time. We take the maximum of Φ(x₁, x₂, 0, t) to be located at σ = 0, where $\sigma \equiv \sqrt{\vphantom{A^A}\smash{\hbox{${x_1^2 + x_2^2}$}}}$. We have assumed that ∂Φ/∂t>0 at σ ≲ σ₀ and ∂Φ/∂t = 0 at σ ≳ σ₀, where σ₀ is the approximate radius in the x₁–x₂ plane of the "potential-energy hills" that appear in the figure. The thick solid line shows the value of Φ along the trajectory of a particle moving in a straight line in the x₁–x₂ plane. Because ∂Φ/∂t>0, the potential-energy hill is shorter when the particle is "climbing up" and higher when the particle is "rolling down." The particle thus experiences a net gain of kinetic energy from "rolling over the hill."

Zoom In Zoom Out Reset image size

We now estimate the rate at which ions are heated by AW or KAW fluctuations with λ_⊥ ∼ ρ. The condition λ_⊥ ∼ ρ is intended to encompass structures with λ_⊥ = 0.5ρ, which we invoke below when discussing Equation (24). However, we neglect fluctuations with λ_⊥ ≫ ρ or λ_⊥ ≪ ρ, as noted above. Although we are interested in stochastic ion orbits, we can still define an effective guiding-center position,

$$\begin{equation} {\bf R} = {\bf r} + \frac{{\bf v} \times \hat{\bf {b}}}{\Omega }, \end{equation} \tag{ 10 }$$

where $\hat{\bf {b}} = {\bf B}/B$ and r is the ion's instantaneous position. When ε ≪ 1, the particle gyrates smoothly about position R. As ε increases toward 1, the particle's motion becomes more complicated, but the particle remains within a distance ∼ρ of position R. Taking the time derivative of Equation (10) and using the equation dv/dt = (q/m)(E + v × B/c), we obtain the equation

$$\begin{equation} \frac{d{\bf R}}{dt} = v_\parallel \hat{\bf {b}} + \frac{c{\bf E} \times {\bf B}}{B^2} + \cdots, \end{equation} \tag{ 11 }$$

where the ellipsis (⋅⋅⋅) represents terms proportional to derivatives of B, which we ignore in our approximate treatment. During a single cyclotron period, an ion passes through a small number of uncorrelated fluctuations or "structures" of transverse scale ∼ρ. Within different structures, the vector cE × B/B² has a similar magnitude (∼δv_ρ) but points in different directions. The time average of cE × B/B² over a single cyclotron period is thus somewhat smaller than, but of order, δv_ρ. The time Δt required for an ion's guiding center to move a distance ρ is thus approximately

$$\begin{equation} \Delta t \sim \frac{\rho }{\delta v_\rho }. \end{equation} \tag{ 12 }$$

In writing Equation (12), we have assumed that the gyroscale fluctuations do not oscillate on a timescale ≪Δt, and we continue to make this assumption in the analysis to follow. Each time the particle moves a distance ρ perpendicular to B₀, it encounters different and uncorrelated gyroscale electromagnetic fields. Thus, dR/dt decorrelates after a time Δt, and the particle's guiding center undergoes a random walk in space with diffusion coefficient ∼ρ²/Δt.

Similarly, when ε is sufficiently large that the ion's motion becomes stochastic, the value of dH/dt decorrelates after a time Δt, and the particle undergoes a random walk in energy. In contrast, as shown in the Appendix, as ε → 0 the interaction between ions and gyroscale electrostatic-potential structures is not a Markov process; instead, changes in H are correlated over long times and to leading order in ε are reversible and bounded. Returning to the stochastic case, we define $\overline{ \partial \Phi /\partial t}$ to be the rms value of ∂Φ/∂t associated with fluctuations with λ_⊥ ∼ ρ. The rms change in H during a time Δt is then

$$\begin{equation} \Delta H \sim q\,\overline{ \frac{\partial \Phi }{\partial t}}\,\Delta t. \end{equation} \tag{ 13 }$$

An ion undergoing stochastic motion can gain kinetic energy in the same way as the particle illustrated in Figure 1. If the ion spends a time Δt localized within a flux tube of cross-sectional area ∼ρ² and length ∼|v_∥|Δt, it exits this flux tube in a random direction. Thus, if $\partial \tilde{\Phi }/\partial t$ is on average positive during this time interval within the flux tube, it does not follow that the ion will move to a region of larger $\tilde{\Phi }$ after a time Δt, where $\tilde{\Phi }$ is the electrostatic potential associated with fluctuations with λ_⊥ ∼ ρ. On the contrary, the change in $\tilde{\Phi }$ along the ion's path is only loosely correlated with the average change in $\tilde{\Phi }$ within the flux tube. As a result, the change in the ion's kinetic energy during a time Δt is of the same order of magnitude as the change in its total energy given in Equation (13).⁵ Because ∇Φ is nearly perpendicular to B, and because the ion's guiding center moves perpendicular to B by a distance of order λ_⊥ ∼ ρ during a time Δt, the ion's perpendicular kinetic energy K_⊥ = mv²_⊥/2 changes by an amount of order

$$\begin{equation} \Delta K_\perp \sim \Delta H \end{equation} \tag{ 14 }$$

during a time Δt. We discuss the parallel kinetic energy following Equation (25) below. We define an effective frequency ω_eff for gyroscale fluctuations through the equation

$$\begin{equation} \overline{ \frac{\partial \Phi }{\partial t}} = \omega _{\rm eff}\, \delta \Phi _\rho. \end{equation} \tag{ 15 }$$

For example, if the gyroscale fluctuations consist of waves with a single frequency ω, then ω_eff = ω. With the use of Equations (5) and (12), we can rewrite Equation (14) as

$$\begin{equation} \Delta K_\perp \sim m v_\perp \omega _{\rm eff} \rho. \end{equation} \tag{ 16 }$$

The kinetic-energy diffusion coefficient D_K ∼ (ΔK_⊥)²/Δt is then given by

$$\begin{equation} D_K \sim m^2 v_\perp ^2 \omega _{\rm eff}^2 \,\delta v_\rho \, \rho. \end{equation} \tag{ 17 }$$

When a single ion undergoes kinetic-energy diffusion, the ion has an equal likelihood of gaining or losing kinetic energy during each "random-walk step" of duration Δt. On the other hand, if a large population of ions undergoes kinetic-energy diffusion, and if the phase-space density f of ions is a monotonically decreasing function of K_⊥, then the average value of K_⊥ increases steadily in time. To distinguish between properties of individual particles and rms quantities within a distribution, we define v_⊥i to be the rms perpendicular velocity of the ions, which is related to the perpendicular ion temperature T_⊥ by the equation

$$\begin{equation} v_{\perp {\rm i}} = \sqrt{\frac{2 k_{B} T_\perp }{m}}. \end{equation} \tag{ 18 }$$

We also define the rms ion gyroradius,

$$\begin{equation} \rho _{\rm i} = \frac{v_{\perp {\rm i}}}{\Omega }. \end{equation} \tag{ 19 }$$

We define δv_i to be the rms amplitude of the fluctuating fluid velocity at λ_⊥ ∼ ρ_i, and we set

$$\begin{equation} \varepsilon _{\rm i} = \frac{\delta v_{\rm i}}{v_{\perp {\rm i}}}. \end{equation} \tag{ 20 }$$

For protons, we define δv_p (δB_p) to be the fluctuating fluid velocity (magnetic field) at λ_⊥ ∼ ρ_p, and define

$$\begin{equation} \varepsilon _{\rm p} = \frac{\delta v_{\rm p}}{v_{\perp {\rm p}}}. \end{equation} \tag{ 21 }$$

The timescale for the average value of K_⊥ in a distribution of ions to double is then roughly

$$\begin{equation} t_{\rm i} \sim \frac{m^2 v_{\perp {\rm i}}^4 }{D_{K{\rm i}}}, \end{equation} \tag{ 22 }$$

where D_Ki is the value of D_K for ions with v_⊥ = v_⊥i and ρ = ρ_i. The perpendicular ion heating rate per unit mass is then Q_⊥ ∼ v²_⊥i/t_i, or

$$\begin{equation} Q_\perp \sim \omega _{\rm eff,i}^2 \,\delta v_{\rm i}\, \rho _{\rm i}, \end{equation} \tag{ 23 }$$

where ω_eff,i is the value of ω_eff at ρ = ρ_i.

We now consider what determines the value of ω_eff in anisotropic AW or KAW turbulence. If the turbulence is driven at an "outer scale" L₀ that is ≫ρ, the advection or "sweeping" of structures with λ_⊥ ∼ ρ by the outer-scale velocity fluctuations leads to rapid time variations in Φ at a fixed point in space. On the other hand, these large-scale velocity fluctuations advect both the ions and the small-scale structures in the electric and magnetic fields. Thus, if one considers ions within a flux tube of radius ∼ρ and length ≪L₀, and if one transforms to a frame of reference moving with the average velocity of that flux tube, then the rapid time variations resulting from large-scale advection disappear. This indicates that large-scale sweeping does not control the rate of ion heating or the value of ω_eff in Equation (23). On the other hand, electrostatic-potential structures at scale λ_⊥ ≃ 0.5ρ are advected by velocity fluctuations at the same scale, and there is no frame of reference in which the velocities at λ_⊥ ≃ 0.5ρ vanish at all points along an ion's gyro-orbit. This advection by velocity fluctuations with λ_⊥ ≃ 0.5ρ causes $\overline{ \partial \Phi /\partial t}$ to have a value of ∼δΦ_ρδv_ρ/ρ, which gives

$$\begin{equation} \omega _{\rm eff} \sim \frac{\delta v_{\rho }}{\rho }, \end{equation} \tag{ 24 }$$

where we have neglected factors of order unity, such as the ratio between δv_ρ and the rms amplitude of the velocity fluctuation at λ_⊥ ≃ 0.5ρ. Put another way, the advection of electrostatic-potential structures at λ_⊥ ∼ 0.5ρ, which are rooted in the electron fluid, leads to a partial time derivative of Φ that ions can feel, and which energizes ions through the process illustrated in Figure 1. We note that in "imbalanced" (or cross-helical) AW turbulence, in which the majority of the waves propagate either parallel to B₀ or anti-parallel to B₀, the energy cascade time for the majority waves can greatly exceed ω⁻¹_eff, since the majority waves are cascaded by the smaller-amplitude waves propagating in the opposite direction. Nevertheless, even for imbalanced turbulence, the arguments leading to Equation (24) continue to hold.

As discussed following Equation (7) and in the Appendix, when ε is sufficiently small, the changes in H remain correlated (and largely reversible) over long times, so that the perpendicular heating rate is strongly reduced relative to our estimate in Equation (23). To account for this, we introduce a multiplicative suppression factor onto the right-hand side of Equation (23) of the form exp(−c₂/ε_i). We also add an overall coefficient c₁ to the right-hand side of Equation (23) to account for the various approximations we have made. Both c₁ and c₂ are dimensionless constants whose values depend upon the nature of the fluctuations (e.g., whether the fluctuations are waves or turbulence, the type of turbulence, etc.) and the shape of the ion velocity distribution. Substituting Equation (24) into Equation (23), we obtain

$$\begin{equation} Q_\perp = \frac{c_1(\delta v_{\rm i})^3}{\rho _{\rm i}} \,\exp \left(-\,\frac{c_2}{\varepsilon _{\rm i}}\right). \end{equation} \tag{ 25 }$$

We emphasize that for protons in low-β plasmas ε_p ≃ β^−1/2δB_p/B₀ and thus ε_p can approach unity even if δB_p/B₀ remains ≪1.

We have restricted our analysis to AWs and KAWs with λ_⊥ ∼ ρ ≳ ρ_p and ω⁻¹ ≳ Δt ∼ λ_⊥/δv_ρ. This condition on the wave frequency implies that the parallel wavelengths of such fluctuations satisfy the inequality λ_∥ ≳ ρv_A/δv_ρ ≫ λ_⊥. If β < 1, then λ_∥>v_∥Δt, and the change in an ion's parallel kinetic energy K_∥ = mv²_∥/2 during a time Δt due to the parallel electric field E_∥ is ΔK_∥ ∼ qE_∥v_∥Δt. When m_e/m_p < β < 1 and λ_⊥>ρ_p, E_∥/E_⊥ ∼ ρ²_p/(λ_⊥λ_∥), where $m_{\rm {e}}$ is the electron mass (Hollweg 1999) and E_⊥ is the electric-field component perpendicular to B. Thus, ΔK_∥ is ≲v_∥/v_A times the value of ΔK_⊥ in Equation (16). For thermal ions in low-β plasmas, v_∥ ≪ v_A. Thus, when ε_i is sufficiently large that stochastic heating is important, stochastic heating leads primarily to perpendicular ion heating rather than parallel heating. For AWs/KAWs in low-β plasmas, the parallel component of the magnetic mirror force is much less than qE_∥ (Hollweg 1999) and thus does not affect our conclusions regarding anisotropic heating at β ≪ 1.

2.1. Test-particle Simulations of Proton Heating

To test the above ideas, we have numerically simulated test-particle protons interacting with a spectrum of randomly phased KAWs. The protons' initial locations are chosen randomly from a uniform distribution within a volume encompassing many wavelengths perpendicular and parallel to B₀. The protons' initial velocities are drawn randomly from an isotropic Maxwellian distribution of temperature T_p. For each particle, we solve the equations

$$\begin{equation} \frac{d {\bf x} }{dt} = {\bf v} \end{equation} \tag{ 26 }$$

and

$$\begin{equation} \frac{d{\bf v}}{dt} = \frac{q}{m}\left({\bf E} + \frac{{\bf v} \times {\bf B} }{c}\right) \end{equation} \tag{ 27 }$$

using the Bulirsch–Stoer method (Press et al. 1992). We take ${\bf B} = B_0 \hat{\bf {z}} + {\bf B}_1$, where B₀ is constant. We take E and B₁ to be the sum of the electric and magnetic fields from 162 waves with randomly chosen initial phases, with two waves at each of 81 different wavevectors. At each wavevector, there is one wave with ω/k_z>0 and a second wave with ω/k_z < 0. This second wave has the same amplitude as the first, so that there are equal fluxes of waves propagating in the +z- and −z-directions. The 81 different wavevectors consist of nine wavevectors at each of nine different values of k_⊥, denoted k_⊥j. The k_⊥j can be expressed in terms of ρ_p. In particular, the values ψ_j = ln(k_⊥jρ_p) are uniformly spaced between −4/3 and 4/3; i.e., ψ_j = −4/3 + j/3, with j = 0, 1, ..., 8. We regard the values ψ_j as corresponding to cell centers in a uniform grid in ψ = ln(k_⊥ρ_p), with grid spacing Δψ = 1/3. The middle three grid cells, with j = 3, 4,  and 5, thus correspond to an interval of width unity in ln(k_⊥) space centered on k_⊥ρ_p = 1. We define the rms amplitudes of the gyroscale velocity and magnetic-field fluctuations δv_p and δB_p in our simulations by taking the rms values of the E × B velocity and magnetic-field fluctuation resulting from the KAWs in these middle three grid cells. At each k_⊥j we include nine different values of the azimuthal angle ϕ in k space, ϕ_l = 2πl/9, where l = 0, 1, ..., 8. At each k_⊥j there is only a single value of k_∥, which we denote k_∥j. We choose k_∥4 so that the frequency at k_⊥ = k_⊥4 and k_∥ = k_∥4 equals k_⊥4δv_ρ. The linear frequency of our gyroscale KAWs is thus comparable to the value of ω_eff given in Equation (24) for KAW turbulence at k_⊥ρ_p ∼ 1. We then set

$$\begin{equation} \frac{k_{\parallel j}}{k_{\parallel 4}} = \left\lbrace \begin{array}{@{}ll} (k_{\perp j}/k_{\perp 4})^{2/3} & {\rm if}\;0 \le j<4 \\ \ms (k_{\perp j}/k_{\perp 4})^{1/3} & {\rm if}\;4 < j \le 8 \end{array}. \right. \end{equation} \tag{ 28 }$$

Our formula for k_∥j at j < 4 is chosen so that the wave periods are comparable to the energy cascade timescales in the critical-balance theory of Goldreich & Sridhar (1995), while the formula for j>4 is chosen so that the wave periods match the nonlinear timescales in the critical-balance theory of Cho & Lazarian (2004). All waves at the same k_⊥ have the same amplitude, and (since there are the same number of waves at each k_⊥j) we take the amplitude of the magnetic-field fluctuation in each wave to be ∝k^−1/3_⊥ for k_⊥ρ_p < 1 and ∝k^−2/3_⊥ for k_⊥ρ_p>1, again motivated by the theories of Goldreich & Sridhar (1995) and Cho & Lazarian (2004).

The relative amplitudes of the different components of E and B₁ for each wave are taken from the two-fluid theory of Hollweg (1999). To apply this theory, we choose plasma parameters that are characteristic of coronal holes. In particular, we set β_e = 8πnk_BT_e/B²₀ = 0.003, v_A = 0.003c, and T_e = 0.5T_p, where n is the electron number density (equal to the proton number density), and T_e is the electron temperature.

Using the above procedures, we have carried out seven simulations with different values for the overall normalization of the wave amplitudes, with δB_p/B₀ ranging from 4.8 × 10⁻³ to 1.9 × 10⁻². Given the polarization properties of KAWs and our method for constructing the wave spectra, the value of δv_p/v_A is 1.19 times the value of δB_p/B₀ in each simulation. The wave frequencies reach their maximum values in the largest δB_p/B₀ simulation. In this simulation, ω = 0.29Ω_p at k_⊥ρ_p = 1, and ω = 0.82Ω_p at the maximum value of k_⊥ρ_p, which is 3.79. Although this maximum frequency is close to Ω_p, the cyclotron resonance condition ω − k_∥v_∥ = lΩ_p (where l is any integer) is not satisfied by thermal protons, because the parallel thermal speed of the protons is only 0.055v_A and |k_∥v_∥/ω| ⩽ |v_∥|/v_A. For most of the waves in these simulations, ω ≪ Ω_p.

We determine the perpendicular proton heating rate per unit mass Q_⊥p in the simulations by plotting 〈v²_⊥〉 versus time, fitting this plot to a straight line to determine (d/dt)〈v²_⊥〉, and then setting Q_⊥p = 0.5(d/dt)〈v²_⊥〉, where 〈...〉 indicates an average over the 10³ particles in each simulation. When fitting the plot of 〈v²_⊥〉 versus time, we ignore the first 10 cyclotron periods, because during the first couple of gyroperiods the particles undergo a modest apparent heating as they "pick up" some portion of the E × B velocity of the waves. We find that after 〈v²_⊥〉 increases by between 20% and 40%, the heating rate starts to decrease for two reasons. First, the small-v_⊥ part of the velocity distribution flattens, after which this part of the distribution is no longer heated as effectively. Second, as 〈v²_⊥〉 increases, ε_p decreases. We neglect this later stage of weaker heating when constructing our fits to the 〈v²_⊥(t)〉 plots, so that the measured heating rates correspond to Maxwellian distributions. (For the smallest values of δv_p, we do not observe a second stage of weaker heating, because the test-particle velocity distributions do not change very much during the simulations, which last 10⁴Ω⁻¹_p.) We illustrate this procedure in Figure 2 for a run with δv_p/v_⊥p = 0.15. In this case, we determine Q_⊥p from the slope of the long-dashed line.

Zoom In Zoom Out Reset image size

In Figure 3, we plot the values of Q_⊥p for several different values of ε_p. Each × in this figure corresponds to a separate simulation with a different value of δv_p but the same initial proton temperature. The solid line is the proton heating rate from Equation (25) with c₁ = 0.75 and c₂ = 0.34; that is,

$$\begin{equation} Q_{\perp {\rm p}} = \frac{0.75 (\delta v_{\rm p})^3}{\rho _{\rm p}} \,\exp \left(-\,\frac{0.34}{\varepsilon _{\rm p}}\right). \end{equation} \tag{ 29 }$$

We expect the constants c₁ and c₂ to be fairly insensitive to variations in β_e, T_p/T_e, and v_A/c, at least within the range of parameters relevant to low-β solar-wind streams, in which case Q_⊥p depends on the plasma parameters primarily through the explicit ρ_p and ε_p terms in Equation (29). We have partially tested this hypothesis by carrying out seven additional simulations, corresponding to three additional choices of plasma parameters, which we denote "cases 1, 2, and 3." In case 1, v_A = 0.003c and T_e = 0.5T_p as in the original simulations, but β_e is increased by a factor of 5 to the value 1.5 × 10⁻². In case 2, v_A = 0.003c as in the original simulations, but β_e is increased by a factor of 4 to the value 1.2 × 10⁻², and T_p/T_e is reduced by a factor of 4 to the value 0.5. In case 3, β_e = 0.003 and T_e = 0.5T_p as in the original simulations, but v_A/c is increased by a factor of 10 to the value 0.03. In all of these additional simulations, the number of waves, the values of ψ_j, the procedure for determining k_∥j, and the wavenumber scaling of the wave amplitudes are the same as in the original simulations, and 500 test particles are used. Figure 4 shows the perpendicular proton heating rates in these additional simulations, as well as the value of $Q_{\perp \rm p}$ in Equation (29) for the three choices of plasma parameters. Equation (29) agrees reasonably well with the simulation results for all three cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The values of c₁ and c₂ in Equation (29) presuppose the presence of a broad spectrum of AWs and KAWs bracketing the perpendicular wavenumber k_⊥ = (ρ_p)⁻¹, encompassing at a minimum the range 0.3 ≲ k_⊥ρ_p ≲ 3. A spectrum of at least this width is probably present in the solar wind, the only uncertainty being the value of the dissipation wavenumber beyond which the wave power spectrum decreases exponentially with increasing k_⊥. If the simulations described in this section are repeated without the smallest three values of k_⊥ and without the largest three values of k_⊥ (keeping the wave amplitudes fixed at the middle three values of k_⊥), then the proton orbits become less stochastic, and Q_⊥p decreases significantly. (The exact amount by which Q_⊥p decreases depends upon the value of ε_p.) We have omitted waves at k_⊥ρ_p < 0.26 and k_⊥ρ_p>3.8, but we expect that waves at such scales do not have a strong effect on perpendicular ion heating in the solar wind, provided ω is sufficiently small that the cyclotron resonance condition can not be satisfied. It is possible that in some cases strongly turbulent fluctuations with k_⊥ρ_p ≫ 1 and nonlinear timescales ∼Ω⁻¹_p could heat ions through a broadened cyclotron resonance, but a detailed investigation of this process is beyond the scope of this study.

We reiterate that the values of c₁ and c₂ in Equation (29) are not universal, but instead depend on the type of fluctuations that are present. In strong turbulence (as opposed to randomly phased waves), the value of c₂ may be smaller than in our simulations (indicating stronger heating), because a significant fraction of the cascade power may be dissipated in coherent structures in which the fluctuating fields are larger than their rms values (Dmitruk et al. 2004). Other differences between AW turbulence and the wave fields in our simulations may also affect the values of c₁ and c₂. For example, in AW turbulence, the electric and magnetic fields vary most rapidly in directions perpendicular to the local magnetic field (Cho & Vishniac 2000). In contrast, in our simulations, the wave fields vary most rapidly in the direction perpendicular to the global mean field, which differs from the local magnetic field. Also, in AW turbulence, the E × B velocity associated with the large-scale fluctuations advects not just the particles but also the small-scale AW fluctuations, an effect that is not captured in our simulations. Because of these differences, the values of c₁ and c₂ that describe our test-particle simulations provide only a preliminary estimate of the values corresponding to solar-wind turbulence.

The lower solid-line curve in Figure 2 is a plot of 〈v²_∥〉 versus time in the simulation with ε_p = 0.15, β_e = 0.003, v_A = 0.003c, and T_e = 0.5T_p. During the interval 10 < Ω_pt < 2200, the increase in 〈v²_∥〉 is about one-fourth the increase in 〈v²_⊥〉. However, most of the increase in 〈v²_∥〉 is an artifact of our numerical method, which equates the parallel electric fields of the waves with the z-component of the electric field in the simulation, and the perpendicular electric field of the waves with the x- and y-components of the electric field in the simulation. The local magnetic field in our simulations, however, is not parallel to the z-axis, but instead has nonzero x- and y-components resulting from the magnetic-field fluctuations. As a result, part of the perpendicular wave electric field is converted into a parallel electric field in the simulation, artificially enhancing the parallel electric field seen by the particles. To eliminate this effect, we have repeated this simulation replacing the local electric field E seen by each particle with the adjusted electric field ${\bf E}^\prime = {\bf E} + \hat{\bf {b}} (E_z - \hat{\bf {b}} \,{\bm \cdot}\, {\bf E})$, where $\hat{\bf {b}}= {\bf B}/B$ and B is the local value of the magnetic field. In this new simulation, the parallel electric field $\hat{\bf {b}} \cdot {\bf E}^\prime$ is the sum of the parallel electric fields of the individual waves in the simulation and does not include any contribution from the perpendicular electric fields of the individual waves. The value of 〈v²_∥〉 in this modified simulation, shown as a dashed line in Figure 2, does not increase significantly during the course of the simulation (in fact it decreases slightly), consistent with our argument above that parallel heating is weak when β ≪ 1.

2.2. Proton Heating at k_⊥ρ_p ∼ 1 as a Fraction of the Turbulent Cascade Power

The cascade power per unit mass at k_⊥ρ_p ∼ 1, which we denote Γ, depends upon whether the turbulence is "balanced" or "imbalanced," where balanced (imbalanced) turbulence involves equal (unequal) fluxes of waves propagating parallel to B₀ and anti-parallel to B₀. In balanced KAW turbulence,

$$\begin{equation} \Gamma = C_{K}^{-3/2} \left(\frac{ \delta v_{\rm p}}{\rho _{\rm p}}\right) \left(\frac{\delta B_{\rm p}}{B_0}\right)^2 v_{A}^2, \end{equation} \tag{ 30 }$$

where C_K is a dimensionless constant (Howes et al. 2008a). It can be inferred from the numerical simulations of Howes et al. (2008b) that C_K = 2.0 (G. Howes 2009, private communication). In the simulations of Section 2.1, δB_p/B₀ = 0.84δv_p/v_A, and we make the approximation that this same ratio is characteristic of KAW turbulence in general. Combining Equations (29) and (30), we obtain

$$\begin{equation} \frac{Q_{\perp {\rm p}} }{\Gamma } = 3.0 \,\exp \left({-}\frac{0.34}{\varepsilon _{\rm p}}\right). \end{equation} \tag{ 31 }$$

We expect that C_K depends only weakly on β, T_p/T_e, and v_A/c, so that the numerical constants 3.0 and 0.34 in Equation (31) are relatively insensitive to the plasma parameters, at least within the range of parameters found in low-β solar-wind streams. Equation (31) implies that perpendicular proton heating by KAWs with k_⊥ρ_p ∼ 1 absorbs ⩾1/2 of the cascade power at k_⊥ρ_p ∼ 1 when ε_p exceeds

$$\begin{equation} \varepsilon _{\rm crit} = 0.19 \;. \end{equation} \tag{ 32 }$$

The cascade power in imbalanced AW turbulence is smaller than in balanced AW turbulence with the same total fluctuation energy, because the AW energy cascade requires interactions between oppositely propagating waves (Iroshnikov 1963; Kraichnan 1965). At k_⊥ρ_p ∼ 1, KAWs propagating in the same direction can interact nonlinearly with one another, but the importance of such interactions relative to interactions between oppositely propagating waves is not well known. Despite this uncertainty, we expect that if AW/KAW turbulence is imbalanced at k_⊥ρ_p ∼ 1, then the cascade power at k_⊥ρ_p ∼ 1 is less than in Equation (30). On the other hand, it is unlikely that imbalance strongly affects Q_⊥p if δv_p is held fixed (except for particles with v_∥ ∼ ±v_A, as discussed in Section 2.4). We thus expect perpendicular proton heating to absorb at least 50% of the cascade power at k_⊥ρ_p ∼ 1 in imbalanced turbulence when ε_p exceeds a value that is smaller than 0.19.

2.3. Proton Heating versus Electron Heating by KAWs with k_⊥ρ_p ∼ 1

Stochastic proton heating removes energy from KAW fluctuations with k_⊥ρ_p ∼ 1, resulting in an effective damping rate for these fluctuations, which we denote γ_p. The value of γ_p is given by the relation

$$\begin{equation} 2 \gamma _{\rm p} {\cal E}_{\rm w} = Q_{\perp {\rm p}}, \end{equation} \tag{ 33 }$$

where ${\cal E}_{\rm w}$ is the energy per unit mass of the KAW fluctuations at k_⊥ρ_p ∼ 1. The factor 2 in Equation (33) is included to make γ_p analogous to a linear wave damping rate, in the sense that the rate at which linear waves lose energy is twice the product of the damping rate and the wave energy. To estimate the value of γ_p in AW/KAW turbulence, we use the test-particle calculations in Section 2.1 for a spectrum of randomly phased KAWs. We take ${\cal E}_{\rm w}$ to be the energy per unit mass of the full spectrum of waves in these simulations. (This choice leads to a conservative estimate of γ_p, since the damping is likely concentrated in the subset of the waves with k_⊥ρ_p ≳ 1.) On the other hand, we continue to define (δv_p)² as the mean-square E × B velocity associated with KAWs with values of k_⊥ lying within a logarithmic interval of width unity centered on k_⊥ρ_p = 1. With these definitions, ${\cal E}_{\rm w} = 2.1 (\delta v_{\rm p})^2$ in all of the simulations in Section 2.1. Combining Equations (29) and (33), we obtain

$$\begin{equation} \gamma _{\rm p } = 0.18 \varepsilon _{\rm p} \Omega _{\rm p} \,\exp \left({-}\frac{0.34}{\varepsilon _{\rm p}}\right). \end{equation} \tag{ 34 }$$

In low-β plasmas, small-amplitude KAWs with k_⊥ρ_p = 1 and ω ≪ Ω_p undergo electron Landau damping but negligible linear proton damping (Quataert 1998; Gruzinov 1998; Gary & Nishimura 2004). Using the numerical method described by Quataert (1998) and Howes et al. (2008a), we numerically solve the full hot-plasma dispersion relation to find the electron damping rate γ_e of KAWs with k_⊥ρ_p = 1 and ω ≪ Ω_p for a range of values of k_∥, T_p/T_e, β_p, and v_A/c, where β_p = 8πnk_BT_p/B²₀. We find that if m_e/m_p ≪ β_e ≪ 1, v_A ≪ c, and 0.1 ≲ T_p/T_e ≲ 10, then the damping rate at k_⊥ρ_p = 1 is well fit by the formula γ_e = 9.5 × 10⁻³(T_e/T_p)^1/2β^−1/2_p|k_∥v_A|, or equivalently

$$\begin{equation} \gamma _{\rm e} = 9.5 \times 10^{-3}\varepsilon _{\rm p} \chi ^{-1} \left(\frac{ T_{\rm e}}{\beta _{\rm p}T_{\rm p}}\right)^{1/2} \Omega _{\rm p}, \end{equation} \tag{ 35 }$$

where χ ≡ k_⊥δv_p/|k_∥v_A|. In some theories of strong MHD turbulence χ ∼ 1 (Goldreich & Sridhar 1995; Boldyrev 2006). This condition, sometimes referred to as critical balance, may characterize AW/KAW fluctuations in coronal holes and the solar wind at k_⊥ρ_p ∼ 1. On the other hand, if the frequencies of the waves launched by photospheric motions are sufficiently small, then AW/KAW turbulence at a heliocentric distance of a few solar radii may be more "two-dimensional" than in critical-balance models, with smaller values of k_∥ and a larger value of χ.

Combining Equations (34) and (35), we obtain

$$\begin{equation} \frac{\gamma _{\rm p}}{\gamma _{\rm e}} = 19\, \chi \, \left(\frac{\beta _{\rm p}T_{\rm p}}{ T_{\rm e}}\right)^{1/2} \exp \left(-\,\frac{0.34}{\varepsilon _{\rm p}}\right). \end{equation} \tag{ 36 }$$

The ratio γ_p/γ_e approximates the ratio of the proton heating rate to the electron heating rate resulting from KAW fluctuations at k_⊥ρ_p ∼ 1 in the low-β conditions present in coronal holes and the near-Sun solar wind. (At β ≳ 1, linear KAW damping on the protons becomes important, increasing the proton heating rate.) We note that if the damping timescales γ⁻¹_p and γ⁻¹_e are both much longer than the energy cascade time at k_⊥ρ_p ∼ 1, then most of the fluctuation energy will cascade past the proton-gyroradius scale to smaller scales.

2.4. How the Heating Rate Depends on q, m, β, and v_∥/v_A

If we re-run our simulations, keeping only waves with ω/k_z>0, and consider a thermal distribution of test-particle protons with a nonzero average velocity equal to $v_{A} \hat{z}$, then the perpendicular heating rate is strongly reduced. This is because the electric field of an AW (or KAW with λ_⊥ ∼ ρ_p) vanishes (or is strongly reduced) in a reference frame moving at speed v_A in the same direction as the wave along the background magnetic field. This effect may explain the observation that the perpendicular heating of α particles in the solar wind is reduced when the differential flow velocity of α particles relative to protons (in the anti-Sunward direction) approaches v_A (Kasper et al. 2008), at least in regions where anti-Sunward propagating KAWs dominate over Sunward-propagating KAWs.

If we hold δB_p/B₀ fixed but increase β_p to 1, then the perpendicular proton heating rate is dramatically reduced, because ε_p = δv_p/v_⊥ ∼ β^−1/2δB_p/B₀ decreases by a large factor. On the other hand, the protons in these β_p ∼ 1 simulations undergo significant parallel heating, consistent with results from linear theory (Quataert 1998) and recent test-particle simulations of ions propagating in numerically simulated MHD turbulence (Lehe et al. 2009).

If we re-run our simulations but use O⁺⁵ ions instead of protons (but with the same temperature as the protons), then the perpendicular heating rate is much larger. This is in large part because ε is larger for O⁺⁵ (and other heavy ions) than for protons at the same temperature, a point to which we return in Section 3. Another reason for enhanced heavy-ion heating can be seen from Equation (22). We rewrite this equation with the aid of Equation (24), increasing t_i by exp(c₂/ε_i) for the same reasons that we reduced Q_⊥ by this same factor in Equation (25), to obtain

$$\begin{equation} t_{\rm i} \sim \frac{v_{\perp {\rm i}}^2 \rho _{\rm i}}{(\delta v_{\rm i})^3} \, \exp \left(\frac{c_2}{\varepsilon _{\rm i}}\right). \end{equation} \tag{ 37 }$$

In a number of theories of MHD turbulence, the ratio (δv_ρi)³/ρ_i is relatively (or completely) insensitive to the value of ρ_i, provided ρ_i is in the inertial range of the turbulence. On the other hand, for ion species at equal temperatures, v²_⊥i is inversely proportional to the ion mass. Thus, even aside from the exponential factor in Equation (37), the heating timescale is shorter for heavier ions than for lighter ions at the same temperature.

Finally, if we repeat the simulations of Section 2.1 for test-particle ions with ρ_i ≫ ρ_p, and with values of k_⊥ρ_i centered on 1 so that the gyroscale fluctuations are now AWs, we recover similar values for the perpendicular heating rate per unit mass. Stochastic perpendicular ion heating thus does not require the particular polarization properties of KAWs, but operates for both KAWs and AWs, as we have argued in our heuristic derivation of Equation (25).

2.5. Lack of Perpendicular Heating by AWs with k_⊥ρ ≪ 1

As shown in the previous section, the stochastic ion heating rate is a strongly increasing function of $\varepsilon _{\rm i} = \delta v_{\rm i}/v_{\perp \rm i}$. For fixed turbulence properties, the value of ε_i depends upon the ion charge q = Ze, the ion mass m = Am_p, and the perpendicular ion temperature T_⊥. For example, if we take the rms amplitude of the turbulent velocity fluctuation at transverse scale λ_⊥ to be given by

$$\begin{equation} \delta v_{\lambda _\perp } = \alpha v_{A} \left(\frac{\lambda _\perp }{L_0}\right)^a \end{equation} \tag{ 38 }$$

for ρ_i < λ_⊥ < L₀, where α and a are dimensionless constants and L₀ is the outer scale or driving scale of the turbulence, then

$$\begin{equation} \varepsilon _{\rm i} = \alpha \left(\frac{T_{\rm p}}{T_\perp \mu _{\rm p} \beta _{\rm p}}\right)^{(1-a)/2} \frac{A^{(1+a)/2}}{Z^a}\left(\frac{d_{\rm p}}{L_0}\right)^a, \end{equation} \tag{ 39 }$$

where d_p = v_A/Ω_p is the proton inertial length, β_p = 8πn_pk_BT_p/B², n_p is the proton density, T_p is the proton temperature, and μ_p is the mean molecular weight per proton; that is, the mass density is μ_pn_pm_p, and the Alfvén speed is $B/\sqrt{4\pi \mu _{\rm p} n_{\rm p} m_{\rm p}}$. If the velocity power spectrum P^(v)_k is $\propto k_\perp ^{-c_3}$ for L⁻¹₀ < k_⊥ < ρ⁻¹_i, then

$$\begin{equation} a = \frac{c_3 - 1}{2}. \end{equation} \tag{ 40 }$$

To investigate the possible role of stochastic ion heating in coronal holes and the fast solar wind, we evaluate Equation (39) as a function of heliocentric distance r using observationally constrained profiles for the density, temperature, and field strength. We take n_p to be given by Equation (4) of Feldman et al. (1997), which describes coronal holes out to several solar radii, plus an additional component proportional to r⁻²:

$$\begin{equation} n_{\rm p}(r) = \left(\frac{3.23 \times 10^8}{d^{15.6}} + \frac{2.51 \times 10^6}{d^{3.76}} + \frac{1.85 \times 10^5}{d^2}\right)\;{\rm cm}^{-3}, \end{equation} \tag{ 41 }$$

where d = r/R_☉. Equation (41) gives n_p = 4 cm⁻³ at 1 AU. We set

$$\begin{equation} T_{\rm p } = 3\times 10^6\;{\rm K} \cdot \left[\frac{ 1 - (2/3) \exp (-d/1.5)}{(1 + 0.1d)^{0.8}}\right], \end{equation} \tag{ 42 }$$

which leads to a proton temperature that is ∼2 × 10⁶ K at r ≲ 2 R_☉, and ∼2.5 × 10⁵ K at 1 AU. We take the magnetic field strength to be (Hollweg & Isenberg 2002)

$$\begin{equation} B_0 = \left[\frac{1.5(f_{\rm max} - 1)}{d^6} + \frac{1.5}{d^2}\right]\;{\rm Gauss}, \end{equation} \tag{ 43 }$$

with f_max (the super-radial expansion factor) equal to 5. We determine the rms amplitude of the fluctuating wave velocity at the outer scale, $\delta v_{L_0} = \alpha v_{A}$, as a function of r using the analytical model of Chandran & Hollweg (2009), which describes the propagation of low-frequency AWs launched outward from the Sun, taking into account both non-WKB reflection and the cascade and dissipation of wave energy. In particular, we set $\delta v_{L_0}$ equal to the value of δv_rms plotted with a solid line in Figure 6 of Chandran & Hollweg (2009; the curve corresponding to their "extended model"). We take L₀ to be 10⁴ km at the coronal base (the limit d → 1 in Equation (43)), and to be proportional to B^−1/2.

We consider three different values for the spectral index c₃: 5/3, 3/2, and 6/5. The value c₃ = 5/3 is suggested by in situ measurements of magnetic-field fluctuations in the solar wind (Matthaeus & Goldstein 1982; Bruno & Carbone 2005), as well as some theoretical and numerical studies of MHD turbulence (Goldreich & Sridhar 1995; Cho & Vishniac 2000; Cho et al. 2002; Müller & Biskamp 2000; Haugen et al. 2004). The value c₃ = 3/2 is motivated by a different set of theoretical and numerical studies (Maron & Goldreich 2001; Müller & Grappin 2005; Boldyrev 2006; Mason et al. 2008; Perez & Boldyrev 2009), as well as recent in situ observations of the velocity power spectrum in the solar wind (Podesta et al. 2007; Podesta & Bhattacharjee 2010). The third value, c₃ = 1.2, follows from recent numerical simulations of reflection-driven AW turbulence in coronal holes and the fast solar wind (Verdini et al. 2009). In these simulations, c₃ = 1.2 at r < 1.2 R_☉, and c₃ gradually increases toward 5/3 with increasing r.

In Figure 5, we plot ε_i for H⁺, He⁺⁺, and O⁺⁵ assuming μ_p = 1.2. Although alpha particles and minor ions are observed to be hotter than protons in the fast solar wind, we have set all the ion temperatures equal to T_p to investigate the relative heating rates of different ion species that start out at the same temperature. Figure 5 illustrates the general point that ε_i depends strongly on the spectral index c₃. In particular, decreasing c₃ by 28% from 5/3 to 1.2 increases ε_i by a factor of >10 at all radii shown for all three ion species. Because Q_⊥ depends strongly on ε_i, Q_⊥ is extremely sensitive to the value of c₃. A second general point illustrated by Figure 5 is that when c₃ is fixed, ε_i depends only weakly on r for 2 R_☉ < r < 1 AU. As a result, given our assumptions, a large radial variation in ε within this range of r requires a radial variation in the spectral index c₃. A third point illustrated by Figure 5 is that ε_i is significantly larger for He⁺⁺ and O⁺⁵ than for protons at the same temperature. For example, at equal temperatures, protons and alpha particles have the same gyroradius, and ε_α = 2ε_p, where ε_α is the value of ε_i for alpha particles. Because of the strong dependence of Q_⊥ on ε_i, the stochastic heating rate per unit volume for alpha particles may in some cases be comparable to the corresponding rate for protons, even though Helium comprises only ∼20% of the mass in the solar wind. In addition, the comparatively large value of ε_i for O⁺⁵ may explain why O⁺⁵ ions are observed to be so much hotter than protons in the coronal holes (Kohl et al. 1998; Li et al. 1998; Esser et al. 1999; Antonucci et al. 2000), and likewise for other minor ions.

Zoom In Zoom Out Reset image size

When an ion interacts with turbulent AWs and/or KAWs, and when the amplitudes of the fluctuating electromagnetic fields at λ_⊥ ∼ ρ are sufficiently large, the ion's orbit becomes chaotic, and the ion undergoes stochastic perpendicular heating. When β ≲ 1, the parameter that has the largest effect on the heating rate is ε = δv_ρ/v_⊥, where δv_ρ is the rms amplitude of the velocity fluctuation at λ_⊥ ∼ ρ. In the limit ε → 0, the ion's magnetic moment is nearly conserved, and perpendicular ion heating is extremely weak. On the other hand, as ε increases toward unity, magnetic moment conservation is violated, and stochastic perpendicular heating becomes increasingly strong.

Using phenomenological arguments, we derive an analytic formula for the perpendicular heating rate Q_⊥ for different ion species in plasmas with β ≲ 1. This formula (Equation (25)) contains two dimensionless constants, c₁ and c₂, whose values depend on the nature of the fluctuations (e.g., waves versus turbulence, the slope of the power spectrum) and the shape of the ion velocity distribution. Using test-particle simulations, we numerically evaluate these constants for the case in which a Maxwellian distribution of protons interacts with a spectrum of random-phase AWs and KAWs at perpendicular wavenumbers in the range 0.264 < k_⊥ρ_p < 3.79, where ρ_p is the rms proton gyroradius in the background magnetic field B₀. The particular form of the wave power spectrum that we choose for these simulations is motivated by the critical-balance theories of Goldreich & Sridhar (1995) and Cho & Lazarian (2004). For this case, c₁ = 0.75 and c₂ = 0.34. The proton heating rate Q_⊥p can be compared to the cascade power Γ that would be present at k_⊥ρ_p ∼ 1 in balanced (see Section 2.2) AW/KAW turbulence with the same value of δv_p. When c₁ = 0.75 and c₂ = 0.34, the ratio Q_⊥p/Γ exceeds 1/2 when ε_p>ε_crit = 0.19, where ε_p is the value of ε for thermal protons.

Our expression for Q_⊥p/Γ (Equation (31)) may differ from the value of Q_⊥p/Γ in the solar wind for two main reasons. First, in strong AW/KAW turbulence (as opposed to randomly phased waves), a significant fraction of the cascade power may be dissipated in coherent structures in which the fluctuating fields are larger than their rms values (Dmitruk et al. 2004). Proton orbits in the vicinity of such structures are more stochastic than in average regions, and thus c₂ may be smaller in AW/KAW turbulence than in our test-particle simulations, indicating stronger heating. The perpendicular heating rate is very sensitive to the value of c₂/ε_p; our test-particle simulations are consistent with Q_⊥p/Γ being ∝exp(−c₂/ε_p). Thus, decreasing c₂ leads to a large increase in Q_⊥p/Γ when ε_p < c₂. Decreasing c₂ also decreases ε_crit, the value of ε_p at which Q_⊥p/Γ = 1/2 in balanced turbulence; it follows from Equations (25) and (30) that if C_K ≃ 2, δB_p/B₀ ≃ 0.84δv_p/v_A, and c₁ ≃ 1, then ε_crit ≃ c₂/2. The second reason that our estimate of Q_⊥p/Γ may differ from the value in the solar wind is that our formula for Γ does not take into account imbalance (see Section 2.2). If the turbulence is imbalanced at k_⊥ρ_p ∼ 1, then Γ is reduced relative to the balanced case. Thus, in imbalanced turbulence, $Q_{\perp \rm p}/\Gamma = 1/2$ at a value of ε_p that is smaller than ε_crit. We define ε* to be the value of ε at which $Q_{\perp \rm p}/\Gamma = 1/2$ regardless of whether the turbulence is balanced or not; thus, ε* = ε_crit for balanced turbulence, and ε* < ε_crit for imbalanced turbulence.

When β ≪ 1, stochastic proton heating by AW/KAW turbulence at k_⊥ρ_p ∼ 1 increases T_⊥ much more than T_∥. In contrast, linear proton damping of KAWs with ω ≪ Ω_p and k_⊥ρ_p ∼ 1 leads almost entirely to parallel heating, and is only significant when the proton thermal speed is ≳v_A; i.e., when β_p ≳ 1 (Quataert 1998). If we assume that (nonlinear) stochastic heating and linear wave damping are the only dissipation mechanisms for low-frequency AW/KAW turbulence, then we arrive at the following conclusions about how the cascade power in AW/KAW turbulence is partitioned between parallel and perpendicular heating, and between protons and electrons.

• 1.  
  If β_p ≪ 1 and ε_p ≪ ε*, then proton heating is negligible and electrons absorb most of the cascade power.
• 2.  
  If β_p ≪ 1 and ε_p ≳ ε*, then parallel proton heating is negligible, and AW/KAW turbulence leads to a combination of electron heating and perpendicular proton heating.
• 3.  
  If β_p ∼ 1 and ε_p ≪ ε*, then perpendicular proton heating is negligible, and AW/KAW turbulence results in a combination of electron heating and parallel proton heating.
• 4.  
  If β_p ∼ 1 and ε_p ≳ ε*, then perpendicular proton heating, parallel proton heating, and electron heating each receives an appreciable fraction of the cascade power.

To investigate the dependence of ε_i (the value of ε for thermal ions) on heliocentric distance r for different ion species in the fast solar wind, we adopt a simple analytic model for the radial profiles of the solar-wind proton density, ion temperatures, and magnetic field strength. We then apply the analytical model of Chandran & Hollweg (2009), which describes the radial dependence of the rms amplitudes of AWs at the outer scale L₀ of the turbulence, taking the velocity power spectrum P^(v)_k to be $\propto k_\perp ^{-c_3}$ for L⁻¹₀ < k_⊥ < ρ⁻¹_i. The resulting values of ε_i for protons, alpha particles, and minor ions depend strongly on c₃. However, for a fixed value of c₃, ε_i is relatively insensitive to r for 2 R_☉ < r < 1 AU. At equal temperatures, alpha particles and minor ions have larger values of ε_i than do protons, which may help to explain why minor ions are much hotter than protons in coronal holes. However, we are not yet able to determine with precision the perpendicular heating rates of different ion species as a function of r because of the uncertainties in the values of c₂ and c₃ in the solar wind, and because of the large sensitivity of the heating rates to these quantities. Future investigations into the value of c₂ for non-random-phase AW/KAW turbulence and the value of c₃ as a function of r will be particularly important for developing a more complete picture of stochastic ion heating in the solar wind.

We thank Greg Howes for providing us with the numerical value of the constant C_K that appears in Equation (30). This work was supported in part by DOE grant DE-FG02-07-ER46372, by NSF grant AGS-0851005, by NSF-DOE grants AST-0613622 and AGS-1003451, and by NASA grants NNX07AP65G and NNX08AH52G. E.Q. was supported in part by NSF-DOE grant PHY-0812811; NSF grant ATM-0752503; the David and Lucille Packard Foundation; and the Miller Institute for Basic Research in Science, University of California Berkeley.

Note added in proof. A pioneering study of stochastic heating by drift Alfvén waves was recently brought to our attention (McChesney et al. 1987). These authors studied stochastic heating numerically and in laboratory experiments, and discovered a threshold wave amplitude for stochastic heating that is equivalent to the condition ε ≃ 1 when applied to gyro-scale fluctuations.

In this Appendix, we consider the interaction between ions and low-frequency, two-dimensional (k_∥ = 0) electrostatic fluctuations with k_⊥ρ ∼ 1. We assume that ε ≪ 1, neglect magnetic-field fluctuations, and show that the leading-order non-vanishing terms in dH/dt are unable to cause secular perpendicular ion heating. We set ${\bf B} = B_0 \hat{\bf {z}}$, where B₀ is a constant. The time derivative of the ion's guiding-center position, defined in Equation (10), is then given by

$$\begin{equation} \frac{d{\bf R}}{dt} = v_z \hat{\bf {z}} + \frac{c {\bf E} \times \hat{z}}{B_0}. \end{equation} \tag{ A1 }$$

Since ε ≪ 1, the particle's orbit in the xy plane during a single gyroperiod is approximately a circle of radius ρ = v_⊥/Ω. We assume that Φ varies slowly in time, on a timescale of ∼ε⁻¹Ω⁻¹, with ∂Φ/∂z = 0. We introduce two related forms of "gyro-averages." First, if h is some physical property of a particle, such as its energy or guiding-center velocity, then we define the gyro-average of h to be

$$\begin{equation} \langle h(t) \rangle = \frac{\Omega }{2\pi } \int _{t- \pi /\Omega } ^{t+\pi /\Omega } h(t_1) \,dt_1. \end{equation} \tag{ A2 }$$

Second, for a general function of position and time g(r, t) satisfying ∂g/∂z = 0, we define the gyro-average of g(r, t) for particles with perpendicular velocity v_⊥ and guiding center R to be given by

$$\begin{equation} \langle g({\bf r},t)\rangle _{{\bf R}, v_\perp } \equiv \frac{1}{2\pi } \int _0^{2\pi } g \mbox{\bf (}{\bf R} + {\bf s}(\theta), t \mbox{\bf)}\, d\theta, \end{equation} \tag{ A3 }$$

where ${\bf s} = \hat{\bf {x}} \rho \cos (\theta) + \hat{\bf {y}}\rho \sin (\theta)$ is the vector illustrated in Figure 6.

Zoom In Zoom Out Reset image size

To simplify the notation, we define

$$\begin{equation} \overline{ g}({\bf R}, t) \equiv \langle g({\bf r},t)\rangle _{{\bf R}, v_\perp }, \end{equation} \tag{ A4 }$$

where the dependence of $\overline{ g}$ on v_⊥ is not explicitly written. If g varies slowly in time at a fixed point in space (e.g., on the timescale ε⁻¹Ω⁻¹), then $\overline{ g}({\bf R}, t)$ is (to leading order in ε) equivalent to a time average over one cyclotron period of g(r, t) evaluated at the position r(t) of a particle with guiding center R:

$$\begin{equation} \overline{ g}({\bf R},t) = \frac{\Omega }{2\pi }\int _{t-\pi /\Omega }^{t+\pi /\Omega } g\mbox{\bf (}{\bf r}(t_1), t_1\mbox{\bf)} \,dt_1. \end{equation} \tag{ A5 }$$

Thus, if we take the gyro-average of the "particle property" dR/dt in Equation (A1) using Equation (A2), we find that

$$\begin{equation} \left\langle \frac{d{\bf R}}{dt} \right\rangle = v_{\rm z} \hat{\bf {z}} + \frac{c}{B_0} \langle {\bf E} ({\bf r}, t)\rangle _{{\bf R}, v_\perp }\times \hat{\bf {z}}. \end{equation} \tag{ A6 }$$

We consider electrostatic fluctuations with ∂A/∂t = 0, and thus, E = −∇Φ. Omitting the explicit time dependence of Φ and $\overline{ \Phi }$ to simplify the notation, we can write the gyro-average of ∂Φ/∂x as

$$\begin{eqnarray} &&\left\langle \frac{\partial \Phi }{\partial x} \right\rangle _{{\bf R},v_\perp } = \lim _{\delta \rightarrow 0 } \frac{\overline{\Phi }({\bf R} + \hat{\bf {x}}\delta) - \overline{ \Phi }({\bf R})}{\delta } = \frac{\partial \overline{ \Phi }}{\partial x^\prime }, \end{eqnarray} \tag{ A7 }$$

where $\hat{\bf {x}}$ is a unit vector in the x-direction, and ∂/∂x' denotes a partial derivative with respect to the x-component of the guiding-center position R. Equation (A6) can thus be rewritten as

$$\begin{equation} \left\langle \frac{d{\bf R}}{dt} \right\rangle = v_{\rm z} \hat{\bf {z}} - \frac{c}{B_0} \nabla ^\prime \overline{ \Phi }\times \hat{\bf {z}}, \end{equation} \tag{ A8 }$$

where ∇' indicates a gradient with respect to the coordinates of the guiding-center position R. Since we have assumed ∂/∂z = 0, Equation (A8) implies that $\langle d{\bf R}/dt \rangle \cdot \nabla ^\prime \overline{ \Phi } = 0$.

We now integrate Equation (9) for an integral number of cyclotron periods, from t_a to t_b = t_a + Nδt, where δt = 2π/Ω and N ⩾ ε⁻¹. We define t₀ = t_a + δt/2 and t_j = t_j−1 + δt for any integer j. Since we have assumed that ∂A/∂t = 0, the integral of Equation (9) can be written as

$$\begin{equation} H(t_{\rm b}) - H(t_{\rm a}) = q\sum _{j=0}^{N-1}\int _{t_j-\delta t/2}^{t_j+\delta t/2} \frac{\partial \Phi }{\partial t}\,dt. \end{equation} \tag{ A9 }$$

In analogy to Equation (A7), it is straightforward to show that $\overline{ \partial \Phi /\partial t} = (\partial /\partial t) \overline{\Phi }$. We can thus rewrite Equation (A9) as

$$\begin{equation} H(t_{\rm b}) - H(t_{\rm a}) = q\sum _{j=0}^{N-1} \,\frac{\partial \overline{\Phi } }{\partial t} \mbox{\bf (}{\bf R}(t_j),t_j\mbox{\bf)}\, \delta t. \end{equation} \tag{ A10 }$$

The timescale on which $\overline{\Phi } \mbox{\bf (}{\bf R}(t),t\mbox{\bf)}$ changes by a factor of order unity is ε⁻¹δt. This is because Φ changes slowly in time at a fixed point in space, k_⊥ρ ∼ 1, and dR/dt ∼ εv_⊥. As a result, $\partial \overline{ \Phi }/\partial t$ is approximately constant within each time interval of duration δt. The right-hand side of Equation (A10) is therefore a discrete approximation of the integral of $q\partial \overline{ \Phi }/\partial t$ from t_a to t_b, with a fractional error of order ε, so that

$$\begin{equation} H(t_{\rm b}) - H(t_{\rm a}) = q\int _{t_{\rm a}}^{t_{\rm b}} \frac{\partial }{\partial t}\, \overline{\Phi } \mbox{\bf (}{\bf R}(t),t\mbox{\bf)}\, dt + \cdots, \end{equation} \tag{ A11 }$$

where the ellipsis (⋅⋅⋅) represents corrections that are of higher order in ε. The right-hand side of Equation (A11) can be rewritten in terms of the total time derivative of Φ, yielding

$$\begin{eqnarray} H(t_{\rm b}) - H(t_{\rm a}) &=& q\int _{t_{\rm a}}^{t_{\rm b}} \frac{d}{d t}\, \overline{\Phi } \mbox{\bf (}{\bf R}(t),t\mbox{\bf)}\, dt \: \nonumber \\ &&-\, q \int _{t_{\rm a}}^{t_{\rm b}} \frac{d{\bf R}}{dt} \cdot \nabla ^\prime \overline{\Phi } \mbox{\bf (}{\bf R}(t),t\mbox{\bf)}\, dt + \cdots\nonumber\\ \end{eqnarray} \tag{ A12 }$$

Since $\nabla ^\prime \overline{ \Phi }$ is nearly constant during a single time interval of duration δt, the second integral on the right-hand side of Equation (A12) satisfies the relation

$$\begin{eqnarray} \int _{t_{\rm a}}^{t_{\rm b}} \frac{d{\bf R}}{dt} \cdot \nabla ^\prime \overline{\Phi } \mbox{\bf (}{\bf R}(t),t\mbox{\bf)}\, dt &=& \sum _{j=0}^{N-1}\left(\int _{t_j-\delta t/2}^{t_j+\delta t/2} \frac{d{\bf R}}{dt}\,dt \right)\nonumber\\ &&\cdot \nabla ^\prime \overline{\Phi } \mbox{\bf (}{\bf R}(t_j),t_j\mbox{\bf)} + \cdots \end{eqnarray} \tag{ A13 }$$

The integral within parentheses on the right-hand side of Equation (A13) is equivalent to 〈dR/dt〉 δt evaluated at t = t_j. From Equation (A8), $\langle d{\bf R}/dt\rangle \cdot \nabla ^\prime \overline{ \Phi } = 0$. Thus, the right-hand side of Equation (A13) and the second integral on the right-hand side of Equation (A12) vanish to leading order in ε. Equation (A12) thus becomes

$$\begin{equation} H(t_{\rm b}) - H(t_{\rm a}) = q \overline{ \Phi }\mbox{\bf (}{\bf R}(t_{\rm b}),t_{\rm b}\mbox{\bf)} - q \overline{ \Phi }\mbox{\bf (}{\bf R}(t_{\rm a}),t_{\rm a}\mbox{\bf)} + \cdots \end{equation} \tag{ A14 }$$

The absolute value of the right-hand side of Equation (A14) remains ≲qδΦ_ρ, regardless of how large the interval (t_b − t_a) becomes. Thus, to leading order in ε, there is no secular change in the particle energy H, consistent with the near-conservation of the first adiabatic invariant in the small-ε, small-ω/Ω limits.