Stochastic proton heating by kinetic-Alfv\'en-wave turbulence in moderately high-$\beta$ plasmas

Stochastic heating refers to an increase in the average magnetic moment of charged particles interacting with electromagnetic fluctuations whose frequencies are smaller than the particles' cyclotron frequencies. This type of heating arises when the amplitude of the gyroscale fluctuations exceeds a certain threshold, causing particle orbits in the plane perpendicular to the magnetic field to become stochastic rather than nearly periodic. We consider the stochastic heating of protons by Alfv\'en-wave (AW) and kinetic-Alfv\'en-wave (KAW) turbulence, which may make an important contribution to the heating of the solar wind. Using phenomenological arguments, we derive the stochastic-proton-heating rate in plasmas in which $\beta_{\rm p} \sim 1-30$, where $\beta_{\rm p}$ is the ratio of the proton pressure to the magnetic pressure. (We do not consider the $\beta_{\rm p} \gtrsim 30$ regime, in which KAWs at the proton gyroscale become non-propagating.) We test our formula for the stochastic-heating rate by numerically tracking test-particle protons interacting with a spectrum of randomly phased AWs and KAWs. Previous studies have demonstrated that at $\beta_{\rm p} \lesssim 1$, particles are energized primarily by time variations in the electrostatic potential and thermal-proton gyro-orbits are stochasticized primarily by gyroscale fluctuations in the electrostatic potential. In contrast, at $\beta_{\rm p} \gtrsim 1$, particles are energized primarily by the solenoidal component of the electric field and thermal-proton gyro-orbits are stochasticized primarily by gyroscale fluctuations in the magnetic field.


Introduction
In the mid-twentieth century several authors published hydrodynamic models of the solar wind that imposed a fixed temperature at the coronal base and took thermal conduction to be the only heating mechanism (e.g., Parker 1958Parker , 1965Hartle & Sturrock 1968;Durney 1972). These models were unable to explain the high proton temperatures and fast-solar-wind speeds observed at a heliocentric distance r of 1 astronomical unit (au) for realistic values of the coronal temperature and density, indicating that the fast solar wind is heated primarily by some mechanism other than thermal conduction. Parker (1965) and Coleman (1968) proposed that Alfvén waves (AWs) and AW turbulence provide this additional heating. Support for this suggestion can be found in the many spacecraft observations of AW-like turbulence in the solar wind (see Belcher 1971;Tu & Marsch 1995;Bale et al. 2005), remote observations of AW-like fluctuations in the solar corona (see Tomczyk et al. 2007;De Pontieu et al. 2007), and the agreement between AW-driven solar-wind models and solar-wind temperature, density, and flow-speed profiles (Cranmer et al. 2007;Verdini et al. 2010;Chandran et al. 2011;van der Holst et al. 2014).
AWs oscillate at a frequency ω = k v A , where k (k ⊥ ) is the component of the wave vector k parallel (perpendicular) to the background magnetic field, B 0 , v A = B 0 / 4πn p m is the Alfvén speed, n p is the proton number density, and m is the proton mass. † In AW turbulence, interactions between counter-propagating AWs cause AW energy to cascade from larger to smaller scales. This energy cascade is anisotropic, in the sense that the small-scale AW "eddies," or wave packets, generated by the cascade vary much more rapidly perpendicular to the magnetic field than along the magnetic field (e.g., Shebalin et al. 1983;Goldreich & Sridhar 1995;Cho & Vishniac 2000;Horbury et al. 2008;Podesta 2013;Chen 2016). As a consequence, within the inertial range (scales larger than the thermal-proton gyroradius ρ th and smaller than the outer scale or driving scale), ω ≪ Ω, where Ω is the proton cyclotron frequency. At k ⊥ ρ th ∼ 1, the AW cascade transitions to a kinetic-Alfvén-wave (KAW) cascade (Schekochihin et al. 2009).
Studies of the dissipation of low-frequency (ω ≪ Ω), anisotropic, AW/KAW turbulence based on linear wave damping (e.g., Quataert 1998;Howes et al. 2008) conclude that AW/KAW turbulence leads mostly to parallel heating of the particles (i.e., heating that increases the speed of the thermal motions along B). On the other hand, perpendicular ion heating is the dominant form of heating in the near-Sun solar wind (Esser et al. 1999;Marsch 2006;Cranmer et al. 2009;Hellinger et al. 2013). This discrepancy suggests that AW/KAW turbulence in the solar wind dissipates via some nonlinear mechanism (e.g. Dmitruk et al. 2004;Markovskii et al. 2006;Lehe et al. 2009;Schekochihin et al. 2009;Chandran et al. 2010;Servidio et al. 2011;Lynn et al. 2012;Xia et al. 2013;Kawazura et al. 2018). This suggestion is supported by studies that find a correlation between ion temperatures and fluctuation amplitudes in solar-wind measurements and numerical simulations (e.g., Wu et al. 2013;Hughes et al. 2017;Grošelj et al. 2017;Vech et al. 2018).
In this paper, we consider one such nonlinear mechanism: stochastic heating. In stochastic proton heating, AW/KAW fluctuations at the proton gyroscale have sufficiently large amplitudes that they disrupt the normally smooth cyclotron motion of the protons, leading to non-conservation of the first adiabatic invariant, the magnetic moment (McChesney et al. 1987;Johnson & Cheng 2001a;Chaston et al. 2004;Fiksel et al. 2009a;Xia et al. 2013). Chandran et al. (2010) used phenomenological arguments to derive an analytical formula for the stochastic-heating rate at β p 1, where β p is the ratio of the proton pressure to the magnetic pressure (see (2.3)). In Section 2 we use phenomenological arguments to obtain an analytic formula for the proton-stochastic-heating rate in low-frequency AW/KAW turbulence when β p ∼ 1 − 30. We limit our analysis to β p 30, since KAWs become non-propagating at k ⊥ ρ th = 1 at larger β p values (see Appendix A and Hellinger & Matsumoto (2000); Kawazura et al. (2018); Kunz et al. (2018)). In Section 3 we present results from simulations of test particles interacting with a spectrum of randomly phased AWs/KAWs, which we use to test our analytic formula for the stochastic-heating rate. Throughout this paper, we focus on perpendicular proton heating rather than parallel proton heating. Stochastic heating can in principle augment the parallel proton heating that results from linear damping of † We neglect the mass density of electrons and heavy ions. AW/KAW turbulence at β p 1, but we leave a discussion of this possibility to future work.

Stochastic proton heating by AW/KAW turbulence at the proton gyroscale
A proton interacting with a uniform background magnetic field B 0 and fluctuating electric and magnetic fields δE and δB undergoes nearly periodic motion in the plane perpendicular to B 0 if δE and δB are sufficiently small or L/ρ is sufficiently large, where L is the characteristic length scale of δE and δB, ρ = v ⊥ /Ω is the proton's gyroradius, v ⊥ is the component of the proton's velocity v perpendicular to the magnetic field, Ω = qB 0 /mc is the proton gyrofrequency, m and q are the proton mass and charge, and c is the speed of light. When (1) the proton's motion in the plane perpendicular to B 0 is nearly periodic and (2) Ωτ ≫ 1, where τ is the characteristic time scale of δE and δB, the proton's magnetic moment µ = mv 2 ⊥ /2B 0 is almost exactly conserved (Kruskal 1962).
Perpendicular heating of protons (by which we mean a secular increase in the average value of µ) requires that one of the above two conditions for µ conservation be violated. For example, Alfvén/ion-cyclotron waves can cause perpendicular proton heating via a cyclotron resonance if Ωτ ∼ 1 (Hollweg & Isenberg 2002). Alternatively, low-frequency AW/KAW fluctuations can cause perpendicular proton heating if their amplitudes at k ⊥ ρ ∼ 1 are sufficiently large that the proton motion in the plane perpendicular to B 0 becomes disordered or "stochastic" (McChesney et al. 1987;Johnson & Cheng 2001b;Chaston et al. 2004;Fiksel et al. 2009b).
We focus on this second type of heating, stochastic heating, and on "thermal" protons, where w ⊥ = 2k B T ⊥p /m and w = 2k B T p /m are the perpendicular and parallel thermal speeds, T ⊥p and T p are the perpendicular and parallel proton temperatures, k B is Boltzmann's constant, and ρ th = w ⊥ /Ω is the thermal-proton gyroradius. We restrict our attention to the contribution to the stochastic-heating rate from turbulent AW/KAW fluctuations with λ ∼ ρ th k ⊥ ρ th ∼ 1, (2.2) where λ is the length scale of the fluctuations measured perpendicular to the background magnetic field, and to (2.4) As mentioned above and discussed further in Appendix A, KAWs at k ⊥ ρ th = 1 become non-propagating at significantly larger values of β p (see also Hellinger & Matsumoto 2000;Kawazura et al. 2018;Kunz et al. 2018). For simplicity, we assume that which implies that and that T e ∼ T p , where T e is the electron temperature. We also assume that where δB ρ is the rms amplitude of the magnetic fluctuations with λ ∼ ρ th , and that the fluctuations are in critical balance (Goldreich & Sridhar 1995), which implies that where l is the correlation length of the gyroscale AW/KAW fluctuations measured parallel to the background magnetic field, and δv ρ is the rms amplitude of the E × B velocity of the AW/KAW fluctuations with λ ∼ ρ th . Since the linear and nonlinear time scales are comparable in the critical-balance model, we take the ratios of the amplitudes of different fluctuating variables to be comparable to the ratios that arise for linear AW/KAWs at k ⊥ ρ th ∼ 1, which, given (2.2) and (2.3), implies that where δB ρ and δB ⊥ρ are, respectively, the rms amplitudes of the components of the fluctuating magnetic field parallel and perpendicular to B 0 (TenBarge et al. 2012). Equations (2.2) through (2.9) imply that (2.10)

Stochastic motion perpendicular to the magnetic field
To understand how gyroscale AW/KAW fluctuations modify a proton's motion, we cannot use the adiabatic approximation (Northrop 1963), which assumes λ ≫ ρ. Nevertheless, we can still define an effective guiding center (2.11) whereb = B/B. This effective guiding center is always a distance ρ from the particle's position r and is, at any given time, the location about which the particle attempts to gyrate under the influence of the Lorentz force. We find it useful to focus on R rather than r because the motion of R largely excludes the high-frequency cyclotron motion of the proton. Upon taking the time derivative of (2.11) and making use of the relations dr/dt = v and dv/dt = (q/m)(E + v × B/c), we obtain can be found by substituting the right-hand side of (2.12) into the right-hand side of (2.13), which yields (2.14) We now estimate each term on the right-hand side of (2.14). Since we are considering only gyroscale fluctuations †, we take the first term on the right-hand side of (2.14) to satisfy the relation To estimate the second and third terms on the right-hand side of (2.14), we take which is satisfied by the majority of particles. The time derivative of the field strength along the particle's trajectory is (2.17) As outlined above, our assumption of critical balance implies that λ ≪ l and ω ≪ v ⊥ /ρ ∼ w ⊥ /ρ th . The second term on the right-hand side of (2.17) is thus much larger than either the first or third terms, and dB dt ∼ w ⊥ δB ρ ρ th .
( 2.18) The second term on the right-hand side of (2.14) thus satisfies which is larger than the first term on the right-hand side of (2.14) by a factor of ∼ β 1/2 p , given (2.6) and (2.9).
For the moment, we assume that the second term on the right-hand side of (2.14) is the dominant term; we discuss the third term in more detail below. If the second term is dominant, then dR (2.20) During a single cyclotron period 2π/Ω, a proton passes through an order-unity number of uncorrelated gyroscale AW/KAW eddies, and the values of (dR/dt) ⊥ within different gyroscale eddies are uncorrelated. If (dR/dt) ⊥ is small compared to w ⊥ , then a proton undergoes nearly circular gyromotion. However, if |(dR/dt) ⊥ | is a significant fraction of w ⊥ , then a proton and its guiding center will move in an essentially unpredictable way, and the proton's orbit will become stochastic rather than quasi-periodic. Given (2.9), is a significant fraction of unity. We illustrate how the value of δ affects a proton's motion in figure 1. We compute the particle trajectories shown in this figure by numerically integrating the equations of motion for protons interacting with randomly phased AWs and KAWs. We present the details of our numerical method and more extensive numerical results in Section 3. In the numerical calculation shown in the left panel of figure 1, δ = 0.03, and the proton's motion in the plane perpendicular to B 0 is quasi-periodic. In the numerical calculation shown in the right panel of figure 1, δ = 0.15, and the proton trajectory is more disordered or random.
We now consider the third term on the right-hand side of (2.14). The instantaneous value of this term is comparable to the instantaneous value of the second term given (2.9) and (2.16), but the third term is less effective at causing guiding-center displacements over time for the following reason. Because of (2.7), the time t required for v to change by a factor of order unity is ≫ Ω −1 . If we integrate the third term on the right-hand side of (2.14) from t = 0 to t = t f , where Ω −1 ≪ t f ≪ t , we can treat v as approximately constant in (2.14), obtaining There is, however, no secular change in the value ofb at the proton's location; the magnetic-field unit vector merely undergoes small-amplitude fluctuations about the direction of the background magnetic field. Thus, over time, the guiding-center displacements caused by the third term on the right-hand side of (2.14) are largely reversible and tend to cancel out. The third term is thus less effective than the second term at making proton orbits stochastic. When the stochasticity parameter δ defined in (2.21) exceeds some threshold, the motion of a thermal proton's guiding center in the plane perpendicular to B 0 is reasonably approximated by a random walk. To estimate the time step of this random walk, we begin by defining the cyclotron average of (dR (2.23) As stated above, during a single cyclotron period a proton's motion projected onto the plane perpendicular to B 0 carries the proton through an order-unity number of uncorrelated gyroscale AW/KAW "eddies." For simplicity, we take the amplitude and direction of each vector term on the right-hand side of (2.14) to be approximately constant within any single gyroscale eddy and the values of these vector terms within different eddies to be uncorrelated. This makes v R approximately equal to the average of some order-unity number of uncorrelated vectors of comparable magnitude. The amplitude of this average is comparable to the instantaneous value of |(dR/dt) ⊥ |. Thus, given (2.6), (2.9), and ( 2.20), Because we are considering the effects of just the gyroscale AW/KAW eddies, v R decorrelates after the proton's guiding center has moved a distance ∼ ρ th in the plane perpendicular to B 0 , which takes a time Thermal protons thus undergo spatial diffusion in the plane perpendicular to B 0 with a spatial diffusion coefficient ( 2.27) The time required for a particle to wander a distance ∼ ρ th perpendicular to the background magnetic field is thus comparable to the time required for the particle to traverse the parallel dimension of a gyroscale AW/KAW eddy.

Energy diffusion and heating
The total energy of a proton is given by its Hamiltonian, where Φ is the electrostatic potential, p is the canonical momentum, and A is the vector potential. From Hamilton's equations, The second term on the right hand side of (2.29) is qv · E s , where E s = −c −1 ∂A/∂t is the solenoidal component of the electric field. Equation (46) of Hollweg (1999) gives the ratio of E s to the magnitude of the irrotational component of the electric field |∇Φ| for In their treatment of stochastic heating at β p 1, Chandran et al. (2010) neglected the second term on the right-hand side of (2.29), because this term makes a small contribution to the heating rate when β p is small. Here we focus on the effects of E s and make the approximation that dH dt ∼ qv · E s . (2.31) We show in Appendix B that the irrotational part of the electric field contributes less to the heating rate than does the solenoidal part when β p 1. As a proton undergoes spatial diffusion in the plane perpendicular to the background magnetic field, the electromagnetic field at its location resulting from gyroscale AW/KAW fluctuations decorrelates on the time scale ∆t given in (2.27). Within each time interval of length ∼ ∆t, the proton energy changes by an amount δH (which can be positive or negative), and the values of δH are uncorrelated within successive time intervals of length ∆t. As a consequence, the proton undergoes energy diffusion.
To estimate the rms value of δH, which we denote ∆H, we adopt a simple model of a proton's motion, in which the proton's complicated trajectory is replaced by a repeating two-step process. In the first step, the proton undergoes circular cyclotron motion in the plane perpendicular to B 0 for a time ∆t. In the second step, the proton is instantly translated a distance ρ th in some random direction perpendicular to B 0 . † In this simple model, a proton undergoes N ∼ Ω∆t ∼ (v ⊥ /ρ th ) × (l/v ) ∼ l/ρ th ≫ 1 circular gyrations in the plane perpendicular to B 0 during a time ∆t. Integrating (2.31) for a time ∆t, we obtain where r(t) is the proton's position at time t. Since ∆t ∼ β −1/2 p ρ th /δv ρ , when β p 1, the time ∆t for a particle to diffuse across one set of gyroscale eddies is shorter than or comparable to the linear or nonlinear time scale ρ th /δv ρ of those eddies. We thus approximate the right-hand side of (2.32) by setting E s (r(t), t) = E s (r(t), 0) and rewrite (2.32) as where the line integral is along the proton's path during one complete circular gyration in the plane perpendicular to B 0 , the surface integral is over the circular surface S of radius ρ th enclosed by the gyration, and we have used Faraday's Law ∇ × E = (−1/c)∂B/∂t. The surface S is perpendicular to B 0 , and dS is anti-parallel to B 0 (antiparallel rather than parallel since q > 0). The rms value of δH thus satisfies the orderof-magnitude relation where is the nonlinear frequency of the gyroscale fluctuations. Upon setting q/c = Ωm/B, N = Ω∆t, and ρ 2 th = w 2 ⊥ /Ω 2 in (2.34), we obtain (2.36) Although we are in the process of estimating the rate at which µ changes over long times, our estimate of ∆H is comparable to the value that would follow from µ † For simplicity, our model of proton motion neglects motion parallel to B. This approximation is to some extent justified by (2.27), which implies that a proton is unable to escape from an eddy of length l by motion along the magnetic field in a time shorter than ∆t. However, we return to the issue of parallel motion in Section 2.3. conservation: ∆H ∼ µ∆B ∼ (mw 2 ⊥ /B)ω eff δB ρ ∆t, where ∆B ∼ ω eff δB ρ ∆t is the rms amplitude of the change in the magnetic flux through the proton's Larmor orbit, divided by πρ 2 , during the time ∆t in which the proton is (in our simple two-step model of proton motion) undergoing continuous, circular, cyclotron motion. This correspondence highlights an alternative interpretation of the stochastic-heating process at β p 1. In the guiding-center approximation, when v 2 ⊥ increases by some factor because of E s , the field strength at the particle's guiding center increases by approximately the same factor, essentially because of Faraday's law. This proportionality underlies µ conservation. In stochastic heating, the same proportionality is approximately satisfied during a single time interval ∆t, but the proton is then stochastically transported to a neighboring set of gyroscale eddies, in which the field strength is not correlated with the field strength at the proton's original location. The proton thus "forgets" about what happened to the field strength at its original location and gets to keep the energy that it gained without "paying the price" of residing in a higher-field-strength location. In this way, spatial diffusion perpendicular to B breaks the connection between changes to v 2 ⊥ and changes to B that arises in the ρ/λ → 0 limit. In our simple model, the energy gained by a proton is in the form of perpendicular kinetic energy, because we neglect the parallel motion of protons. (We do not preclude the possibility of parallel stochastic heating, but we do not consider it further here.) The perpendicularkinetic-energy diffusion coefficient D K is thus ∼ ∆H 2 /∆t, or where we have used (2.27) to estimate ∆t and (2.9) to set δB ρ ∼ δB ρ . A single proton undergoing a random walk in energy can gain or lose energy with equal probability during a time ∆t. However, if a large number of thermal protons (e.g., with an initially Maxwellian distribution) undergo energy diffusion, then on average more protons will gain energy than lose energy, leading to proton heating. The heating time scale τ h is the characteristic time for the perpendicular kinetic energy of a thermal proton to double, τ h ∼ (mw 2 ⊥ ) 2 /D K , and the perpendicular-heating rate per unit mass is (2.39) To account for the uncertainties introduced by our numerous order-of-magnitude estimates, we multiply the right-hand side of (2.39) by an as-yet-unknown dimensionless constant σ 1 . As δv ρ → 0, dµ/dt decreases faster than any power of δv ρ (Kruskal 1962). To account for this "exponential" µ conservation in the small-δv ρ limit, we follow Chandran et al. (2010) by multiplying the right-hand side of (2.39) by the factor exp(−σ 2 /δ), where σ 2 is another as-yet-unknown dimensionless constant, and δ is defined in (2.21). For comparison, the stochastic-heating rate per unit mass found by Chandran et al.
(2010) when β p 1 is and the dimensionless constants c 1 and c 2 serve the same purpose as those in (2.40). As discussed by Chandran et al. (2010) for the case of c 1 and c 2 , we expect the constants σ 1 and σ 2 to depend on the nature of the fluctuations. For example, at fixed δv ρ , we expect stronger heating rates (i.e., larger σ 1 and/or smaller σ 2 ) from intermittent turbulence than from randomly phased waves (Chandran et al. 2010;Xia et al. 2013;Mallet et al. 2018), because, in intermittent turbulence, most of the heating takes place near coherent structures in which the fluctuations are unusually strong and in which the proton orbits are more stochastic than on average.

Orbit stochasticity from parallel motion
In Section 2.1, we focused on proton motion perpendicular to B. However, motion along the magnetic field can also produce stochastic motion in the plane perpendicular to B 0 (see, e.g., Hauff et al. 2010). In particular, the perpendicular magnetic fluctuations at the scale of a proton's gyroradius perturb the direction ofb. These perturbations, when fed into the first term on the right-hand side of (2.12), v b , cause the proton's guiding center R to acquire a velocity perpendicular to B 0 of where δB ⊥ρ is the component of δB (from gyroscale fluctuations) perpendicular to B 0 at the proton's location. The value of u ⊥ varies in an incoherent manner in time, with a correlation time ∼ Ω −1 . If u ⊥ is a significant fraction of v ⊥ , then u ⊥ will cause a proton's orbit in the plane perpendicular to B 0 to become stochastic. This leads to an alternative high-β p stochasticity parameter, Asδ increases towards unity, proton orbits become stochastic. For thermal protons with v ⊥ ∼ v and ρ ∼ ρ th ,δ is equivalent to δ in (2.21), which was based upon the parallel magnetic-field fluctuation δB ρ (even though we set δB ρ ∼ δB ρ in (2.21))). The contribution of parallel motion to orbit stochasticity thus does not change our conclusions about the rate at which thermal protons are heated stochastically. However, the contribution of parallel motion to orbit stochasticity should be taken into account when considering the ability of stochastic heating to produce superthermal tails, because in AW turbulence the perpendicular (parallel) magnetic fluctuation at perpendicular scale λ, denoted δB ⊥λ (δB λ ), is an increasing (decreasing) function of λ when λ is in the inertial range. Orbit stochasticity through the interaction between parallel motion and δB ⊥ρ could thus contribute to the development of superthermal tails when β p 1. An investigation of superthermal tails, however, lies beyond the scope of this paper.

Numerical Test-Particle Calculations
To test the phenomenological theory developed in Section 2, we numerically track test-particle protons interacting with a spectrum of low-frequency randomly phased AWs and KAWs. The initial particle positions are random and uniformly distributed within a cubical region of volume (100d p ) 3 , where d p = v A /Ω is the proton inertial length. The initial velocity distribution is an isotropic Maxwellian with proton temperature T p . To trace each particle, we solve the equations of motion, using the Boris method (Boris 1970) with a time step of 0.01Ω −1 .

Randomly phased waves
The code used to implement the AW/KAW spectrum is similar to the code used by Chandran et al. (2010). The magnetic field is B = B 0ẑ + δB, where B 0 is constant. We take E and δB to be the sum of the electric and magnetic fields of waves at each of 81 different wave vectors, with two waves of equal amplitude at each wave vector, one with ω/k z < 0 and the other with ω/k z > 0. † The initial phase of each wave is randomly chosen.
The 81 wave vectors correspond to nine evenly spaced values of the azimuthal angle in k space (in cylindrical coordinates aligned with B 0 ) at each of nine specific values of k ⊥i : i ∈ [0, . . . , 8]. The values of k ⊥i are evenly spaced in ln(k ⊥ )-space, with ln(k ⊥i ρ th ) = −4/3 + i/3. The middle three cells, in which i = 3, 4, and 5, have a combined width of unity in ln(k ⊥ )-space, centred at precisely k ⊥ ρ th = 1. We computationally evaluate δv ρ and δB ρ via the rms values of the E × B velocity and δB that result from the waves in just these middle three cells.
There is one value of k ≡ |k z | at each k ⊥i , denoted k i . We determine k 4 by setting the linear frequency at k 4 equal to k ⊥4 δv ρ . At other values of k ⊥ , we set : 4 < i 8. (3.2) The exponents 2/3 and 1/3 in (3.2) are chosen to match the scalings in the criticalbalance models of Goldreich & Sridhar (1995) and Cho & Lazarian (2004), respectively. We take the individual wave magnetic-field amplitudes to be proportional to k −1/3 ⊥ and k −2/3 ⊥ for k ⊥ ρ th < 1 and k ⊥ ρ th > 1, respectively, in order to match the same two criticalbalance models. (All the waves at the same value of k ⊥i have the same amplitudes.) We determine the wave frequency and relative amplitudes of the different components of the fluctuating electric and magnetic fields using Hollweg's (1999) two-fluid analysis of linear KAWs, setting where T e and T p are the (isotropic) electron and proton temperatures. We do not expect the particular choices in (3.3) to have a large effect on our results, but choose those values to facilitate a direct comparison to the previous numerical results of Chandran et al. (2010). Since we take T ⊥p = T p , we set (3.4)

A Note on the Electric Field
Following Lehe et al. (2009), we correct the electric field because the magnetic field (including its fluctuations, i.e., B = B 0ẑ + δB) in the simulation is not orientated along the z-axis. The simulation, however, equates the parallel and perpendicular components of the electric field to the parallel and perpendicular components of the wave electric field that would arise if the magnetic field were aligned exactly on the z-axis. The result is a numerical addition of perpendicular electric field terms to the parallel electric field, which, in turn, causes non-physical parallel heating. This may be seen in figure 2 of Chandran et al. (2010). To fix this, we replace the sum of the individual wave electric fields described in Section 3.1, which we denote E wave , with the modified electric field E = E wave +b(ẑ · E wave −b · E wave ).

Perpendicular Heating
We perform numerical test-particle calculations at five different values of β p , in particular β p = {0.006, 0.01, 0.1, 1, 10}. For each β p value, we carry out a test-particle calculations for five different values of δ (or, equivalently, five different values of ε). Each calculation returns the value of v 2 ⊥ and v 2 as functions of time. We show two examples in figure 2. The slope of the best-fit line for each v 2 ⊥ (t) curve determines the perpendicular-heating rate per unit mass Q ⊥ = (1/2)(d/dt) v 2 ⊥ , where · · · indicates an average over the 10 4 or 10 5 particles in the simulation. (We use more particles in simulations with smaller ε and δ because the heating rates are smaller in these simulations, and the extra particles increase the signal-to-noise ratio.) We fit the v 2 ⊥ (t) curves during the time interval (t i , t f ), where t i = 20π/Ω and t f is the smaller of the following two values: 10 4 Ω −1 and the time required for v 2 ⊥ to increase by ≃ 30%. We do not include the first ten cyclotron periods when calculating Q ⊥ , because it takes the particles a few cyclotron periods to adjust to the presence of the waves, during which time there is typically strong transient heating. (As figure 2 shows, the test particles undergo parallel heating as well as perpendicular heating, as was found previously by Xia et al. (2013) in simulations of test particles interacting with reduced magnetohydrodynamic turbulence at β p = 1.) The perpendicular-heating rates in our test-particle calculations are shown in figure 3. The agreement between our numerical results and (2.40) suggests that the approximations used to derive this equation are reasonable.
The lower-left panel of figure 3 shows that at β p < 1, Q ⊥ /(Ωw 2 ) is a function of ε alone, consistent with the fact that (2.41) can be rewritten in the form The top-right panel of figure 3 shows that at β p 1, Q ⊥ /(Ωv 2 A ) is a function of δ alone, consistent with the fact that (2.40) can be rewritten as We note that in our model of randomly phased KAWs (Chandran et al. 2010), and thus As a consequence, if we adopt the best-fit values of σ 1 , σ 2 , c 1 , and c 2 , then the value of Q ⊥ at β p = 1 in (2.40), which matches our test-particle calculations quite well, exceeds the value that would follow from (2.41) at β p = 1. A similar phenomenon was found by Xia et al. (2013) in numerical simulations of test particles interacting with strong reduced magnetohydrodynamic turbulence.
To obtain a fitting formula that can be used to model stochastic heating at large β p , small β p , and β p ≃ 1, we use (3.4) and (3.10) to rewrite the low-β p heating rate in (3.7) in terms of δ and v A . We then add the low-β p heating rate to the high-β p heating rate in (3.8), obtaining (3.11) The first term on the right-hand side dominates at β p 1 in part because σ 1 ≃ 6.5c 1 . Figure 3: Numerical results for the perpendicular-heating rate per unit mass, Q ⊥ , for protons interacting with randomly phased AWs and KAWs. Top-Left: β p < 1, and Q ⊥ normalized by Ωv 2 A . Top right: β p 1 and Q ⊥ normalized by Ωv 2 A . Bottom left: β p < 1 and Q ⊥ normalized by Ωw 2 , where w is the proton thermal speed defined in (3.4); the numerical results for all three β p values (0.006, 0.01, and 0.1) are within the bars shown. Bottom right: β p 1 and Q ⊥ normalized by Ωw 2 . In the left two panels, the solid lines are plots of (2.41) for the best-fit values of c 1 and c 2 given in (3.5). In the right two panels, the solid lines are plots of (2.40) for the best-fit values of σ 1 and σ 2 in (3.6).
The second term on the right-hand side dominates at β p ≪ 1. Figure 4 shows that (3.11) is consistent with our numerical results. This figure also illustrates how, at fixed δB ρ /B 0 , the stochastic heating rate increases as β p decreases.
As mentioned above, stochastic heating becomes more effective as the fluctuations become more intermittent (Xia et al. 2013;Mallet et al. 2018). The randomly phased waves Figure 4: The data points reproduce the numerical results from figure 3 for β p = 0.01, 0.1, 1.0, and 10. The dotted, long-dashed, solid, and short-dashed lines plot (3.11) for, respectively, β p = 0.01, β p = 0.1, β p = 1, and β p = 10. The solid and short-dashed lines are difficult to distinguish because they are nearly on top of each other. in our test-particle simulations are not intermittent, but gyroscale fluctuations in space and astrophysical plasmas generally are (see, e.g., Mangeney et al. 2001;Carbone et al. 2004;Salem et al. 2009;Chandran et al. 2015;Mallet et al. 2015). Further work is needed to determine how the best-fit constants in (3.5) and (3.6) depend upon the degree of intermittency at the proton gyroradius scale. Until this dependency is determined, some caution should be exercised when applying (3.11) to space and astrophysical plasmas. For reference,  found that lowering c 2 to ≃ 0.2 led the heating rate in (2.41) to be consistent with the proton heating rate and fluctuation amplitudes inferred from measurements of the fast solar wind from the Helios spacecraft at r = 0.3 au. However, if c 2 = 0.33, then the heating rate in (2.41) is too weak to explain the proton heating seen in the Helios measurements.

Summary
In this paper we use phenomenological arguments to derive an analytic formula for the rate at which thermal protons are stochastically heated by AW/KAW turbulence at k ⊥ ρ th ∼ 1. We focus on β p ∼ 1 − 30. Smaller values of β p were considered by Chandran et al. (2010). At larger values of β p , KAWs at k ⊥ ρ th ∼ 1 become nonpropagating, and some of the scalings we have assumed do not apply. At β p ∼ 1 − 30, the motion of a proton's effective guiding center is dominated by the interaction between the proton and gyroscale fluctuations in the magnetic field, whose amplitude is denoted δB ρ . As δB ρ /B 0 increases from infinitesimal values towards unity, the proton motion in the plane perpendicular to B 0 becomes random (stochastic), leading to spatial diffusion, and this spatial diffusion breaks the strong correlation between changes in a proton's perpendicular kinetic energy and the magnetic-field strength at the proton's location that normally gives rise to magnetic-moment conservation. The interaction between the proton and the electric field then becomes a Markov process that causes the proton to diffuse in energy. This energy diffusion leads to heating. At β p ∼ 1−30, it is the solenoidal component of the electric field that dominates the heating.
The analytic formula that we derive for the stochastic heating rate Q ⊥ contains two dimensionless constants, σ 1 and σ 2 , whose values depend upon the nature of the AW/KAW fluctuations that the proton interacts with (e.g., randomly phased waves or intermittent turbulence). We numerically track test particles interacting with randomly phased AWs and KAWs and find that our analytic formula for Q ⊥ agrees well with the heating rate of these test particles for the choices σ 1 = 5.0 and σ 2 = 0.21. We note that previous work has shown that for fixed rms amplitudes of the gyroscale fluctuations, stochastic heating is more effective when protons interact with intermittent turbulence than when protons interact with randomly phased waves (Chandran et al. 2010;Xia et al. 2013;Mallet et al. 2018). The reason for this is that in intermittent turbulence, most of the heating occurs near coherent structures, in which the fluctuation amplitudes are larger than average and in which the particle orbits are more stochastic than on average.
Our work leaves a number of interesting questions unanswered. Two such questions are how the energy-diffusion coefficient depends on energy at β p ∼ 1 − 30 and how the proton distribution function evolves in the presence of stochastic heating. (For a discussion of the low-β p case, see Klein & Chandran (2016).) We also have not addressed the question of how stochastic heating changes as β p is increased to values 30 and KAWs at k ⊥ ρ th ∼ 1 become non-propagating, or how the stochastic heating rate for minor ions depends upon minor-ion mass, charge, and average flow speed along B 0 in the proton frame. (For a discussion of the low-β p case, see Chandran et al. (2013).) Previous studies have compared observationally inferred heating rates in the solar wind with the low-β p stochastic-heating rate in (2.41) derived by Chandran et al. (2010), finding quantitative agreement at r = 0.3 au assuming c 2 ≃ 0.2  and qualitative agreement at r = 1 au (Vech et al. 2017). However, it is not yet clear whether the stochastic heating rate in (2.40) agrees with solar-wind measurements in the large-β p regime. In addition, stochastic heating at β p 1 could trigger temperatureanisotropy instabilities, which could in turn modify the rate(s) of perpendicular (and parallel) proton heating. Future investigations of these questions will be important for determining more accurately the role of stochastic heating in space and astrophysical plasmas.
We thank L. Arzamasskiy, Y. Kawazura, M. Kunz, E. Quataert, and A. Schekochihin for valuable discussions. This work was supported in part by NASA grants NNX15AI80, NNX16AG81G, NNS16AM23G, NNX17AI18G, and NNN06AA01C, and NSF grant PHY-1500041. D. Verscharen acknowledges the support of STFC Ernest Rutherford Fellowship ST/P003826/1. values, KAWs are non-propagating throughout an interval of k ⊥ ρ th values centered on unity that broadens to both larger and smaller values as β p increases (Kawazura et al. 2018).
Appendix B. Stochastic heating by the electrostatic potential at β p 1 In Section 2 we considered the rms change to a thermal proton's energy ∆H resulting from the solenoidal component of the electric field E s during the particle residence time ∆t within one set of gyroscale eddies. We also evaluated the contribution of E s to the stochastic heating rate Q ⊥ . Here, we show that the contribution to Q ⊥ from E s is larger than the contribution from the electrostatic potential Φ when β p 1.
We assume that the rms amplitude of the potential part of the electric field at k ⊥ ρ th is comparable to the rms amplitude of the total gyroscale electric-field fluctuation, δE ρ , which in turn is ∼ δv ρ B 0 /c. As discussed by Chandran et al. (2010), the contribution of the time-varying electrostatic potential to ∆H is ∆H potential ∼ qω eff ∆Φ ρ ∆t, (B 1)