Energy transfer, pressure tensor and heating of kinetic plasma

Kinetic plasma turbulence cascade spans multiple scales ranging from macroscopic fluid flow to sub-electron scales. Mechanisms that dissipate large scale energy, terminate the inertial range cascade and convert kinetic energy into heat are hotly debated. Here we revisit these puzzles using fully kinetic simulation. By performing scale-dependent spatial filtering on the Vlasov equation, we extract information at prescribed scales and introduce several energy transfer functions. This approach allows highly inhomogeneous energy cascade to be quantified as it proceeds down to kinetic scales. The pressure work, $-\left( \boldsymbol{P} \cdot \nabla \right) \cdot \boldsymbol{u}$, can trigger a channel of the energy conversion between fluid flow and random motions, which is a collision-free generalization of the viscous dissipation in collisional fluid. Both the energy transfer and the pressure work are strongly correlated with velocity gradients.


I. INTRODUCTION
The classical energy cascade scenario is of great importance in explaining the heating of corona and solar wind [1][2][3][4][5]: In these applications, one envisions that significant amounts of energy reside in large-scale fluctuations. Nonlinear interactions cause a cascade that transports energy to smaller scales where dissipation occurs. At dissipation scales kinetic processes that absorb these energy fluxes, produce temperature enhancements. Inhomogeneities generated in this process may in turn be responsible for large scale flows, while also producing populations of energetic particles. Many systems, including low-collisionality astrophysical plasmas, may be well described by fluid theory at large scales. For plasmas, including electromagnetic fields, this large scale description would usually be taken to be some form of magnetohydrodynamics (MHD). The present paper is devoted to extending ideas about turbulence cascade into the deep kinetic range, so that we may develop a better understanding of turbulence cascade in low collsionality or even collsionless plasma.
It is reasonable to assume that MHD remains a credible approximation for a kinetic plasma at scales large enough to be well separated from kinetic effects. For that range, scale filtered analysis [6] shows how the fluid flow energy cascades almost conservatively from large to small scales, despite not being strictly an invariant of the MHD system. Then an important subsequent question asks how energy transfer proceeds down to kinetic scales as various kinetic processes come to the fore. The present paper addresses statistical properties of energy transfer across scales, recognizing the possible significance of energy cascade in explaining the heating and acceleration of the wind, and many other properties.
In this work we will avoid adopting familiar approaches that rely heavily on linear theory of waves, instabilities and damping rates [7][8][9], or on weak turbulence approaches that require a leading order description in terms of linear modes. Instead we consider the full Vlasov-Maxwell system, and employ a filtering approach that is familiar in hydrodynamics [10] and large-eddy simulation [11][12][13][14][15] communities but less used in kinetic plasma. Examining filtered equations for energy transfer, we can assess the relative importance of different transfer terms at all scales ranging from MHD to electron scales.
The present approach provides extensions of what fluid models tell us about the plasma cascade. In the context of plasma applications, MHD simulations adopt an ad hoc model of dissipation (e.g., viscous and resistive dissipation), rather than engaging the details of the small dynamics that make up the plasma dissipation range. The turbulence in most astrophysical contexts, on the other hand, is typically of weak collisionality, and frequently modeled as collisionless, and thus collisional (viscous and resistive) dissipation at small scales cannot emerge immediately. While various specific processes may contribute to conversion of energy from fields into random degrees of freedom, for example, wave-particle interactions (WPI) [16][17][18][19][20] and processes associated with coherent structures (CS) [21][22][23][24][25][26], are likely ingredients, but nevertheless an explicit dissipation function cannot at this moment be defined clearly for a collisionless system.
Lacking such an explicit form for dissipation, Wan et al. [27] considered a surrogate dissipation measure related to the work done by the electromagnetic field on the plasma particles. Recent studies in compressible MHD turbulence [6,28] demonstrated that, apart from collisional dissipation, the pressure dilatation, −p∇ · u, can trigger an alternative channel of the conversion between kinetic and internal energy. Accordingly, one could expect that there might be an analogous role of the pressure tensor in collisionless plasma. In fact one expects pressure to be influential in at least several ways. On the one hand, the pressure term in anisotropic compressible turbulence moderates the competition and balance between two energy redistributive processes, i.e., return-to-isotropy [29,30] and kinetic-potential (internal) energy equipartition [31][32][33][34]. On the other hand, the pressure tensor in kinetic plasmas plays a very important role in the force balance equations as well as in the generalized Ohm's law near neutral lines [35,36]. Here we show that the global energy exchange between fluid flow and particles (i.e., kinetic and thermal energies), derived from the Vlasov equation, is bridged immediately by the collaboration of pressure tensor and strain stress (i.e., velocity gradient). This possibly provides a new perspective on the collisionless dissipation mechanism and on the collisionless plasma cascade in general. The possible importance of pressure work in generating internal energy has been brought up in Ref. [37] (hereafter, Paper I). Here we extend that study and explore this novelty in a more comprehensive and detailed way.

II. GLOBAL ENERGY CONVERSION
Standard manipulation of the Vlasov equation yields macroscopic equations for plasma particles of type α in a collisionless plasma: Here ρ α = m α n α represents the mass density; m α is the mass of particles of species α; n α is the number density; u α gives the fluid flow (bulk) velocity; n α q α represents the charge density; dv is the total (average and random) kinetic energy; Decomposing the total energy E α into average and random parts facilitates the understanding of energy converting processes. On defining the species fluid flow energy as Computing the inner product of Eq. 2 with u α results in the fluid flow energy equation: Subtracting Eq. 4 from Eq. 3 we obtain a time evolution equation for the random kinetic energy, Using the Maxwell curl equations, the equation governing electromagnetic energy, E m = 1 8π (B 2 + E 2 ), can be written as: where j = α j α is the total electric current density, and j α = n α q α u α is the electric current density of species α. Under certain boundary conditions, e.g., periodic, integrating Eqs. 4, 5, and 6 over the whole volume, we can have where · · · denotes a space average over the entire volume. values of − (P α · ∇) · u α and j α · E are not sign-definite. Therefore, there are two possible directions of energy conversion. · · · denotes the space average over the entire volume. − (P α · ∇) · u α , is the only term converting fluid flow energy into thermal (random) energy, while the term, j α · E , represents the conversion between fluid flow and electromagnetic energies. The pressure work seems to be a more straightforward measure of heating rate when compared with electromagnetic work [27,38]. At present, we cannot rule out the possibly strong correlation between the work done by pressure and work done by the electric field.
For example, for a generalized Ohm's law, or the electron momentum equation, in the limit of massless electrons in collisionless plasma, we find that j e · E = − (P e · ∇) · u e .
Notwithstanding that the pressure work is a general property in various fluid systems, seldom have studies investigated the role of pressure in modifying the thermal (random) energy in a turbulent kinetic plasma (however, see, e.g., Birn and Hesse [39], Birn et al. [40], Birn and Hesse [41]). In the realm of observations this is more or less due to the intractability of calculating velocity gradient from single spacecraft datasets and until recently, a lack of high cadence determination of the pressure tensor. These complication have led most observational studies of solar wind turbulence to rely on high cadence magnetic field data, which is generally much more accessible. Even in simulation studies, accurate determination of the pressure tensor is challenging, requiring either large numbers of particles in PIC codes, or the use of computationally demanding Eulerian Vlasov simulations. Here we will use numerical simulations to explore the role of pressure tensor in heating of kinetic plasma in detail.

III. SIMULATION DETAILS
The fully kinetic particle-in-cell (PIC) simulation employed here spans from macroscopic fluid scales to kinetic scales. It is expedient for studies of energy transfer and dissipation.
The simulation was performed using P3D code [42] in 2.5D geometry (three components of dependent field vectors and a two-dimensional spatial grid). Number density is normalized to the reference number density n r (=1 in this simulation), mass to proton mass m i (=1 in this simulation), and magnetic field to B r (=1 in this run). Length is normalized to the ion inertial length d i , time to the ion cyclotron time Ω −1 i , and velocity to the reference Alfvén speed v Ar = B r / (4πm i n r ) 1/2 . Parameters of the PIC simulation are listed in Table   I. It is conducted in periodic boundary conditions in both directions. The initial density, temperature and out-of-plane magnetic field B 0 are uniform. The initial v and b fluctuations are transverse to B 0 ("Alfvén mode") and excited for prescribed wavenumbers with specified spectra and cross helicity. More details about the simulation can be found in Ref. Wu et al. [43]. TABLE I: Simulation parameters: box size L, grid points N 2 , mass ratio m i /m e , proton beta β i , electron beta β e , out of plane uniform magnetic field B 0 , the number of particles of each species per grid ppg and correlation scale λ c .
We analyze statistics near the time of maximum root mean square(r.m.s.) electric current density (i.e., tΩ i = 206.25), shown by the dashed lines in Fig. 2(a). By this time the turbulence has fully developed and it has generated a complex pattern, characterized by sheet-like current structures spanning a range of scales as seen in Fig. 2(b). In order to remove particle noise at grid scales, we apply a low-pass Fourier filter to the field data for kλ d > 1 (Debye length λ d ) prior to statistical analyses. Demonstration of the efficacy of this filtering technique will be documented elsewhere [44].

IV. ROLE OF PRESSURE TENSOR
A somewhat surprising result [37] further documented here is the correlation of pressure effects and dynamically appearing coherent structures associated with intermittency. In order to clarify the correlation between pressure tensor and intermittent structures, we split the pressure work, − (P · ∇) · u, into two parts as follows. We suppress the subscript α for simplicity.
Decomposing the pressure tensor, P = (P ij ), gives where p = 1 3 P ii is the scalar pressure; a sum on repeated indexes is implied; δ ij is the Kronecker delta; Π ij is the remaining part with the scalar pressure subtracted from the pressure tensor, that is, the deviatoric pressure tensor.
The intrinsic decomposition of ∇u into symmetric and anti-symmetric parts gives are the strainrate and rotation-rate tensors, respectively; ω i is the vorticity; θ = S ii is the dilatation; 3 θδ ij is the traceless strain-rate tensor. Using Eqs. 10 and 11, we can obtain One can see that −pθ, accounting for compressive effects, is of the same form as the pressure dilatation in compressible MHD turbulence. The new term in kinetic plasma in comparison with MHD is −Π ij D ij , the double contraction of deviatoric pressure tensor and traceless strain-rate tensor. This term which will be called "Pi-D" hereafter is a salient feature of the present paper.
At this point we may recall that in the continuum formulation leading to the Navier-Stokes equations, with strong collisions, −Π ij actually is present, but rarely written in this way. Instead, in a Chapman-Enskog development, this term is equated with a viscous stress, which can then be expressed in terms of velocity gradient [45,46]. In Sec. IV A, we will check whether the "Pi-D" term shares some properties with viscous dissipation.
average over about an electron gyroperiod is used in computing E · j α to eliminate very high frequency oscillations.
To clarify the effect of the pressure tensor, we plot in Fig. 3   Q ω = 1 4 ω 2 / ω 2 , and the normalized second invariant of the traceless strain-rate tensor (bottom row), Left panels: electrons (reproduced from Paper I for completeness); and right panels: ions.

A. Energy conversion related to coherent structures
Various studies based on numerical simulations [24,25,27,[47][48][49][50] and solar wind data [23,[51][52][53][54] support the idea that enhanced kinetic activity, such as temperature anisotropy, heating, particle acceleration, and departures from Maxwellian velocity distributions in general, all of which commonly observed in astrophysical and laboratory plasmas, are strongly inhomogeneous. These effects are associated typically with coherent structures such as magnetic structures. Indeed, intense kinetic activity is often found in the general proximity to strong gradients, including not only magnetic, but also strong density and velocity gradients [16,[55][56][57][58][59][60]. In particular, Refs. [58,60,61] find that heating is correlated with both vorticity and current density, but more strongly with vorticity.
Here we are interested in looking at velocity gradients and their interaction with the pressure tensor, in view of the important role of these quantities in energy conversion in kinetic plasma, as seen in Eq. 5. We base our diagnostics on the geometric invariants of the relevant second-order tensorial quantities, an approach extensively employed in hydrodynamics to describe flow patterns [62][63][64][65]. Based on the decomposition in Eq. 11, the normalized second invariants of the traceless strain-rate matrix (D ij ) and the rotation-rate matrix (Ω ij ) are Q D = 1 2 D ij D ij / 2D ij D ij and Q ω = 1 4 ω 2 / ω 2 , respectively. We can also define a similar quantity for the electric current density, say, Q j = 1 4 j 2 / j 2 . Figs. 4 and 2(b) show contours of Q ω , Q D and Q j , which are found to be non-uniformly distributed in space. Moreover, comparison of Figs. 4, 2(b) with Fig. 3(a) reveals greatly similar patterns of Pi-D, Q ω , Q D and Q j , as also shown in [60]. This indicates, as also described in Paper I, that these intermittent structures might be sites of enhanced energy conversion, which is consistent with various recent results [21-26, 57, 61].
We highlight the possible correlations of Pi-D, Q ω , Q D and Q j by plotting them on cuts along the X direction. A sample of the absolute values of "Pi-D" terms along the X direction with Y 35d i is shown in Fig. 5. One sees the spatial distributions of "Pi-D" terms are evidently bursty, suggesting a connection to the spatial intermittency of the turbulence. A useful intermittency measure is given by the partial variance of increments (PVI) [48], PVI(f ) = |∆f | √ |∆f | 2 , where ∆f = f (s + ∆s) − f (s). Here we choose a small scale lag, ∆s 0.2d i = d e . The PVI series of bulk velocities for electrons (u e ) and ions (u i ) and magnetic field (b) in Fig. 5 behave quite similar to "Pi-D", with peak values in the vicinity of high "Pi-D". Events with high "Pi-D" on the series are selected by labeling with sequential numbers in the figure, and also displayed along with Q ω , Q D and Q j in Fig.   6. We might expect that intense kinetic activities are associated with current sheets, in particular, with high values of current density. In Fig. 6, however, the energy conversion through "Pi-D" is more correlated with Q ω and Q D in comparison with Q j . This can be seen  [56], the kinetic heating of protons might be a "viscous like" process instead of a magnetic process. A possible clue comes from recent work [59] that describes in detail how the eigenvalues of the velocity gradient tensor contribute to the time derivative of the pressure tensor (see, e.g., Eq.(15) of Ref. [59]). In the light of this, one can expect "Pi-D" to be correlated with Q D more strongly in comparison with Q ω , as is the case for ions in the right panels of Fig. 6.
The net effect of the deviatoric pressure work in our case is found to increase internal energy as was done by dissipation in collisional fluids. This provides a possible pathway when collisional relaxation is either absent or weak in kinetic plasma to understand how part of the fluid flow energy is converted into thermal (random) energy. To arrive at this perspective, there is no need to specify a particular mechanism, such as reconnection heating, Landau damping, cyclotron damping, or stochastic heating. Since the above commentary is based on the Vlasov equation alone, any process that contributes to the net "Pi-D" interaction is contributing to thermal energy increase. We have also not attempted to address the question of what physical processes make such energy exchanges irreversible. Hence we cannot comment on how much of the energy converted into random motions can settle down to heat permanently. With these caveats in mind, we defer these questions to future study.

B. Joint PDFs of coherent structures
Several types of coherent structures emerge in turbulent flow, such as Q D reflecting straining motions, Q ω corresponding to rotation and Q j related to magnetic discontinuities.
All of them can interact with one another. Fig. 7 shows the joint PDFs P (x, y).
Generally, in many hydrodynamic turbulent flows, the dominant structures are found to be tube-like structures, like vortex tubes with concentrated enstrophy (mean square vorticity), whereas sheet-like structures between these tubes are regions of dissipation. It is therefore expected that strain rate is not correlated with rotation rate. In hydrodynamic turbulence away from walls [67,68], the joint PDF of Q D versus Q ω is spread very broadly.
However, here the joint PDFs P (Q ω , Q D ) in Fig. 7 are dominated by a population near the Q D = Q ω line, which demonstrate the strong correlation between these two quantities. It also indicates a frequently occurring class of sheet-like rather than tube-like structures, a feature consistent with many visualizations of MHD turbulence [69]. We also find that the correlation between Q D and Q ω for electrons (see Fig. 7(a.1)) is stronger than that for ions (see Fig. 7(a.2)). Meanwhile, Fig. 6 confirms that the curves of Q D and Q ω for electrons almost overlap, while those for ions are more dispersed.
The joint PDFs of Q j versus Q D and Q ω , shown in Fig. 7, are spread rather broadly, indicating weak pointwise correlation between these quantities. By examining the curl of the Lorentz force, one may deduce [60,70] that the vorticity structures tend to form on the flanks of strong current structures. Therefore, the vorticity distribution does not correlate exactly with the current, but are somewhat offset in space, but nevertheless located in the vicinity of the current. This can be also observed in Fig. 6.

V. ENERGY TRANSFER ACROSS SCALES
Now that the role of pressure tensor in global energy conversion is established, we turn to the hierarchy of scales and seek to find clues to how energy cascade proceeds from MHD scales down to the kinetic scales. A low-pass spatial filter approach [71] eliminates information at small scales without affecting remaining information at large scales, thus providing a powerful technique to study cross-scale energy transfer.
Extension of the filtering approach to kinetic plasma is technically easy. More derivation FIG. 7: Joint PDFs P (x, y) of the normalized second invariants of rotation-rate, traceless strain-rate tensors, and current density, i.e., Q ω = 1 4 ω 2 / ω 2 , Q D = 1 2 D ij Dij/ 2D ij D ij , and Q j = 1 4 j 2 / j 2 , for both electrons (left panel) and ions (right panel). The correlation between Q ω and Q D is strong, for both electrons and ions, indicating sheet-like vortex structures. Meanwhile, this correlation for electrons is stronger than that for ions. The correlation between Q ω (or Q D ) and Q j is weak. details can be found in App. A. Note that for conciseness of notation, the filtering scale is not written explicitly unless necessary for clarity. The equation of filtered fluid flow energy, E f α =ρ αũ 2 α /2, can be written as The meaning of each term can be understood referring to Eqs. 7-9 as: , is the flux of the fluid flow energy transfer across scales (if the filter scale is set to zero, these terms vanish); Φ uT α = − P α · ∇ ·ũ α is the fluid flow energy converted into thermal (random) energy due to pressure work; Λ ub α = −q αnα E ·ũ α is the rate of fluid flow energy conversion into electromagnetic energy, i.e., electromagnetic work done on the fluid (seen more clearly in the filtered equation for electromagnetic energy in App. B).

A. Energy fluxes
We show energy fluxes, Π uu α , Φ uT α and Λ ub α varying with filtering scales in Fig. 8. The positive Π uu α transfers fluid flow energy from large to small scales (i.e., a forward cascade) due to the interaction of sub-grid scales with large scales, such as the sub-scale stresses The flow energy cascade proceeds from the largest scales (here we calculate up to the correlation scale) and transfer is sustained down to electron inertial scale d e . The classical theory [72] in incompressible hydrodynamics suggests an energy cascade where energy is transferred from large to small scales at a constant rate. This issue has been studied in compressible HD and MHD, and the SGS energy flux is found to be approximately constant; see [6,10,73,74]. The idealized notion of a constant cascade rate is not applicable here, likely due to the limited scale separation, equivalent to a relatively low Reynolds number.
Unlike the term Π uu α associated with interactions between scales > and < , the terms Φ uT α = − P α · ∇ ·ũ α and Λ ub α = −q αnα E ·ũ α incorporate information only from scales > . Therefore, they are cumulative quantities. Note that Φ uT α is vanishingly small In contrast, Λ ub α , the filtered contribution to −j α · E , is fairly constant over kinetic scales, and the observed increases are concentrated at a few times ion inertial scale d i .
Therefore contributions to the energy conversion between fluid flow and electromagnetic fields mainly result from the large scales in this simulation. Simulation with larger size is required to conclude if this transfer is macroscopic. We schematically represent this qualitative analysis in Fig. 9.

B. Location of energy transfer enhancement
In analogy with the pressure work Pi-D in Fig. 3(a) and coherent structures in Figs. 4, we portray spatial contours of energy flux Π uu α in Fig. 10. It is apparent that the energy flux is highly inhomogeneous. Comparing these figures, it is apparent that there is quite a significant coincidence between the coherent structures and the sites of enhanced energy transfer. Arguments have been made in favor of the hierarchy of coherent structures (see [50]). It also appears clearly in Fig. 10. In comparing the top row with ∼ 2d e and the second invariant of the traceless strain-rate tensor Q D = 1 2 D ij Dij/ 2D ij D ij is supported by the substantially similar patterns of contour maps and cuts along the X direction, which corroborates again the idea suggested in Paper I. In collisional fluid, Q D characterizes the viscous dissipation. The energy transformation by −Π ij D ij in our case behaves like viscous dissipation. We note that there is not any known quantitative theory to explain this strong dependence rigorously as far as we are aware. As described, for example, by Del Sarto et al. [59] and Cerri et al. [75], velocity shear provides an effective mechanism that can cause an initially isotropic pressure to become nongyrotropic, which can trigger a rich variety of inhomogeneous instabilities and nonlinear effects [57]. This might be of some help to understand the result. Further confirmation, based on different codes, initial conditions, parameters and so on, is needed to assert any degree of universality of this result. Once it is clarified, we might find a pathway (perhaps a "viscous-like collisionless dissipation") by which various processes act in concert and finally lead to intermittent heating. It is noteworthy that this result is in some ways complementary to the recent finding [27,38] that in several types of kinetic simulations of plasma turbulence, conditional averages of j ·E based on local values of j behave very similarly to what would be expected if j ∝ E, which is exact for MHD, but does not emerge in any obvious way for a kinetic plasma. These results represent a significant challenge to understand based on theory, and also provide guidance for developing phenomenological theories of plasma turbulence and dissipation.
Here we also introduced several energy transfer functions using filtering approaches to analyze their statistical properties. This approach affords a path to understand how classical cascade theory is extended from MHD scales down to kinetic scales. One can envision that energy exchange between fluid flow and magnetic and electric fields through E · j α occurs at large scales, and a part of fluid flow energy is converted into internal energy through − (P α · ∇) · u α at small scales. These two conversions are bridged by the fluid flow energy cascade at moderate scales.
Another important aspect of plasma turbulence and energy conversion that we have described here is the association between energy cascade, conversion processes and coherent structures. It transpires that energy transfer and conversion are in general inhomogeneous localized processes, correlated with several types of structures, such as velocity gradients (i.e., symmetric straining and vortex) and current sheets. Conversion into random energy, in effect, dissipation, is more strongly correlated with velocity gradients. This situation is actually very complex, involving a high correlation between symmetric straining and vorticity, which requires special conditions such as sheet-like structures. Conversely a weak correlation is found between current sheets and vorticity (or symmetric straining), even though these quantities are concentrated in juxtaposed regions of space. Clearly, more work is required to reveal the dynamical interactions among these coherent structures. For this, we might gain some inspiration from Refs. [60,70] or the evolution equation of enstrophy for incompressible MHD, d t 1 2 ω 2 = ω · D · ω + νω · ∆ω + ω · ∇ × (j × B) The 2.5D simulation used in this paper presents some limitations with regard to potential generality of the conclusions. In fact the present analysis might miss some important effects present in 3D geometry (in spite of recently reported qualitatively similarities [38,76,77]).
The intent here is not to accommodate all the dynamical complexity of the astrophysical plasma, but to propose an alternative pathway to study turbulence and heating in a kinetic plasma. Further development of these ideas may be useful, e.g., in theoretical study of coronal heating and acceleration of solar wind. Finally, neither the detailed differences between electron and ion dynamics, nor its implications for heating in space plasmas, are discussed here, but rather deferred for later study.
One can see that these equations can also be derived immediately from the macroscopic Eqs.
Taking the electromagnetic work (i.e., Λ ub α in Eqs. 13 and A10) as a source in Eq. B1 yields which can be written as where J b = c 4π E × B is the spatial transport; Π bb α = −q αnατ e α ·ũ α , whereτ e α = E − E , is the flux of electromagnetic energy across scales due to sub-scale work done by the electric field; Λ ub α = −q αnα E ·ũ α , the rate of fluid flow energy conversion into electromagnetic energy, is the same as those in Eqs. 13 and A10, which is canceled out on combining Eqs. 13 (or A10) and B1. Therefore, the filtered equation for total fluid flow and electromagnetic energy takes the form where J u α = E f αũ α +ρ ατ u α ·ũ α + P α ·ũ α , Π uu α = − (ρ ατ u α · ∇) ·ũ α − q αnατ b α ·ũ α and Φ uT α = − P α · ∇ ·ũ α are defined in the text following Eqs. 13 and A10.