Thermal disequilibration of ions and electrons by collisionless plasma turbulence

Significance Large-scale astrophysical processes inject energy into turbulent motions and electromagnetic fields, which carry this energy to small scales and eventually thermalize it. How this energy is partitioned between ions and electrons is important both in plasma physics and in astrophysics. Here we determine this energy partition via gyrokinetic turbulence simulations and provide a simple prescription for the ion-to-electron heating ratio. We find that turbulence promotes disequilibration of the species: When magnetic energy density is greater than the thermal energy density, electrons are preferentially heated, whereas when it is smaller, ions are. This is a relatively rare example of nature promoting an ever more out-of-equilibrium state in an environment where particle collisions are not frequent enough to equalize the temperatures of the species.

I n many astrophysical plasma systems, such as accretion disks, the intracluster medium, and the solar wind, collisions between ions and electrons are extremely infrequent compared to dynamical processes and even compared to collisions within each species. In the effective absence of interspecies collisions, it is an open question whether there is any mechanism for the system to self-organize into a state of equilibrium between the two species and, if not, what sets the ion-to-electron temperature ratio. This is clearly an interesting plasma-physics question on a fundamental level, but it is also astrophysically important for interpreting observations of plasmas from the heliosphere to the Galaxy and beyond. Historically, the posing of this question 20 y ago in the context of radiatively inefficient accretion flows and in particular of our own Galactic Center, Sagittarius A * (Sgr A * ) [in which preferential ion heating was invoked to explain low observed luminosity (1-3)], has prompted a flurry of research and porting of analytical and numerical machinery developed in the context of fusion plasmas and of fundamental turbulence theories to astrophysical problems (see, e.g., refs. 4-12, but also ref. 13 and references therein for an alternative strand of investigations). In more recent years, heating prescriptions resulting from these investigations have increasingly been in demand for global models aiming to reproduce observations quantitatively (e.g., refs. 14 and 15 and references therein).
In a nonlinear plasma system, turbulence is generally excited by large-scale free-energy sources (e.g., the Keplerian shear flow in a differentially rotating accretion disk), then transferred to ever smaller scales in the position-velocity phase space via a "turbulent cascade," and finally converted into thermal energy of plasma particles via microscale dissipation processes. This turbulent heating is not necessarily distributed evenly between ions and electrons. It may, in principle, lead to either thermal disequilibration or equilibration between ions and electrons, depending on how the ion-to-electron heating ratio changes with the ratio of their temperatures, Ti/Te. Here we determine this dependence-along with the heating ratio's dependence (which turns out to be much more important) on the other fundamental parameter characterizing the thermal state of the plasma, the ratio of the ion-thermal to magnetic energy densities, βi.
This task requires a number of assumptions, many of which are quite simplistic, but are made here to distill what we consider to be the most basic features of the problem at hand. We assume that the large-scale free-energy injection launches a cascade of perturbations that are anisotropic with respect to the direction of the ambient mean magnetic field and whose characteristic frequencies are Alfvénic-we know both from theory (6,16) and detailed measurements in the solar wind (17) that this is what inertial-range turbulence in a magnetized plasma would look like. This means that the particles' cyclotron motion can be averaged out at all spatial scales, all the way to the ion Larmor radius and below. This "gyrokinetic" (GK) approximation (4, 18) leaves out any heating mechanisms associated with cyclotron resonances (because frequencies are low) and with Significance Large-scale astrophysical processes inject energy into turbulent motions and electromagnetic fields, which carry this energy to small scales and eventually thermalize it. How this energy is partitioned between ions and electrons is important both in plasma physics and in astrophysics. Here we determine this energy partition via gyrokinetic turbulence simulations and provide a simple prescription for the ionto-electron heating ratio. We find that turbulence promotes disequilibration of the species: When magnetic energy density is greater than the thermal energy density, electrons are preferentially heated, whereas when it is smaller, ions are. This is a relatively rare example of nature promoting an ever more outof-equilibrium state in an environment where particle collisions are not frequent enough to equalize the temperatures of the species.
shocks (19) (because sonic perturbations are ordered out). The amplitude of the fluctuations is assumed to be asymptotically small relative to the mean field, and thus stochastic heating (20) and any other mechanisms relying on finite-amplitude fluctuations (21)(22)(23)(24)(25) are also absent. Furthermore, we assume that ions and electrons individually are near Maxwellian equilibria, but at different temperatures. This excludes any heating mechanisms associated with pressure anisotropies (26)(27)(28) or significant nonthermal tails in the particle distribution functions (29,30). We note that reconnection is allowed within the GK model, and so the results obtained here include any heating, ion or electron, that might occur in reconnecting sheets spontaneously formed within the turbulent dynamics. [Note, however, that the width of the inertial range that we can afford is necessarily modest. It therefore remains an open question whether reconnecting structures that emerge in collisionless plasma turbulence in extremely wide inertial ranges (31,32) are capable of altering any of the features of ion-electron energy partition reported here.] Although the GK approximation may be viewed as fairly crude [e.g., it may not always be appropriate to neglect high-frequency fluctuations at ion Larmor scales (33)], it does a relatively good job of quantitatively reproducing solar wind observations (5); see ref. 34 for a detailed discussion of the applicability of the GK model to solar wind. In any event, such a simplification is crucial for carrying out multiple kinetic turbulence simulations at reasonable computational cost.
It can be shown that in GK turbulence, Alfvénic and compressive (slow-wave-like) perturbations decouple energetically in the inertial range (6). In the solar wind, the compressive perturbations are energetically subdominant in the inertial range (17), although it is not known how generic a situation this is. [For example, turbulence in accretion flows is mostly driven by the magnetorotational instability (MRI) (35). The partition of compressive and Alfvénic fluctuations in MRI-driven turbulence is an open question.] At low βi, it can be shown rigorously that the energy carried by the compressive cascade will always end up as ion heat. Here we ignore this heating channel and focus on the Alfvénic cascade only, bearing in mind that, at low βi, our results likely represent a lower limit on ion heating [another possible source of additional ion heating of low βi is the stochastic heating (20,25)].

Numerical Approach
An Alfvénic turbulent cascade starts in the magnetohydrodynamic (MHD) inertial range, where ions and electrons move in concert. Therefore, it is not possible to determine the energy partition between species within the MHD approximation. This approximation breaks down and the two species decouple at the ion Larmor scale, k ⊥ ρi ∼ 1, where k ⊥ is the wave number perpendicular to the mean field. At this scale, a certain fraction of the cascading energy is converted into ion heat (via linear and/or nonlinear phase mixing; see below) and the rest continues on as a cascade of "kinetic Alfvén waves" (KAWs), which ultimately heats electrons (6). The transition between these two types of turbulence is well illustrated by the characteristic shape of their spectra, familiar from solar wind measurements at βi ∼ 1 (17) (see Fig. 2, Center).
Thus, the energy partition is decided around the ion Larmor scale, where the electron kinetic effects are not important (at least in the asymptotic limit of small electron-to-ion mass ratio). We may therefore determine this partition within a hybrid model in which ions are treated gyrokinetically and electrons as an isothermal fluid (6). The isothermal electron fluid equations are derived from the electron GK equation via an asymptotic expansion in the electron-to-ion mass ratio (me/mi) 1/2 . This is valid at scales above the electron Larmor radius and so covers a broad range including both the MHD and ion-kinetic (k ⊥ ρi ∼ 1) scales. In this model, there is an assumed separation of timescales between the fluctuations and the mean fields (4), which are parameterized by fixed βi and Ti/Te values over the entire course of the simulation.
Our hybrid GK code (12) [based on AstroGK (8), an Eulerian δf GK code specialized to slab geometry] substantially reduces the cost of nonlinear simulations. It has allowed us to compute the turbulent heating in a proton-electron plasma over a broad parameter range, varying βi from 0.1 to 100 and Ti/Te from 0.05 to 100. Most space and astrophysical plasmas have βi and Ti/Te within this range. Previous GK simulations of this problem (5,(9)(10)(11) were limited to a single point in the parameter space, specifically, (βi, Ti/Te) = (1, 1), because of the great numerical cost of resolving both ion and electron kinetic scales.
To model the large-scale energy injection, we use an oscillating Langevin antenna (36), which excites Alfvén waves (AWs) by driving an external parallel current. We set the driven modes to have the oscillation frequency ωa = 0.9ωA0, the decorrelation rate γa = 0.6ωA0, where ωA0 is the AW frequency at the largest scale, and wave numbers (kx /kx0, ky /ky0, kz /kz0) = (0, 1, ±1) and (1, 0, ±1), where the subscript 0 indicates the smallest wave number in the simulation. The antenna amplitude is set to drive critically balanced turbulence, i.e., to make the nonlinear cascade rate at the driving scale comparable to the linear wave frequency ωA0.
The ions have a fully conservative linearized collision operator, including pitch-angle scattering and energy diffusion (37,38). The collision frequency is chosen to be νi = 0.005ωA0. The ions are thus almost collisionless. Since the scale range covered in our simulations is limited, these "true" collisions are not sufficient to dissipate all of the energy contained in the ion entropy fluctuations, especially at small spatial scales, where the turbulent eddy-turnover rates are higher. Therefore, we use hypercollisions with a collision frequency proportional to (k ⊥ /kmax) 8 , where kmax is the wave number corresponding to the grid scale (5). While the free energy contained in the perturbed ion distribution function is dissipated by these collisional mechanisms, the physical dissipation mechanisms for the sub-Larmor-scale turbulence destined for electron heating are ordered out by the (me/mi) 1/2 expansion. Therefore, we introduce artificial hyperdissipation (hyperviscosity and hyperresistivity) proportional to (k ⊥ /kmax) 8 in the isothermal electron fluid equations to terminate the KAW cascade (see ref. 12 for details). We carefully tune the hypercollisionality and hyperdissipation coefficients to make the artificial dissipation effective only at the smallest scales.

Energy Partition
The main result of our simulations is given in Fig. 1, which shows the dependence of the ratio of the time-averaged ion and electron heating rates Qi/Qe on βi and Ti/Te. Fig. 1, Left shows that Qi/Qe increases as βi increases regardless of Ti/Te. When (βi, Ti/Te) = (1, 1), we find Qi/Qe ≈ 0.6, in good agreement with the result found in the full GK simulation studies that resolved the entire range from MHD to electron kinetic scales (10,11). We find that ions receive more energy than electrons when βi 1 while electron heating is dominant in the low-βi regime.
Low Beta. In the limit βi → 0, our results suggest Qi/Qe → 0, which is physically intuitive: In this regime, the ion thermal speed is much smaller than the Alfvén speed, so ions cannot interact with Alfvénic perturbations and so the cascade of the latter smoothly turns into a sub-Larmor KAW cascade, without any energy being diverted into ions (41). This "smooth" transition is manifest when one examines the energy spectra in this regime (Fig. 2, Left).
The scale where the ion heating occurs is apparent in Fig. 2, Bottom. For low to moderate βi, the ion heating is dominated by grid-scale hyperdissipation. This is consistent with the previous full GK simulation with βi = 1 (9-11), where the ion heating peaked at 20 k ⊥ ρi 30. In contrast, the ion heating for high βi occurs predominantly at large scales, which is revealed in this study (next paragraph).
High Beta. In the opposite limit of high βi, simulations show that Qi/Qe increases and appears to tend to a constant 30 for βi 10.  (40); the scale ρ * at which this starts, defined in the text, is also shown). Clearly, at these resolutions, a definitive determination of spectral slopes is not feasible. Bottom panels show ion heating rate vs. k ⊥ , in units of total injected power (Qtot = Q i + Qe) times ρ i . The uptick in ion heating at the smallest scales is due to ion hyperresistivity and hyperviscosity. We note that halving the box size for the β i = 100 simulation results in only a 10% change to Q i /Qe (which is smaller than the error due to finite-time averaging), suggesting that this result is independent of injection scale.
The physics behind this result are more complicated. In a highβi plasma, AWs are damped at a rate that peaks around k ⊥ ρi ∼ β −1/4 i , where it is comparable to their propagation frequency: Namely, in the limit βi 1, the complex frequency is (4, 28) where ρ * = (3/4π 1/4 √ 2)β 1/4 i ρi. At k ⊥ ρ * > 1, AWs can no longer propagate and at k ⊥ ρ * 1, damping peters out for magnetic perturbations (ω ≈ −i|k |vA/2k 2 ⊥ ρ 2 * ), but becomes increasingly strong for velocity (electric-field) perturbations (ω ≈ −i|k |vA2k 2 ⊥ ρ 2 * ). The situation resembles an overdamped oscillator, with magnetic field in the role of displacement. This means that at k ⊥ ρ * ∼ 1, the MHD Alfvénic cascade is partially damped and partially channeled into a purely magnetic cascade, as is indeed evident in Fig. 2, Right [this resembles the subviscous cascade in high-magnetic Prandtl-number MHD and, similarly to it (40), might be exhibiting a k −1 ⊥ spectrum, arising from nonlocal advection of magnetic energy by ρ * -scale motions]. The magnetic cascade extends some way into the sub-ion-Larmor range, but eventually, at k ⊥ ρi 1, it must turn into a KAW cascade. While the sorts of spectra that we find at βi 1 (Fig. 2, Left and Center) are very similar to what has been observed both in numerical simulations (5,9,10,33,42,43) and in solar wind observations (17) at βi ∼ 1, the high-βi spectra described above have not been seen before and represent an interesting type of kinetic turbulence.
Thus, there is a finite wave-number interval of strong damping around k ⊥ ρ * ∼ 1. In a "critically balanced" turbulence, |k |vA is of the same order as the cascade rate, so this damping will divert a finite fraction of total cascaded energy into ion heat (this is manifest in Fig. 3D). Exactly what fraction it will be is what our numerical study tells us. We do not have a quantitative theory that would explain why Qi/Qe should saturate at the value that we observe numerically (which, based on a resolution study, appears to be converged). Presumably, this is decided by the details of the operation of ion Landau damping in a turbulent environment [a tricky subject (44)(45)(46)] and by the efficiency with which energy can be channeled from the MHD scales into the magnetic cascade below ρ * and the KAW cascade below ρi. In the absence of a definitive theory, Qi/Qe ≈ 30 should be viewed as an "experimental" result.
Relation to Standard Model Based on Linear Damping. It is instructive to compare Qi/Qe obtained in our simulations with the simple theoretical model for the turbulent heating proposed in ref. 7, which has been used as a popular prescription in global disk models (14,15). The model is based on assuming (i) continuity of the magnetic-energy spectrum across the ion-Larmor-scale transition, (ii) linear Landau damping as the rate of free-energy dissipation leading to ion heating, and (iii) critical balance between linear propagation and nonlinear decorrelation rates. As evident in Fig. 1, Left, Inset, the model gives a broadly correct qualitative trend, but produces some noticeable quantitative discrepancies: notably, much lower ion heating at low βi and an absence of the ceiling on Qi/Qe at high βi. This is perhaps not surprising, for a number of reasons. First, the Landau damping rate is not, in general, a quantitatively good predictor of the rate at which linear phase mixing would dissipate free energy in a driven system (47). Indeed, we have found that an approximation such as EQ i (k ⊥ ) ∝ Im ω(k , k ⊥ )EB ⊥ (k ⊥ ) (with ω the linear frequency and k either directly measured or inferred from the critical-balance conjecture) did not reproduce quantitatively the heating spectra shown in Fig. 2  lence in the no-propagation region at k ⊥ ρ * as a nonlocally driven magnetic cascade, choosing rather to smooth the frequency gap between the AWs and KAWs. Third, at low βi, as we are about to see below, the ion heating is controlled by the nonlinear, rather than linear, phase mixing ["entropy cascade" (6,33,48,49)].
Temperature Disequilibration. Apart from the βi dependence, the key finding of our simulations is that Qi/Qe is mostly insensitive to Ti/Te (keeping βi constant; Fig. 1, Right). Some dependence on Ti/Te does exist when βi 1 and Ti/Te is small [for βi 1, this is the "Hall limit" of GK (6)]. This dependence is redistributive: Colder ions are heated a little more. At low βi, most of the energy still goes into electrons, but at βi ∼ 1, the effect might be of some help in restoring some parity between ions and electrons because Qi/Qe > 1 at low Ti/Te and Qi/Qe < 1 at high Ti/Te. Overall, we see that whether ions and electrons are already disequilibrated or not makes relatively little difference to the heating rates-there is no intrinsic tendency in the collisionless system to push the two species toward equilibrium with each other (except at βi ∼ 1). In fact, in the absence of ion cooling and at constant magnetic field, turbulent heating would gradually increase βi and thus push the system toward a state of dominant ion heating and hence hotter ions. Runaway increase of Ti/Te can be envisioned if Te is capped by, e.g., radiative cooling.
Fitting Formula. For a researcher who is interested in using these results in global models (as in, e.g., refs. 14 and 15), here is PHYSICS a remarkably simple fitting formula, which, without aspiring to ultrahigh precision, works quite well over the parameter range that we have investigated (Fig. 1, Right): Phase-Space Cascades One of the more fascinating developments prompted by the interest in energy partition in plasma turbulence has been the realization that, in a kinetic system, we are dealing with a freeenergy cascade through the entire phase space, with energy travelling from large to small scales in both position and velocity space (6,33,44,45,(48)(49)(50)(51)(52)(53)(54). This is inevitable because the plasma collision operator is a diffusion operator in phase space and so the only way for a kinetic system to have a finite rate of dissipation at very low collisionality is to generate small phase-space scales-just like a hydrodynamic system with low viscosity achieves finite viscous dissipation by generating large flow-velocity gradients. The study of velocity-space cascades in kinetic systems is still in its infancy-but advances in instrumentation and computing mean that the amount of available information on such cascades in both real (space) physical plasmas (52) and their numerical counterparts (33,46,54) is rapidly increasing. Let us then investigate the nature of the phase-space cascade in our ion-heating simulations.
In low-frequency (GK) turbulence, there are two routes for the velocity-space cascade: Linear phase mixing, also known as Landau damping (55), produces small scales in the distribution of the velocities parallel to the magnetic field (v ) (47,56), whereas the cascade in the perpendicular velocities (v ⊥ ) is brought about by nonlinear phase mixing, or entropy cascade, associated with particles following Larmor orbits (whose radii are ∝ v ⊥ ) sampling spatially decorrelated electromagnetic perturbations (6,48,49). The latter mechanism switches on at spatial scales for which the Larmor radius is finite, i.e., at k ⊥ ρi 1. While these velocityspace cascades are interesting in themselves as fundamental phenomena setting the structure of plasma turbulence in phase space, they also give us a handle on whether the ion heating tends to be parallel or perpendicular (this could become important if we asked, e.g., toward what kind of pressure-anisotropic states turbulence pushes the plasma).
We use the Hermite-Laguerre spectral decomposition of the gyroaveraged perturbed distribution function g = δf (57), where Hm (x ) and L (x ) are the Hermite and Laguerre polynomials. In this language, higher m and correspond to smaller scales in v and v ⊥ , respectively. Fig. 3 shows the phase-space spectra of the ion entropy [|ĝ| 2 , the contribution of the perturbed ion distribution function to the free energy (6)] for βi = 0.1 and βi = 100 cases with Ti/Te = 1. We see that the distribution of the free energy and, consequently, the nature of its cascade through phase space change with βi.
Low Beta. At low βi, linear phase mixing is suppressed ( Fig. 3C; this is because ions' thermal motion is slow compared to the phase speed of the Alfvénic perturbations), so most of the ion entropy is cascaded simultaneously to large k ⊥ ρi and by nonlinear phase mixing (Fig. 3A) before being thermalized by collisions, giving rise to (perpendicular) ion heating. The Fourier-Laguerre spectrum contains little energy at high when k ⊥ ρi < 1 (because plasma dynamics are essentially drift kinetic at these scales and there is no phase mixing in v ⊥ ), but at k ⊥ ρi > 1 it is consistent with aligning along ∼ (k ⊥ ρi) 2 . This is a manifestation of the basic relationship between the velocity and spatial scales, δv ⊥ /v thi ∼ 1/k ⊥ ρi, that is characteristic of sub-Larmor entropy cascade (6,48,49) (δv ⊥ /v thi ∼ 1/ √ follows from the trigonometric asymptotic of Laguerre polynomials at high ). Similar "diagonal" structure has previously been found in 4D electrostatic GK simulations (58) and in 6D electromagnetic hybrid-Vlasov simulations (33). Note also that for the case (βi, Ti/Te) = (1, 1), ref. 11 compared the contributions to ion heating from the v ⊥ and v parts of the collision operator and also concluded that the nonlinear phase mixing was the dominant process.
High Beta. In contrast, at high βi, most ion entropy is channeled to high m at k ⊥ ρi < 1 (Fig. 3D) by linear phase mixing, as is indeed confirmed by the characteristic m −1/2 slope of the Hermite spectrum (41,47) [Fig. 3E; at low βi, the Hermite spectrum is steeper, implying very little dissipation (44,45)]. These perturbations are then thermalized at high m by collisions. Thus, the preferential heating of ions at high βi is parallel and occurs via ordinary Landau damping. [We make this statement with some caution. The velocity resolution of our simulations is necessarily limited, so our plasma has a certain effective collisional cutoff mc (typically, mc ∼ 10). The order of limits mc → ∞ and βi → ∞ may matter to the system's ability to block linear phase mixing via the stochastic echo effect because the rate at which free energy is transferred from m to m + 1 by linear phase mixing is ∼ |k |v thi / √ m (44,45) whereas the nonlinear advection rate in a critically balanced Alfvénic turbulence is ∼ |k |vA = |k |v thi / √ βi. At the highest values of βi , our simulations have mc < βi, so the effective collisionality may interfere with the echo. If, at infinite resolution (i.e., in an even less collisional plasma than we simulate currently), the echo is restored, ion heating at βi mc may be all via the entropy cascade.]

Discussion
To discuss an example of astrophysical consequences of our findings, let us return briefly to the curious case of low-luminosity accretion flows-most famously, the supermassive black hole Sgr A * at our Galaxy's center. Two classes of theory have been advanced to explain the observed low-luminosity, each corresponding to a distinct physical scenario: The first scenario has Qi/Qe 1 and so most of the thermal energy is deposited into nonradiating ions, which are swallowed by the black hole (1-3); the second scenario has Qi/Qe ∼ 1 but the accretion rate is very small, with most of the plasma being carried away by outflows (59). Determining which of these is closer to the truth is tantamount to identifying the fate of the accreting matter. The low-accretion rate scenario has gradually become more widely accepted (26,60,61), whereas early studies used the high-Qi/Qe scenario (2,62). The value Qi/Qe 30 that we have found for moderately high values of βi is about 10 times larger than the value used today. However, even with this value, the accretion rate must be much smaller than the Bondi rate (figure 1 in ref. 61), given the observational fact that the outflow is present (60,63). Within this scenario, the relative amount of electron heating in the low-βi, central region of the disk turns out to be crucial to enable a detectable jet: Ref. 15 found a radiating jet in global simulations using the linear prescription with very low ion heating (7) and no visible jet with a more equitable heating model (13). Our heating prescription is perhaps closer to ref. 7 in that regard, but not as extreme-it would be interesting to see what effect this has on global models of accreting systems.
On a broader and perhaps more fundamental level, we have shown that turbulence is capable of pushing weakly collisional plasma systems away from interspecies thermal equilibriumdepending on whether βi is high or low, it favors preferential thermalization of turbulent energy into ions or electrons, respectively (although at βi ∼ 1, there is some tendency to restoration of species equality). This is a relatively rare example of turbulence failing to promote Le Chatelier's principle and instead causing a disequilibration of a collisionless system.