TURBULENCE AND PROTON–ELECTRON HEATING IN KINETIC PLASMA

The nature of the nonlinear energy cascade and its effects on thermal properties of turbulent plasma are central problems in understanding the heating of the solar corona, the origin of the solar wind, and numerous other problems in space physics (Matthaeus & Velli 2011). The same fundamental physics issues are broadly relevant in turbulent astrophysical plasmas (Lazarian et al. 2012) and enter in diverse circumstances including molecular clouds, where turbulent heating influences star formation (Mac Low 1999) and cooling flows where turbulent fluctuations regulate heat transport (Chandran & Cowley 1998; Banerjee & Sharma 2014). A standard turbulence scenario inherited from hydrodynamics and magnetohydrodynamics (MHD) at high Reynolds numbers emphasizes a power-law inertial range spectrum that describes the steady distribution of energy across scales, a matter of great importance, e.g., in the resonant scattering of energetic particles (Jokipii 1966). However, that the inertial range acts as a near-lossless conduit that connects energy-containing large scales to the dissipative small scales is of equal importance. This energy transfer across scales in effect "completes" the circuit that allows a system far from equilibrium to relax toward thermal equilibrium. In sufficiently large hydrodynamic and MHD systems, it is widely recognized that the energy-containing eddies control the rate of transfer into the inertial range and therefore control the dissipation rate, on average. For a kinetic plasma, the theoretical understanding is less developed and perhaps more controversial, but there is also accumulating evidence that large-scale reservoirs of energy exert and regulate the rate of relaxation, fueling the the response of the plasma physics at much smaller scales, leading eventually to dissipation (Karimabadi et al. 2013). The present study further elucidates the kinetic plasma cascade process, focusing on the same two energy fluxes: the flux into the inertial range and the flux into dissipation channels. We proceed based on empirical evidence from simulations, extending an earlier preliminary report (Wu et al. 2013). Progress on these issues will lead to a better understanding of heating and other crucial effects of turbulence in space and astrophysical applications.

The standard von Karman approach to computing turbulent decay rates (de Kármán & Howarth 1938) may be adapted to MHD (Hossain et al. 1995; Wan et al. 2012). A similarity-law decay rate, written in terms of the Elsässer amplitudes ${Z}_{+}$ and ${Z}_{-}$ and their respective similarity length scales ${\lambda }_{+}$ and ${\lambda }_{-}$, may be written as $\epsilon ={C}_{+}{Z}_{+}^{3}/{\lambda }_{+}+{C}_{-}{Z}_{-}^{3}/{\lambda }_{-}$. ${C}_{+}$ and ${C}_{-}$ are von Karman constants that may depend on other parameters. A simplified version, often employed when normalized cross helicity ${\sigma }_{{\rm{c}}}=({Z}_{+}^{2}-{Z}_{-}^{2})/({Z}_{+}^{2}+{Z}_{-}^{2})$ is near zero, is $\epsilon \approx {C}_{{\rm{K}}}\tfrac{{Z}^{3}}{\lambda }$ containing a single correlation scale λ and a single von Karman constant C_K. Adopting the latter for simplicity, the decay rate at large scales is related to the dissipation and the rate of heating as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {Z}^{2}}{\partial t}={C}_{{\rm{K}}}(\ast )\displaystyle \frac{{Z}^{3}}{\lambda }={Q}_{p}+{Q}_{e}.\end{eqnarray} \tag{ 1 }$$

Here, protons and electrons are the only species present, and Q_p and Q_e are their respective heating functions, or local rates of addition of internal energy.

Note that we write the von Karman "constant" as ${C}_{{\rm{K}}}(\ast )$ to emphasize that it is expected to depend on parameters such as cross helicity, plasma beta, and the ratio of kinetic to magnetic energy as well as Reynolds number or system size. Figure 1 shows the total heating rate ${Q}_{p}+{Q}_{e}$, obtained as the average rate of temperature increase for each species, as a function of von Karman decay rate ${Z}^{3}/{\lambda }_{{\rm{c}}}$, for each run in Table 1. To compute the (time-averaged) heating rates from the simulations, we evaluated the change in the temperature over a period of five nonlinear times and divide by that time interval to give thermal energy input per nonlinear time. That is, e.g., for protons, ${Q}_{p}=({E}_{\mathrm{th},p}(t=5{\tau }_{\mathrm{nl}})-{E}_{\mathrm{th},p}(t=0))/5{\tau }_{\mathrm{nl}}$ (where τ_nl is the global nonlinear time; see below), and similarly for electrons. Within each family of runs a linear relationship is apparent, indicating a residual dependence on other parameters in the von Karman constant, which can be accommodated in the parametric dependence of ${C}_{{\rm{K}}}(\ast )$ as anticipated above. (Such parameters are constant within each family.) This result further strengthens the conclusions of Wu et al. (2013) in affirming the applicability of the von Karman average decay rate. In that earlier work the approach was somewhat different, examining the temporal dependence of the fluctuation energy decay rate normalized to its von Karman value, the full version of which is given just above Equation (1). Note that the runs other than Orszag-Tang runs employed in the present study (see Table 1) have low cross-helicity (${\sigma }_{{\rm{c}}}\sim 0$) and comparable length scales ${\lambda }_{+}\sim {\lambda }_{-}$ for the two Elsässer fields; therefore, the simplified form of given in Equation (1) is relevant.

Zoom In Zoom Out Reset image size

We now direct our attention to the question of how the proton heating responds to the strength of the cascade. More specifically, we will examine how the ratio of heating rates ${Q}_{p}/{Q}_{e}$ varies with the total cascade rate, as measured by the total heating rate ${Q}_{p}+{Q}_{e}$. Before presenting new results, we call attention to several suggestions made in the recent literature.

An important clue comes from the observational study of solar wind data (Smith et al. 2006), which found steeper sub-proton scale slopes of the magnetic energy spectrum in cases that have larger estimated cascade rates. The suggestion was that a more robust turbulence cascade at larger MHD scales leads to greater dynamical response by thermal protons. A subsequent solar wind observational study (Matthaeus et al. 2008) bolstered this conclusion by examining the influence of the cascade rate on Taylor microscale ${\lambda }_{{\rm{T}}}$, which is a measure of the root mean square curvature of the magnetic field. For smaller cascade rates, ${\lambda }_{{\rm{T}}}$ can actually be smaller than the ion inertial scale d_i, while for large cascade rates ${\lambda }_{{\rm{T}}}\gt {d}_{i}$. To interpret this, one may note that the hydrodynamic ordering of Taylor and dissipation scale η is always ${\lambda }_{{\rm{T}}}\gt \eta$. Then, if one associates the "steepening scale" of the solar wind magnetic spectrum with d_i, the suggestion emerges that the magnetic field is more greatly smoothed by processes active at scales near d_i when the cascade rate is larger.

Based on the above reasoning, the working hypothesis becomes that dissipation is enhanced near protons scales when the total dissipation rate is greater. This is basically because larger cascade rate in a Kolmogorov-style cascade increases the small scale fluctuations. In particular, the typical component magnetic increments $\delta {B}_{{\ell }}$ at scale ℓ behave as $\delta {B}_{{\ell }}\sim {\epsilon }^{1/3}{{\ell }}^{1/3}$, so it is apparent that for ion scales ${\ell }={d}_{i}$, the fluctuations will increase as ${\epsilon }^{1/3}$. This causes a greater perturbation to proton orbits, which will increase the heating of protons. It also leaves less energy in the sub-proton scale cascade, resulting in less heating of electrons. At large fluctuation amplitude we expect that the relative heating of protons should saturate at some value. We may formulate this expectation as

$$\begin{eqnarray}&&\displaystyle \frac{{Q}_{p}}{{Q}_{e}}\sim {\epsilon }^{1/3}({\rm{small}}\,\epsilon );\,{\rm{while}}\,{Q}_{p}+{Q}_{e}\sim \epsilon ,\end{eqnarray} \tag{ 2 }$$

the former of which will be subject to saturation at large amplitude, while the second relation remains accurate on average. These relations are directly testable in kinetic plasma simulation. Using the same ensemble of simulations described in Table 1, we tested this hypothesis; the results are illustrated in Figure 2. There is considerable scatter in this figure, in part due to the fact that these runs are from different families of initial conditions. However, it is apparent that there is a general trend for increased relative relative heating of protons compared to heating of electrons in runs that have larger total heating. This result is consistent with the findings of Wu et al. (2013).

Zoom In Zoom Out Reset image size

To have an effective transfer of energy between collective fluid scales and random thermal motions of particles in a magnetized plasma, one must either energize particles in the parallel (to magnetic field) direction or disrupt the perpendicular gyromotion. Parallel energization is most easily accomplished in current channels in which there are low levels of transverse magnetic fluctuations. This process is obviously more effective for electrons, given their smaller gyroradius, and has been documented for test particle electrons in MHD turbulence and in kinetic simulations (Dmitruk et al. 2004; Drake et al. 2006). This type of parallel energization of protons may also occur, but will tend to operate at lower energies. Disruption of perpendicular gyromotion is necessary to break the adiabatic invariance of magnetic moment (see e.g., Chandran et al. 2010). The effectiveness of this disruption depends on the degree to which the magnetic field at the scale of the particle gyro-orbit changes during a gyroperiod. In turbulence, the magnetic field moving with the material elements (here, the ExB drift velocity) is reconfigured over a timescale given by the nonlinear time ${\tau }_{\mathrm{nl}}({\ell })$ at the relevant scale ℓ. Therefore, the fractional reconfiguration in a gyration timescale ${\tau }_{{\rm{c}}}={{\rm{\Omega }}}_{{\rm{c}}}^{-1}=\tfrac{{mc}}{e\tilde{B}}$ (total rms magnetic field $\tilde{B}$ including mean and fluctuations) is given by ${\tau }_{{\rm{c}}}/{\tau }_{\mathrm{nl}}({\ell })$.

The scale-dependent nonlinear time may be evaluated from standard von Karman–Kolmogorov phenomenology. The large-scale or global nonlinear time is ${\tau }_{\mathrm{nl}}=\lambda /Z$, where Z is the mean turbulence amplitude and λ the correlation or energy-containing scale. (Again, for simplicity, we specialize to zero- or low-cross helicity and a single length scale.) Then, one estimates the turbulence amplitude at smaller inertial range scales as ℓ given by $\delta {Z}_{{\ell }}=Z{({\ell }/\lambda )}^{1/3}$, consistent with both the 1941 and 1962 Kolmogorov inertial range scalings. The nonlinear timescale at ℓ is ${\tau }_{\mathrm{nl}}{({\ell })={\ell }/\delta {Z}_{{\ell }}={\tau }_{\mathrm{nl}}({\ell }/\lambda )}^{2/3}$.

The simplest choice for ℓ is the ion inertial scale d_i, which has often been associated with the spectral break observed in the solar wind (Leamon et al. 1998) as well as current sheet thicknesses in plasma simulations (Shay et al. 1998) and in laboratory experiments (Brown et al. 2006). However, there are also reasonable arguments for selecting the thermal gyroradius ${\rho }_{i}$ as the relevant length scale for heating and we will discuss this below. For now, we adopt the hypothesis the heating of protons is sensitive to the nonlinear timescale at the ion inertial scale and in the simulation database we empirically test the hypothesis that

$$\begin{eqnarray}&&\displaystyle \frac{{Q}_{p}}{{Q}_{e}}\sim \displaystyle \frac{{\tau }_{{\rm{c}}}}{{\tau }_{\mathrm{nl}}({d}_{i})}=\displaystyle \frac{{\tau }_{{\rm{c}}}}{{\tau }_{\mathrm{nl}}}{\left(\displaystyle \frac{\lambda }{{d}_{i}}\right)}^{2/3}.\end{eqnarray} \tag{ 3 }$$

Note that the last factor is related to system size. Furthermore, if one associates the ion inertial scale with the turbulence Kolmogorov dissipation scale η, then effective Reynolds number is ${R}_{\mathrm{eff}}={(\lambda /\eta )}^{4/3}\to {(\lambda /{d}_{i})}^{4/3}$ and Equation (3) is equivalent to ${Q}_{p}/{Q}_{e}=({\tau }_{{\rm{c}}}/{\tau }_{\mathrm{nl}})\sqrt{{R}_{\mathrm{eff}}}$.

Figure 3 shows the scaling of relative proton–electron heating with ${\tau }_{{\rm{c}}}/{\tau }_{\mathrm{nl}}({d}_{i})$ according to the above hypothesis. Note that this ratio was also discussed by Matthaeus et al. (2014) as a useful measure of the influence of turbulence on kinetic processes such as linear Vlasov waves and linearized damping rates, which arguably require relatively undisturbed conditions for attaining their assumptions. Here, we are examining the complementary situation in which the disturbance of gyromotion by turbulent reconfiguration at ion scales is expected to promote proton heating. Note that the data plotted has many PIC simulations, including simulations in (Wu et al. 2013). We clearly observe that ${Q}_{p}/{Q}_{e}$ increases with ${\tau }_{{\rm{c}}}/{\tau }_{\mathrm{nl}}({d}_{i})$, supporting the hypothesis in Equation (3).

Zoom In Zoom Out Reset image size

The statistics examined so far have not specifically emphasized the role of coherent structures in the proton heating processes, yet theory indicates that coherent structures in the magnetic and plasma flow (Dmitruk et al. 2004; Markovskii et al. 2006; Parashar et al. 2011; Dalena et al. 2014; Wan et al. 2015) play an important role, as is confirmed in observations (Retinò et al. 2007; Osman et al. 2011). In particular, substantial fractions of heating are known to occur in boundary layer-like regions within or near proton scale current sheets (Ambrosiano et al. 1988; Dmitruk et al. 2004; Dalena et al. 2012; Drake et al. 2014). These structures typically have one or more dimensions on the order of an ion inertial scale (Dmitruk & Matthaeus 2006; Markovskii & Vasquez 2011; Wang et al. 2013; Makwana et al. 2015), which motivates our evaluation of the nonlinear time at d_i. The dynamical processes that occur in association with these structures include energization by direct (parallel) electric fields; strong scattering and confinement leading to stochastic orbits; betatron effect; and pickup by flows through ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift effect. These effects have been parameterized in various phenomenologies (Ambrosiano et al. 1988; Drake et al. 2006; le Roux et al. 2015) and it is likely that observed transfer of energy into random degrees of freedom, i.e., heating, occurs as a consequence of all of these acting in concert. The net effect of these processes, and notably their localization in structured patterns reflecting the magnetic structure of the turbulence, is evident in the simulations.

To illustrate the spatial intermittency and related properties of the heating, we plot T_p, T_e, ${T}_{p}/{T}_{e}$ and J_z in Figure 4. The Figure is based on run Turb817, which has ${\beta }_{i}={\beta }_{e}=0.1$, ${dZ}/{B}_{0}=0.8$ and $k=\{2,4\}$ (see Table 1).

Zoom In Zoom Out Reset image size

The electrons are hottest at the sites of strongest thin current sheets. Proton hotspots are relatively broader and are concentrated around the strong currents near magnetic "islands" and the interaction region between pairs of islands (Ambrosiano et al. 1988; Dmitruk et al. 2004; Chandran et al. 2010; Shay et al. 2014). (The hotspots where protons are hotter than electrons are mostly the spots where protons are strongly heated. These spots are the sites of strong currents around the islands.)

The above results lead to a perspective on plasma heating in which a von Karman, MHD-like heating rate is supported by highly nonuniform, intermittent dissipation processes that are related to coherent structure that form in or near boundaries at the edges of interacting magnetic flux tubes. There seems to be both a spatial and a temporal character to the conditions that favor strong proton heating. Evidently, proton heating is favored when the timescale for reconfiguration of the favored fluctuations becomes small enough that distortions occur within a gyroperiod. On the other hand the favored spatial scales are those that couple well to the proton population, i.e., the important fluctuations are those at or near the proton kinetic scales. In this way, the coupling of turbulent motions to heating might be viewed as a kind of generalization of resonance conditions; however, here with the interactions being with coherent structures rather than with wave trains and the time variations being of nonlinear origin and not simply oscillatory.

In the above analysis we employed the ion inertial scale d_i as the scale of interest, at which dissipative couplings to random motions occur, based on observations of the spectral breakpoints and current sheet thicknesses seen in solar wind, simulations and laboratory experiments. However, one also might have employed the thermal particle gyroradius ${\rho }_{i}$ as the preferred scale. In that case, the relevant nonlinear timescale would be ${\tau }_{\mathrm{nl}}{({\rho }_{i})={\tau }_{\mathrm{nl}}({d}_{i})({\rho }_{i}/{d}_{i})}^{2/3}={\tau }_{\mathrm{nl}}({d}_{i}){\beta }^{1/3}$, where $\beta ={c}_{{\rm{s}}}^{2}/{V}_{{\rm{A}}}^{2}$ (sound speed c_s; Alfvén speed V_A). Note that with the choice of gyroradius, the theory captures a better fit to the particle population; however, employing the ion inertial scale the theory makes better contact with the expected width of boundary layers, such as current sheet thicknesses. For a more general description of proton–electron heat functions, one might like to anticipate a more general β dependence, of the form ${Q}_{p}/{Q}_{e}=c(\beta )/{{\rm{\Omega }}}_{{\rm{c}}}{\tau }_{\mathrm{nl}}({d}_{i})$.

Based on the above analyses, we can now propose a description of proton heating:

$$\begin{eqnarray}&&{Q}_{p}=\displaystyle \frac{c(\beta ){\tau }_{c}/{\tau }_{{nl}}({d}_{i})}{1+c(\beta ){\tau }_{c}/{\tau }_{{nl}}({d}_{i})}\left[{C}_{K}\displaystyle \frac{{Z}^{3}}{\lambda }\right]\end{eqnarray} \tag{ 4 }$$

and electron heating:

$$\begin{eqnarray}&&{Q}_{e}=\displaystyle \frac{1}{1+c(\beta ){\tau }_{c}/{\tau }_{{nl}}({d}_{i})}\left[{C}_{K}\displaystyle \frac{{Z}^{3}}{\lambda }\right]\end{eqnarray} \tag{ 5 }$$

that are in accord with the theoretical and empirical descriptions given in the sections above. These prescriptions satisfy

• 1.  
  a von Karman scaling of the total heating rate ${Q}_{p}+{Q}_{e};$
• 2.  
  a scaling of the heating ratio ${Q}_{p}/{Q}_{e}$ with the total heating rate;
• 3.  
  a scaling of the proton–electron heat function ratio with the timescale ratio ${\tau }_{{\rm{c}}}/{\tau }_{\mathrm{nl}}({\ell })$, the latter evaluated at a scale ℓ in the range of d_i or ${\rho }_{i}$ controlled by the function $c(\beta )$.

Kinetic simulations (Servidio et al. 2012; Karimabadi et al. 2013) as well as test particle studies (Dmitruk et al. 2004) have shown that temperature anisotropies occur preferentially in the vicinity of coherent structures. That issue has not been examined here, although we fully expect that association is also present in the present simulations. These simulations were carried out in 2.5D (see Table 1) to maximize available spatial resolution. Although fully 3D studies would be preferred, we note that recent study Wan et al. (2015, 2016) showed a surprising degree of similarity between dissipation statistics in 2.5 and 3D runs of several types. It is also noteworthy that our simulations clearly exhibit coherent structures and intermittency (see Figure 4) but nonetheless are of limited system size (and effective Reynolds numbers). Therefore, we would expect that in still larger simulations the role of heating in or near coherent structures would be even greater.

This research was supported by NSF AGS-1063439, AGS-1156094 (SHINE), AGS-1460130 (SHINE), NASA grants NNX14AI63G (Heliophysics Grand Challenge Theory), and the Solar Probe Plus science team (ISIS/SWRI subcontract No. D99031L).

All the simulations used in this Letter were performed using the fully electromagnetic particle in cell code P3D (Zeiler et al. 2002) in periodic 2.5D (3 component vectors on a 2D grid) geometry. Simulations in the PAPJ15 family were also used to study the approach to MHD-like behavior in kinetic plasmas in (Parashar et al. 2015). Table 1 does not list simulations used in (Wu et al. 2013). Runs in the OTVdbB family are Orszag–Tang vortex runs with varying $\delta B/{B}_{0}$ amplitudes. Runs in the k13, k24lb and k2sb families were performed specifically for this study. These were initialized with MHD-like v and b fluctuations (solenoidal velocity and uniform density) excited at the wavenumbers within the k band given in the table with a specified spectrum, so that the initial energy is concentrated near the largest excited wavenumbers. The energy then cascades to fill in the complete wavenumber range available and produce turbulence. Typical parameters for these runs (except for Wu et al. 2013 runs) are given in Table 1.