CURRENT SHEETS AND COLLISIONLESS DAMPING IN KINETIC PLASMA TURBULENCE

Turbulence plays an important role in space and astrophysical plasmas by mediating the transfer of energy from large-scale motions to the small scales at which the turbulence can be dissipated. A major unsolved problem is the identification of the physical mechanisms that dissipate the small-scale turbulent motions, ultimately converting the turbulent energy to plasma heat. The dynamics at the dissipative scales are typically weakly collisional in diffuse astrophysical plasmas, such as the solar wind, so the mechanisms responsible for the dissipation and plasma heating are described by kinetic plasma physics. Two mechanisms have been proposed to be the dominantly involved in the dissipation process for plasma turbulence: collisionless wave-particle interactions (Howes et al. 2008a; Schekochihin et al. 2009; TenBarge et al. 2013a) and dissipation in small-scale current sheets (Dmitruk et al. 2004; Markovskii & Vasquez 2011; Matthaeus & Velli 2011; Osman et al. 2011; Servidio et al. 2011; Wan et al. 2012; Karimabadi et al. 2013).

In weakly collisional plasmas, it is well known that wave-particle interactions lead to significant collisionless damping of the linear kinetic wave modes. In turbulent astrophysical plasmas, it has been proposed that the fluctuations at perpendicular scales smaller than the ion Larmor radius, k_⊥ρ_i ≳ 1, have properties typical of kinetic Alfvén waves (Howes et al. 2008a, 2008b; Howes 2008; Schekochihin et al. 2009; Salem et al. 2012), and will therefore suffer collisionless damping. In the case of the solar wind, it has been suggested (Howes et al. 2008a, 2011a, 2011b; Schekochihin et al. 2009) that electron Landau damping dominates the dissipation of these turbulent electromagnetic fluctuations at k_⊥ρ_i ≳ 1. Free energy transferred conservatively to the particle distribution functions by wave-particle interactions is ultimately thermalized by arbitrarily weak collisions through the action of an entropy cascade in phase space (Schekochihin et al. 2009).

On the other hand, a number of recent studies focusing on the intermittent structures that inherently develop in plasma turbulence have suggested that dissipation dominantly occurs in coherent structures, in particular, small-scale current sheets (Dmitruk et al. 2004; Markovskii & Vasquez 2011; Matthaeus & Velli 2011; Osman et al. 2011; Servidio et al. 2011; Wan et al. 2012; Karimabadi et al. 2013). The kinetic physical mechanism by which dissipation occurs in current sheets has not been clearly elucidated. Hybrid kinetic-ion and fluid-electron simulations in two-dimensional (2D) suggest stochastic perpendicular ion heating due to demagnetization in current sheets (Parashar et al. 2009; Markovskii & Vasquez 2011), 2D and three-dimensional (3D) particle-in-cell simulations of reconnection suggest the acceleration of electrons by parallel electric fields (Drake et al. 2003; Pritchett & Coroniti 2004; Egedal et al. 2008, 2009, 2010, 2012) and/or Fermi acceleration (Drake et al. 2006), and 2D gyrokinetic simulations suggest linear phase mixing/Landau damping (Loureiro et al. 2013). Temperature measurements in the near-Earth solar wind have been used to both support (Osman et al. 2011, 2012) and refute (Borovsky & Denton 2011) the proposal that plasma heating dominantly occurs in current sheets.

In this Letter, we present a wave-driven, 3D gyrokinetic turbulence simulation at scales smaller than the ion Larmor radius that self-consistently generates small-scale current sheets. We find that the current sheet filling fraction is well correlated with the electron heating rate in the simulation. Yet, Ohmic dissipation is negligible, and the measured electron heating rate by scale is well reproduced by assuming dissipation is entirely associated with electron Landau damping, suggesting the possibility that the current sheets are damped collisionlessly by resonant wave-particle interactions and correspond to regions of local heating.

The simulation was performed with the Astrophysical Gyrokinetics code, AstroGK (Numata et al. 2010), which solves the equations of gyrokinetics (Frieman & Chen 1982; Howes et al. 2006). Collisions are treated using a fully conservative, linearized, and gyroaveraged collision operator (Abel et al. 2008; Barnes et al. 2009). The simulations are driven at the simulation domain scale with an oscillating Langevin antenna coupled to the component of the vector potential parallel to the equilibrium magnetic field, $\mathbf {B}_0=B_0 \hat{\mathbf {z}}$ (TenBarge et al. 2013b).

The simulation models a proton and electron plasma with a realistic mass ratio m_i/m_e = 1836, β_i = 1, and T_i/T_e = 1, where $\beta _i =v_{{\rm ti}}^2/v_A^2$, v_A is the Alfvén speed, and $v_{{\rm ti}} = \sqrt{2T_i/m_i}$ is the ion thermal speed. The simulation employs a periodic domain of size $L_\perp ^2 \times L_z$, elongated along the straight, uniform equilibrium magnetic field B₀. Relevant parameters are k_⊥ρ_i ∈ [5, 105], k_zρ_i/ ∈ [1, 16], A₀/ρ_iB₀ = 0.2, ν_iρ_i/v_ti = 0.2, and ν_eρ_i/v_ti = 0.5, where ρ_i = v_ti/Ω_i is the ion Larmor radius, Ω_i is the ion gyrofrequency, = 2πρ_i/L_z ≪ 1 is the gyrokinetic expansion parameter, A₀ is the antenna amplitude, and ν_s is the collision frequency of species s. Time is normalized to the linear frequency of a kinetic Alfvén wave at the simulation domain scale, ω₀ = 3.6k_z0v_A = 3.6ω_A0. Therefore, the corresponding domain scale turn-around time is $\tau _0 = 2\pi /\omega _0 \simeq 1.75\omega _{A0}^{-1}$. Collision frequencies are chosen to prevent build-up of small-scale structure in velocity space but remain small enough not to alter the weakly collisional dynamics: ν_s ≪ ω₀ is satisfied, so the simulation is weakly collisional. A value for ∼ δB/B₀ can be estimated by examining solar wind magnetic energy spectra at our simulation domain scale, kρ_i = 5. Based on spectra available in Alexandrova et al. (2009) and Sahraoui et al. (2009, 2010), ∼ δB/B₀ ∼ 0.01 at kρ_i = 5 in the solar wind. The antenna amplitude is chosen to satisfy critical balance at the domain scale, so the simulation represents critically balanced, strong turbulence. Analysis of a similar simulation (TenBarge et al. 2013a) demonstrates a magnetic energy spectrum in excellent agreement with the empirical form found from a large statistical sample of dissipation range measurements in the solar wind (Alexandrova et al. 2012).

First we demonstrate that the magnetic energy spectrum from our simulation is consistent with measurements in the solar wind and then present the perpendicular wavenumber spectrum of the square of the current density, j²(k_⊥). In Figure 1(a), we plot (solid) the average one-dimensional (1D) trace magnetic energy spectrum E_B(k_⊥), where the total magnetic energy $E_{B}^{({\rm tot})}= \int dk_\perp E_{B}(k_\perp)$. The average is performed over the steady-state evolution of the system, 1.5τ₀ ⩽ t ⩽ 4.1τ₀; error bars represent the variance over the same interval. This spectrum is quantitatively consistent with a large sample of measurements of the dissipation range magnetic energy spectrum in the solar wind (Alexandrova et al. 2012), suggesting that this simulation contains the essential physical ingredients underlying turbulence in the solar wind. Since j = |j| = |(c/4π)∇ × B|, we expect to find j∝kB, and therefore the scaling of j²(k_⊥) = ∫k_⊥dϕdk_∥j²(k) should satisfy the relation $j^2(k_\perp) \propto k_\perp ^3 E_{B}(k_\perp)$, as confirmed by the plotted spectrum (dashed) of j²(k_⊥).

Zoom In Zoom Out Reset image size

In a weakly collisional plasma with a guide magnetic field, current sheets are not expected to form at scales below k_⊥ρ_e ∼ 1. In Figure 2, we plot (a) the parallel current density j_z(x, y) from a perpendicular cut through the simulation domain. To explore the contribution to the current from different scales, we present spatially band-pass filtered data: (b) 5 ⩽ k_⊥ρ_i < 21, (c) 21 ⩽ k_⊥ρ_i < 84, and (d) k_⊥ρ_i ⩾ 84. The large-scale currents visible in panels (a) and (b) are dominated by the driving, which generates upward and downward propagating kinetic Alfvén waves with k_⊥ρ_i = 5. Panel (c), whose filter is approximately centered on the electron gyroradius, 0.49 ⩽ k_⊥ρ_e ⩽ 1.95, shows that this 3D gyrokinetic simulation indeed produces current sheets at the ρ_e scale, consistent with such development in a wide range of plasma turbulence simulations. These electron-scale diffusion regions are highly intermittent, both spatially and temporally, with a typical lifetime τ ≲ 0.1τ₀. The lack of significant current density in panel (d) shows that current sheets do not form at scales k_⊥ρ_e > 2. To confirm that our simulation indeed has sufficient perpendicular resolution to capture the current sheet dynamics, we ran a convergence test with double the resolution in the perpendicular plane. Plotted in Figure 1(b) are the instantaneous 1D trace magnetic energy spectrum for the manuscript resolution (solid) extending to k_⊥ρ_e = 2.5 and the same spectrum after doubling the resolution (dashed) to k_⊥ρ_e = 5 and allowing the simulation to saturate, demonstrating that the magnetic energy spectrum is resolved. The spatially filtered current density for the double resolution simulation shows similar results to Figure 2(d), that no current sheet structure forms at scales k_⊥ρ_e > 2, confirming that our simulation has sufficient perpendicular spatial resolution to capture fully the electron-scale current sheet dynamics.

Zoom In Zoom Out Reset image size

Note that we are using the general definition of current sheet as a discontinuity in the magnetic field. This is consistent with the definition of current sheet used in recent solar wind literature, e.g., Vasquez et al. (2007), Greco et al. (2009), and Osman et al. (2011). We do not attempt to differentiate between current sheets associated with reconnecting magnetic flux and those arising from interfering Alfvén waves.

It is also important to note that the simulation is driven by injecting domain-scale waves, which generate strong turbulence. The current sheets form self-consistently from the cascade of wave-driven turbulence and are not seeded or otherwise initialized to form. Therefore, we conclude that electron-scale current sheets form as a natural consequence of Alfvénic turbulence in this 3D gyrokinetic simulation.

We next examine the relative contribution to the measured electron heating from wave-particle interactions and dissipation in current sheets. The analytical equations for plasma heating (Howes et al. 2006) have been implemented as a diagnostic in AstroGK. Boltzmann's H Theorem states that the entropy increase necessary for irreversible heating requires collisions (Howes et al. 2006), so the collisional heating (plus a small amount of numerical dissipation) is used to measure the heating rate of each plasma species. Over the range of scales simulated, k_⊥ρ_i ∈ [5, 105], little ion heating occurs, so we focus on the electron heating. The collisional electron heating is given by

$$\begin{eqnarray} Q_e& =& -\sum _{\mathbf {k}_\perp } \int ^{L_z}_{-L_z} \frac{dz}{2 L_z} \int d^3 \mathbf {v} \frac{T_{0e}}{F_{0e}} \nonumber\\ & & \left[ h_{e\mathbf {k}_\perp }\left(\frac{\partial h^*_{e\mathbf {k}_\perp }}{\partial t} \right)_{{\rm coll}} + h^*_{e\mathbf {k}_\perp } \left(\frac{\partial h_{e\mathbf {k}_\perp }}{\partial t} \right)_{{\rm coll}}\right], \quad \end{eqnarray} \tag{ 1 }$$

where $h_{e\mathbf {k}_\perp }= h_e(k_x,k_y,z,v_\parallel,v_\perp,t)$ is the non-Boltzmann portion of the perturbed electron distribution function and F_0e is the equilibrium electron distribution function (Howes et al. 2006; Numata et al. 2010). The collisional electron heating as a function of perpendicular wavenumber Q_e(k_⊥) is computed by summing over annular rings in the perpendicular plane such that the total heating is given by Q_e = ∫dk_⊥Q_e(k_⊥). Note that, in the weakly collisional limit, the heating rate is independent of the collision frequency (Howes et al. 2006). For the steady-state evolution of the simulation over 1.5τ₀ ⩽ t ⩽ 4.1τ₀, the heating diagnostics recover total power balance to ≲ 2% (TenBarge et al. 2013a).

To estimate the contribution of current sheet dissipation to the heating rate, we compare the fraction of volume occupied by current sheets to the electron heating rate as a function of time. The volume filling fraction of current sheets is computed as the percentage of the volume with current density j > j_th, with a chosen threshold j_th = j_max/3, where j_max is the maximum current density over all time and space in the simulation. Varying this threshold alters the magnitude of the filling fraction but not the form of its variation with time. In Figure 3(a), we plot the volume filling fraction in percent (black dashed) and the electron collisional heating rate, discussed above, Q_e (red solid) as a function of time—all quantities in the figure have been integrated over the entire simulation domain. In Figure 3(b) is plotted the boxcar averaged (over Δt = 0.13τ₀) power injected into the plasma by the Langevin antenna (magenta dotted) and the total energy of the turbulent fluctuations in the simulation including the magnetic field and kinetic energies $E_{{\rm KAW}} = E_{B_\perp } + E_{B_\parallel } + E_{{\rm KE}}$ (blue dash-dotted), where $E_{{\rm KE}} = \sum _s m_s n_{0s} u_s^2 /2$ and u_s is the fluid velocity of each species.

Zoom In Zoom Out Reset image size

The cross correlations between the electron heating rate and filling fraction (black solid), antenna power (magenta dotted), and E_KAW (blue dash-dotted) are plotted in Figure 3(c). The injected antenna power and the electron heating rate are poorly correlated, 〈max(Corr(Q, j_ant))〉 = 0.52 ± 0.02—the mean and variance are calculated from the present simulation and five other identically prepared simulations employing different random number seeds for the turbulent driving. Similarly, the electron heating rate and total turbulent energy are not well correlated, 〈max(Corr(Q, E_KAW))〉 = 0.78 ± 0.04, suggesting that the heating rate is not simply a function of the driving or magnitude of turbulent energy. On the other hand, the electron heating rate is well correlated with the current sheet filling fraction, 〈max(Corr(Q, n_fill))〉 = 0.91 ± 0.04. The strong correlation suggests that dissipation associated with current sheets plays an important role in heating the electrons.

In Figure 4, we present a plot of the electron collisional heating rate by perpendicular wavenumber, Q_e(k_⊥) (solid), averaged over 1.5τ₀ ⩽ t ⩽ 4.1τ₀, where error bars represent the variance over the interval. The instantaneous shape of the heating curve is similar to the average. This plot shows that the electron heating is nearly constant over all scales, with about half of the total heating occurring at scales k_⊥ρ_e < 1. The turn-down at k_⊥ρ_i > 105 is an artifact due to the diminishing number of Fourier modes in the corner beyond the fully resolved simulation domain.

Zoom In Zoom Out Reset image size

As a function of k_⊥, we may predict the collisionless damping of the turbulent fluctuations by resonant wave-particle interactions in our simulation using Q_wp(k_⊥) = 2γ(k_⊥)E_KAW(k_⊥), where γ is the linear kinetic damping rate of kinetic Alfvén waves (dominated by electron Landau damping). This prediction for the wave-particle interaction heating rate requires integration over k_∥, where parallel is with respect to the local magnetic field and is typically determined via structure functions or wavelets. To avoid the complications of determining the local magnetic field direction, we use frequency as a proxy for the parallel wave vector since ω∝k_∥ for kinetic Alfvén waves (TenBarge & Howes 2012). This prediction for wave-particle damping, plotted in Figure 4 (dotted), admits no free parameters, yet it agrees well with the measured collisional electron heating (solid): the integrated, total predicted electron heating is within 4% of the collisional heating diagnosed in AstroGK. The slight excess of wave-particle damping at 5 < k_⊥ρ_i < 40 and of electron collisional heating at k_⊥ρ_i > 40 is consistent with the action of the electron entropy cascade (Schekochihin et al. 2009). Through the entropy cascade, energy removed by electron Landau damping at 5 < k_⊥ρ_i < 40 is expected to appear as collisional heating at higher wavenumbers.

We also compute the Ohmic heating rate Q_η = ηj², where $\eta = 0.38(4\pi)\nu _{{\rm ei}} d_e^2/c^2$ is the Spitzer resistivity (Spitzer & Härm 1953), ν_ei = ν_e is the electron–ion collision frequency, d_e = c/ω_pe is the electron inertial length, and ω_pe is the electron plasma frequency. The Ohmic heating rate is plotted (dashed) in Figure 4, clearly Q_η ≪ Q_wp ≃ Q_e. Theory predicts negligible Ohmic heating for a weakly collisional plasma, since electron–ion collisions are insufficient to significantly heat the electrons. Therefore, Ohmic dissipation of the current cannot account for the observed electron heating in the simulation, despite the strong correlation between current sheet filling fraction and heating rate.

A puzzling aspect of these results is that, despite the clear correlation between the electron heating rate and the volume filling fraction of current sheets, the electron heating as a function of wavenumber is well predicted assuming that Landau damping is entirely responsible for the electron heating. This unexpected agreement raises the interesting possibility that the dissipation of the current sheets in the simulation occurs entirely via Landau damping.

The solution to this puzzle lies in the relationship between the current and magnetic field, namely j∝kB. This relationship implies that regions of strong current correspond to regions with enhanced small-scale magnetic structure since the current is weighted toward small spatial scales. Since the linear kinetic damping rate increases with wavenumber, regions with enhanced small-scale magnetic structure will also correspond to regions of enhanced wave-particle damping. Therefore, regions of strong current may also correspond to regions of enhanced wave-particle damping.

To test this hypothesis, we apply a high-pass filter to the data presented in Figure 3 retaining only modes with k_⊥ρ_i ⩾ 10—no filter is applied to the antenna current. The result of the filtering is plotted in Figure 5, where all of the cross correlations exceed 0.9, with the exception of the antenna power, which remains poorly correlated.

Zoom In Zoom Out Reset image size

The picture suggested by this simulation is one in which current sheets are self-consistently formed by the interaction of kinetic Alfvén wave-like fluctuations, where each of the fluctuations is dissipating at its Landau damping rate. Therefore, current sheet formation and dissipation is dominated by the evolution of the Alfvénic turbulence, and current sheets correspond to sites of enhanced dissipation and heating, as suggested by recent analyses of solar wind turbulence (Osman et al. 2011, 2012).

The validity of the Landau prescription described in this manuscript is predicated on the distribution function not deviating significantly from Maxwellian. This can be simply tested by examining the electron fluid velocity moments. We find that max(u_ze/v_te) ∼ 0.01 based on the value of determined in Section 2, suggesting that the Landau prescription is indeed valid. A more detailed examination of the electron distribution functions within regions of intense current in our simulations is necessary to confirm the dominance of the dissipation by electron Landau damping.

We find that electron scale current sheets develop self-consistently as a consequence of wave-drive turbulence in a 3D gyrokinetic simulation of the dissipation range, k_⊥ρ_i ∈ [5, 105]. The electron heating rate is well correlated with the volume filling fraction of current sheets, suggesting that the dissipation of current sheets plays an important role in the heating of electrons. However, the electron heating rate as a function of scale is well predicted by assuming that all dissipation is due to collisionless damping of the turbulent fluctuations via Landau resonance with the electrons. In the weakly collisional plasma, Ohmic dissipation of current sheets is negligible. Due to the relationship between the current and magnetic field, significant current highlights regions with enhanced small-scale magnetic structure, which will be collisionlessly damped at a rate greater than surrounding plasma. This suggests that current sheets may correspond spatially to locations of enhanced dissipation and heating, regardless of whether that dissipation is due to collisionless wave-particle interactions or active magnetic reconnection.

The authors thank Homa Karimabadi. This work was supported by NASA grant NNX10AC91G and NSF CAREER grant AGS-1054061, resources of the Oak Ridge Leadership Computing Facility, supported by the Office of Science of the U.S. DOE under Contract No. DE-AC05-00OR22725, and advanced computing resources provided by the NSF, partly performed on Kraken at the National Institute for Computational Sciences.