Energy conversion in turbulent weakly-collisional plasmas: Eulerian Hybrid Vlasov-Maxwell simulations

TL;DR

This paper uses Eulerian Hybrid Vlasov-Maxwell simulations to analyze energy conversion processes in turbulent weakly-collisional plasmas, revealing the role of pressure-strain interactions and local structures in proton kinetic effects.

Contribution

It demonstrates that Eulerian HVM simulations can accurately track energy transfer and conversion in plasma turbulence, providing new insights into pressure-strain interactions and local structure correlations.

Findings

Pressure-strain interaction drives proton temperature growth.

Eulerian simulations are nearly noise-free and accurate.

Strong-current, high-vorticity structures are linked to pressure-strain effects.

Abstract

Kinetic simulations based on the Eulerian Hybrid Vlasov-Maxwell (HVM) formalism permit the examination of plasma turbulence with useful resolution of the proton velocity distribution function (VDF). The HVM model is employed here to study the balance of energy, focusing on channels of conversion that lead to proton kinetic effects, including growth of internal energy and temperature anisotropies. We show that this Eulerian simulation approach, which is almost noise-free, is able to provide an accurate energy balance for protons. The results demonstrate explicitly that the recovered temperature growth is directly related to the role of the pressure-strain interaction. Furthermore, analysis of local spatial correlations indicates that the pressure-strain interaction is qualitatively associated with strong-current, high-vorticity structures, although other local terms -- such as the heat flux -- weaken the correlation. These numerical capabilities based on the Eulerian approach will enable deeper study of transfer and conversion channels in weakly collisional Vlasov plasmas.

Equations

$$\frac{\partial f_p}{\partial t}+\mathbf{v}\cdot \frac{\partial f_p}{\partial \mathbf{x}}+\left(\mathbf{E}+\mathbf{v}\times \mathbf{B}\right)\cdot \frac{\partial f_p}{\partial \mathbf{v}}=0,$$

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \mathbf{E}$$

$$\mathbf{E}=-{\mathbf{u}}_p\times \mathbf{B}+\frac{\mathbf{j}\times \mathbf{B}}{n_p}-\frac{\nabla P_e}{n_e}+\eta \mathbf{j}$$

$$\frac{\partial {\mathcal{E}}_p^f}{\partial t}+\nabla \cdot \left({\mathbf{u}}_p{\mathcal{E}}_p^f+{\mathbf{u}}_p\cdot {\mathbf{P}}_p\right)=\left({\mathbf{P}}_p\cdot \nabla \right)\cdot {\mathbf{u}}_p+n_p{\mathbf{u}}_p\cdot \mathbf{E}$$

$$\frac{\partial {\mathcal{E}}_p^{th}}{\partial t}+\nabla \cdot \left({\mathbf{u}}_p{\mathcal{E}}_p^{th}+{\mathbf{q}}_p\right)=-\left({\mathbf{P}}_p\cdot \nabla \right)\cdot {\mathbf{u}}_p$$

$$\frac{\partial {\mathcal{E}}^m}{\partial t}+\nabla \cdot \left(\mathbf{E}\times \mathbf{B}\right)=-\mathbf{j}\cdot \mathbf{E}$$

$${\mathbf{P}}_p=\frac{1}{2}\int d^{3}v(\mathbf{v}-{\mathbf{u}}_p)(\mathbf{v}-{\mathbf{u}}_p)f_p$$

$${\mathbf{q}}_p=\frac{1}{2}\int d^{3}v(\mathbf{v}-{\mathbf{u}}_p{)}^2(\mathbf{v}-{\mathbf{u}}_p)f_p$$

$$\frac{\partial ⟨{\mathcal{E}}_p^f⟩}{\partial t}=⟨{\mathbf{j}}_p\cdot \mathbf{E}⟩+⟨\left({\mathbf{P}}_p\cdot \nabla \right)\cdot {\mathbf{u}}_p⟩$$

$$\frac{\partial ⟨{\mathcal{E}}_p^{th}⟩}{\partial t}=-⟨\left({\mathbf{P}}_p\cdot \nabla \right)\cdot {\mathbf{u}}_p⟩$$

$$\frac{\partial ⟨{\mathcal{E}}^m⟩}{\partial t}=-⟨\mathbf{j}\cdot \mathbf{E}⟩$$

$$\left({\mathbf{P}}_p\cdot \nabla \right)\cdot {\mathbf{u}}_p=P_p{\theta }_p+{\mathbf{\Pi }}_p:{\mathbf{D}}_p$$

$${ϵ}_p(\mathbf{x},t)=\frac{1}{n_p}\sqrt{\int (f_p-f_{{}_M}{)}^2{d}^{3}v}.$$

$$\Delta ⟨{\mathcal{E}}_p^{th}⟩(t)=-{\int }_0^{t}d\tau ⟨{\mathbf{\Pi }}_p:{\mathbf{D}}_p⟩(\tau )-{\int }_0^{t}d\tau ⟨P_p{\theta }_p⟩(\tau )$$

$$\Delta ⟨{\mathcal{E}}_p^{th}⟩(t)=\frac{3}{2}\left(n_{p,0}\Delta ⟨\delta T_p⟩+\Delta ⟨\delta n_p\delta T_p⟩\right)$$

$$\frac{3}{2}\left[\frac{\partial T_p}{\partial t}+\nabla \cdot \left({\mathbf{u}}_p{T}_p\right)\right]=-\frac{1}{n_p}\left(P_p{\theta }_p+{\mathbf{\Pi }}_p:{\mathbf{D}}_p+\nabla \cdot {\mathbf{q}}_p\right)$$

$$\frac{3}{2}{n}_{p,0}\frac{\partial ⟨T_p⟩}{\partial t}\simeq -⟨{\mathbf{\Pi }}_p:{\mathbf{D}}_p⟩,$$

Full text

I Introduction

The space physics community has become increasingly aware that plasma turbulence has significant impact on observed properties of solar wind and other accessible space plasmas Bruno and Carbone (2005); Matthaeus and Velli (2011). Given that the cascade rate towards smaller scales is measuredMacBride, Smith, and Forman (2008); Sorriso-Valvo et al. (2007); Hadid, Sahraoui, and Galtier (2017) and that an addition of proton internal energy is required in order to explain the observed temperature evolution within the inner heliosphereMatthaeus et al. (1999); Hellinger et al. (2013), it is natural to study the plasma dissipation that must be present to reconcile these observations. Accordingly, kinetic dissipation of turbulence in weakly collisional plasmas has become a very active topic throughout heliospheric physics and several standard approaches have been employed. One approach is to assume that the notion of Ohmic dissipation carries over from the collisional case, and therefore to adopt as a proxy for dissipation the electromagnetic work on particles, namely ${\bm{j}}\cdot{\bm{E}}$ ($\bm{j}$ the electric current density, $\bm{E}$ the electric field) (Sundkvist et al., 2007; Wan et al., 2015). Another definition, especially adopted in simulation works, regards the role of inter-particle collisions, that introduce irreversibility into the system (Howes et al., 2006; TenBarge and Howes, 2013; Pezzi, Valentini, and Veltri, 2016; Pezzi, 2017; Pezzi et al., 2019). Within this viewpoint, dissipation is associated with an entropy growth. Still different approaches study the dissipation related to peculiar mechanisms. Some examples are the dissipation of particular wave-modes, such as kinetic Alfvén(Salem et al., 2012; Chen, Klein, and Howes, 2019) and whistler waves(Chang, Gary, and Wang, 2015; Gary, Hughes, and Wang, 2016), and the heating related -directly or indirectly- to magnetic reconnection(Drake et al., 2009; Servidio et al., 2011, 2012; Osman et al., 2012a, b; Wu et al., 2013). Magnetic reconnection in particular play a decisive role in energizing particles and dissipating energy in turbulent plasmas(Cerri and Califano, 2017).

In each of these approaches there is clearly an opportunity for physical insights, even if each approach may not capture all mechanisms that operate, or their mutual interactions. The concept of dissipation here advanced –and implemented in the analysis of HVM simulation of plasma turbulence– differs in that it is based on the premise that pathways to dissipation may be identified based on separations of the effects of transport, scale transfer and energy conversion. In particular, one may uniquely characterize the conversion into internal energy, although this process is not related to an entropy growth.

Beginning with the HVM system, it is straightforward to derive energy balance relations (Krall and Trivelpiece, 1973), which may be further analyzed to identify the basic channels of energy exchange Yang et al. (2017a, b) and the evolution of the pressure tensorDel Sarto, Pegoraro, and Califano (2016); Del Sarto and Pegoraro (2018). One basic finding, reviewed in more detail below, is that the energy exchanged between the electromagnetic field and the proton bulk motions is due to the interaction of the electric field with the proton charge flux (current). Meanwhile, the flow kinetic energy may be exchanged with the internal energy through the combined action of the pressure-dilatation and the pressure-strain interactions.

In the Vlasov energy balance, one also finds several energy transport terms that –although being quite important in dynamics– do not contribute, on average, to converting energy from one form to another. In particular, transport terms cannot produce internal energy. These developments are based on elementary properties of the Vlasov equation, but their significance for collisionless turbulence has not always been fully recognized. Based on this taxonomy, one may proceed to analyze the contributions to these channels of energy conversion and transfer employing a number of approaches: by global averaging, by examining distribution functions, by computing joint distributions and correlations, and by studying real space maps to identify coherent structures Yang et al. (2017a, b). Some of these analyses have also been implemented in MMS observations in the terrestrial magnetosheath Chasapis et al. (2018). This program of analysis of energy conversion channels is undertaken here by using HVM turbulence simulation. It is notable that, for HVM simulation, it is possible to reconcile a total energy balance with good accuracy, an undertaking not easily achieved using methods such as kinetic PIC. One may also implement scale averaging(Germano, 1992) to examine locality of energy transfer and the relative contribution of conversion processes at different scales Yang et al. (2018); Camporeale et al. (2018). This will be undertaken for HVM in a subsequent study.

The paper is organized as follows. Section II introduces the theoretical model and the numerical setup of the numerical simulation performed here to implement the pressure-strain interaction Pi-D analysis, already adopted in recent studies Yang et al. (2017a, b, 2018); Chasapis et al. (2018). Then, in Sect. III, we provide a general overview of the simulation results, in terms of production of a turbulent scenario. In Sect. IV, we describe several outcomes concerning the application of the Pi-D analysis on our simulation. Finally, in Sect. V, we summarize the presented results.

II Theoretical model and channels of energy conversion

Within the HVM model (Valentini et al., 2007; Servidio et al., 2012; Perrone et al., 2013; Servidio et al., 2015; Valentini et al., 2016; Cerri, Servidio, and Califano, 2017; Perrone et al., 2018), protons are kinetic and electrons are a massless, isothermal fluid. Maxwell’s equations are further approximated by neglecting the displacement current and by assuming charge neutrality. The Faraday’s law and the generalized Ohm’s law are then coupled to the proton Vlasov equation. Dimensionless HVM equations are:

[TABLE]

where $f_{p}({\bm{x}},{\bm{v}},t)$ is the proton distribution function (DF); ${\bm{E}}({\bm{x}},t)$ and ${\bm{B}}({\bm{x}},t)$ are respectively the electric and magnetic fields; ${\bm{j}}=\nabla\times{\bm{B}}$ is the current density; and $n_{p}=\int d^{3}vf_{p}$ and ${\bm{u}}_{p}=\int d^{3}v{\bm{v}}f_{p}/n_{p}$ are respectively the proton number density and bulk speed. The electron pressure $P_{e}=n_{e}T_{e}$ ($k_{{}_{B}}=1$) depends on spatial coordinates through the electron density $n_{e}=n_{p}$, while the electron temperature is only a parameter $T_{e}=T_{p,0}$, since electrons are isothermal. Equations (1–3) have been normalized as follows. Time, velocities and lengths are respectively scaled to the inverse proton cyclotron frequency $\Omega_{cp}^{-1}=m_{p}c/eB_{0}$, to the Alfvén speed $c_{A}=B_{0}/\sqrt{4\pi n_{p,0}m_{p}}$, and to the proton skin depth $d_{p}=c_{A}/\Omega_{cp}$; being $m_{p}$, $e$, $c$, $B_{0}$ and $n_{p,0}$ respectively the proton mass and charge, the light speed, the background magnetic field and the equilibrium proton density.

The method described in the work of Yang et al. (2017a, b), briefly revisited in the following, is based on the energy conservation equations, which –for the HVM model– read as:

[TABLE]

where $\mathcal{E}^{f}_{p}=n_{p}{\bm{u}}_{p}^{2}/2$, $\mathcal{E}^{th}_{p}=\frac{1}{2}\int d^{3}vf_{p}\left({\bm{v}}-{\bm{u}}_{p}\right)^{2}=3n_{p}T_{p}/2$ and $\mathcal{E}^{m}={\bm{B}}^{2}/2$ are respectively the flow, internal, and magnetic energy densities. Note that the proton total energy is $\mathcal{E}_{p}=\frac{1}{2}\int d^{3}vf_{p}{\bm{v}}^{2}=\mathcal{E}^{f}_{p}+\mathcal{E}^{th}_{p}$. Moreover, the kinetic temperature $T_{p}$, implicitly defined as usual through the VDF $2^{nd}$-order moment, does not imply thermal equilibrium.

In Eqs. (4–6), the pressure tensor ${\bm{P}}_{p}$ and the heat flux ${\bm{q}}_{p}$ are defined as follows:

[TABLE]

To appreciate the conversions between different type of energy, Eqs. (4–6) are averaged on the spatial domain to get:

[TABLE]

where ${\bm{j}}_{p}=n_{p}{\bm{u}}_{p}$ ($e=1$ in scaled units) and all divergence terms vanish due to the spatial boundary conditions. Note that the sum of the terms at right-hand sides of Eqs. (9–11) is not zero, but reduces to the electromagnetic work performed on electrons ($\sim\langle{\bm{j}}_{e}\cdot{\bm{E}}\rangle$). This term, for the HVM algorithm, is determined by the hybrid plasma approximation described above, and is not reducible to a simple time derivative. Eqs. (9–11) suggest that the term $\langle{\bm{j}}_{p}\cdot{\bm{E}}\rangle$ controls the transfer between flow and magnetic energy, while $\langle\left({\bm{P}}_{p}\cdot\nabla\right)\cdot{\bm{u}}_{p}\rangle$ transforms flow energy into internal energy (Yang et al., 2018). This last term can be further decomposed as:

[TABLE]

where $P_{p,ij}=P_{p}\delta_{ij}+\Pi_{p,ij}$, $P_{p}=P_{p,ii}/3$, $\theta_{p}=\nabla\cdot{\bm{u}}_{p}$ and $D_{p,ij}=\left(\partial_{j}u_{p,i}+\partial_{i}u_{p,j}\right)/2-\theta_{p}\delta_{ij}/3$. In the last expression, $\delta_{ij}$ is the Kronecker’s delta and $\partial_{i}$ denotes the partial derivatives with respect to the $i$–th spatial coordinate. The two terms at the right-hand side of Eq. (12) can be respectively associated with plasma dilatations and compressions and with the trace-less pressure tensor (including also off-diagonal pressure terms). Hereafter, $P_{p}\theta_{p}$ and ${\bm{\Pi}}_{p}:{\bm{D}}_{p}$ are briefly denoted as the P-$\theta$ and Pi-D terms, respectively. Both terms appear with a negative sign in Eq. (10), thus implying that negative values of $\langle P_{p}\theta_{p}\rangle$ and $\langle{\bm{\Pi}}_{p}:{\bm{D}}_{p}\rangle$ have the effect of increasing the proton internal energy. Moreover, both terms are not positive definite, hence their signs are decisive to discriminate the direction of the energy transfer between flow and internal energy.

The idea adopted in recent studiesYang et al. (2017a, b, 2018) is that the Pi-D term can provide significant insights on the mechanisms that transfer energy towards smaller scales, where energy is eventually dissipated. Indeed, in a collisional plasma (MHD framework), the traceless pressure-tensor terms are usually related to viscous dissipation (Braginskii, 1965). Investigating whether this connection holds also in a weakly-collisional system can be hence useful to understand the nature of dissipation, although –within the HVM model– collisions are not modeled and the system is formally reversible(Pezzi, Valentini, and Veltri, 2016; Pezzi, 2017; Pezzi et al., 2019).

III Overview of the Numerical Simulations

Equations (1–3) are integrated in a 2.5D–3V phase space domain. We retain all velocity space directions, while physical-space vectors are three-dimensional but depend only on the two spatial coordinates $(x,y)$. This double-periodic spatial domain, whose size is $L_{x}=L_{y}=L=2\pi\times 20d_{p}$, is discretized with $N_{x}=N_{y}=512$ grid-points. The velocity domain is discretized with $N_{v_{x}}=N_{v_{y}}=N_{v_{z}}=71$ points in the range $v_{j}=\left[-5v_{th,p},5v_{th,p}\right]\;(j=x,y,z)$. Boundary conditions impose $f(v_{j}>5v_{th,p})=0$, where the proton thermal speed $v_{th,p}=\sqrt{k_{{}_{B}}T_{p,0}/m}$ is directly related to the Alfvén speed through $\beta_{p}=2v_{th,p}^{2}/c_{A}^{2}=2$.

The proton VDF is initially Maxwellian, with uniform unit density. A uniform background out-of-plane magnetic field ${\bm{B}_{0}}=B_{0}{\bm{e}}_{z}$ ($B_{0}=1$) is initially imposed. The equilibrium is perturbed by imposing an initial $2D$ spectrum of Fourier modes, equi-partitioned between proton bulk velocity and magnetic field, and having an amplitude such that $\delta B_{rms}/B_{0}=\delta u_{rms}/c_{A}=1/2$. The energy is injected with random phases and in the wave-number range $k\in[0.1,0.3]$, being $k=mk_{0}$ with $2\leq m\leq 6$ and $k_{0}=2\pi/L$. Neither density fluctuations nor parallel perturbations are introduced ($\delta n_{p}=\delta u_{p,z}=\delta B_{z}=0$). A small resistivity $\eta=1.5\times 10^{-2}$ is accurately introduced to avoid numerical instabilities.

Before focusing on the Pi-D analysis, we conclude this section by describing the simulation results concerning the onset of turbulence at proton inertial scales.

Owing to the presence of strong nonlinearities, the initial relatively large-scale perturbations couple and produce fluctuations at smaller scales. This cascade-like activity can easily be appreciated by examining the evolution of $\langle j_{z}^{2}\rangle$ [Fig. 1(a)], that represents a standard proxy for highlighting the presence of nonlinear couplings. The $\langle j_{z}^{2}\rangle$ peak, occurring at $t^{*}=30\Omega_{cp}^{-1}$, is associated with the most intense turbulent activity in the numerical simulation (Servidio et al., 2012). Note that the time $t^{*}$ corresponds to about $0.6\tau_{nl}$, where the nonlinear time is $\tau_{nl}=l_{c}/\delta u(l_{c})$, $l_{c}$ is the turbulence correlation length and $\delta u(l_{c})$ is the fluctuations amplitude at $l_{c}$. In our case, $l_{c}\simeq L/5$ and $\delta u(l_{c})\simeq\delta u_{rms}$. The generation of small-scale fluctuations can be recognized by looking at the omni-directional (perpendicular) power spectral densities (PSDs), reported in Fig. 1(b), as a function of $kd_{p}$. PSDs have been computed by averaging on the temporal range $t\Omega_{cp}\in[20,40]$ (shaded area in Fig. 1(a)), where $\langle j_{z}^{2}\rangle$ is nearly constant. This time interval corresponds to a quasi-stationary regime of turbulence. For about a decade of scales larger than the proton inertial scales, an inertial range is observed, wherein magnetic fluctuations are dominant and show a power-law-like PSD, which is reminiscent of the $-5/3$ Kolmogorov slope. At smaller scales, the electric field power increases as well as compressibility (Bale et al., 2005).

The real-space picture of out-of-plane electric current $j_{z}(x,y)$ at $t=t^{*}$, is reported in Fig. 1(c). This contour plot confirms an highly structured and inhomogeneous current density pattern, with strong current sheets and magnetic islands. In, or near, these strong current sheets, which themselves are embedded in a dynamic and strongly turbulent region, complex phenomena such as magnetic reconnection, energy dissipation and other kinetic effects (Servidio et al., 2015; Retinò et al., 2007; Osman et al., 2011; Karimabadi et al., 2013) are found. Intense current regions are often in close proximity to regions in which the plasma VDF is far from the equilibrium Maxwellian shape (Drake et al., 2003, 2010; Greco et al., 2012a; Swisdak, 2016; Pezzi et al., 2017). The presence of non-Maxwellian features in the proton VDF can be also associated with the enstrophy velocity-space cascade (Servidio et al., 2017; Cerri, Kunz, and Califano, 2018; Pezzi et al., 2018). To quantify deviations of the proton VDF from the associated Maxwellian $f_{M}$, we here adopt the $\epsilon_{p}$ parameter Greco et al. (2012a):

[TABLE]

where $f_{{}_{M}}$ has the same density, bulk speed and temperature of $f_{p}$. The temporal evolution of $\langle\epsilon_{p}\rangle(t)$, reported in Fig. 1(a), highlights that the generation of non-Maxwellian features grows together with the increase of the current activity (i.e. magnetic field gradients and small-scale fluctuations). The $\epsilon_{p}$ parameter finally saturates at a nearly constant value. The contour plot of $\epsilon_{p}(x,y)$ at $t=t^{*}$, reported in Fig. 1(d), demonstrates a qualitative correlation of $\epsilon_{p}(x,y)$ and $j_{z}(x,y)$.

IV The Pi-D analysis: energy conversion and temperature increase

IV.1 Bookkeeping the internal energy

In order to understand the energy transfer mechanisms occurring in a turbulent plasma, it is useful to apply the Pi-D method on the simulation described in Section II.

By integrating on time Eq. (10) and using Eq. (12), one gets:

[TABLE]

where $\Delta$ denotes the variation with respect to the initial values. Figure 2 reports the temporal evolution of the cumulative integrals of $\langle-{\bm{\Pi}}_{p}:{\bm{D}}_{p}\rangle$; $-\langle P_{p}\theta_{p}\rangle$ ; and $\langle-\left({\bm{P}}_{p}\cdot\nabla\right)\cdot{\bm{u}}_{p}\rangle$. These integrals have been numerically evaluated through the trapezoidal rule. The Pi-D term dominates over the pressure-dilatation term in the internal energy increase, suggesting that the cascade process remains quasi-incompressible. Since Eulerian simulations are almost noise-free, we can also report the evolution of the internal energy $\Delta\langle\mathcal{E}^{th}_{p}\rangle(t)=3\Delta\langle n_{p}T_{p}\rangle/2$. According to Eq. (10), the internal energy variation should be exactly equal to the cumulative integral of $\langle-\left({\bm{P}}_{p}\cdot\nabla\right)\cdot{\bm{u}}_{p}\rangle$. The slight difference between these two quantities is mainly due to the level of accuracy for the energy conservation recovered in the numerical simulation ($\sim 0.5\%$).

It can be also shown that the Pi-D term has a direct effect on the temperature growth. This effect can be highlighted by writing $n_{p}=n_{p,0}+\delta n_{p}$ and $T_{p}=T_{p,0}+\delta T_{p}$, where $\langle\delta n_{p}\rangle=0$ since the total mass is well preserved. Therefore:

[TABLE]

The first term in this expression, namely $3n_{p,0}\Delta\langle\delta T_{p}\rangle/2$, is clearly connected to temperature growth that uncorrelated from density fluctuations. It can be almost superposed on the Pi-D term evolution, the slight discrepancy arguably due to the simulation accuracy. This behavior can be easily explained by deriving, from Eq. (5), the evolution equation for the temperature $T_{p}$:

[TABLE]

By averaging on the spatial domain and by considering a weakly-compressible case as the one analyzed here (density fluctuations always remain quite weak being, at $t=t^{*}$, $\delta n_{p,rms}/n_{p,0}\sim 7.7\times 10^{-2}$), it can be easily deduced that the dominant term in Eq. (16) is the Pi-D one. Hence, one gets:

[TABLE]

that fully clarifies the behavior observed in Fig. 2.

On the other hand, the evolution of the second term in Eq. (15), namely $3\Delta\langle\delta n_{p}\delta T_{p}\rangle/2$, is similar to P-$\theta$, as it can be expected since eventual correlations between density and temperature can be understood in terms of compression and dilatation.

Note also that, in the current simulation, the average temperature increases without showing a saturation regime. The evolution of correlated density-temperature fluctuations instead saturates in the first stages of the numerical simulation. This is probably due to the finite level of compressive-like structures contained into the initial perturbations. This effect may possibly be understood in terms if nearly-incompressible MHD theory (Zank and Matthaeus, 1992). Detailed analysis of the effects of different initial conditions, namely by i) introducing a different level of initial compressibility and by ii) changing the value of $\beta_{p}$, that also regulates the level of compressive activity (Bale et al., 2009), will be the content of a future work.

IV.2 Role of ${\bm{j}}\cdot{\bm{E}}$

In the initial stage of the simulation and up to the peak of the turbulence ($t\lesssim t^{*}$) $\int_{0}^{t}dt\langle-{\bm{j}}\cdot{\bm{E}}\rangle$ is positive, associated with the initial increase of magnetic fluctuation energy at the expense of flow energy (a feature familiar from MHD simulations). This evolution does not compare well with the temperature growth in Fig. 2, suggesting that, despite being adopted frequently as a proxy for heating and dissipation, the electromagnetic work is not globally correlated with the temperature increase. This information is already partially contained in Eqs. (9–11), even though –in principle– there could be transfers between magnetic and internal energy, through the intermediary of the kinetic energy.

For the sake of completeness, in Fig. 2, we also report the cumulative integral of $\langle-{\bm{j}}_{p}\cdot{\bm{E}}\rangle$. Compared to the integral of total $\langle-{\bm{j}}\cdot{\bm{E}}\rangle$, the difference represents the physical work that is done on electrons. It is worth highlighting that here we are focusing on global (averaged) quantities, while ${\bm{j}}\cdot{\bm{E}}$ is locally important to identify sites of kinetic activity and heating occurs Zenitani et al. (2011); Wan et al. (2012); Greco et al. (2012b).

IV.3 Spatial Distributions and Correlations

Figure 3 reports the contour plots of $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ (a), $Q_{\bm{\omega}}={\bm{\omega}}^{2}/4\langle{\bm{\omega}}^{2}\rangle$ (b), $Q_{\bm{j}}={\bm{j}}^{2}/4\langle{\bm{j}}^{2}\rangle$ (c) and $Q_{D}={\bm{D}}_{p}:{\bm{D}}_{p}/4\langle{\bm{D}}_{p}:{\bm{D}}_{p}\rangle$ (d), at the time instant $t=t^{*}$. $Q_{\bm{\omega}}$, $Q_{\bm{j}}$, and $Q_{D}$ estimate the relative strength of rotation, current sheets, and strain, respectively.

Our finding supports the recent works Yang et al. (2017a, b, 2018) where a coarse-grained correlation between intense Pi-D regions as well as strong current structures and and high vorticity and velocity strain regions has been described. These correlations can be also appreciated by looking at Fig. 4, where the conditional averages of $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ with respect to $Q_{\bm{\omega}}$ (a), $Q_{\bm{j}}$ (b) and $Q_{D}$ (c) are reported: stronger Pi-D event regions are found near regions of strong fluid vorticity and strain, with a weaker correlation evident near regions of current enhancement. Although we are not modeling irreversible dissipation, the correlation between Pi-D and $Q_{D}$ –usually related to viscosity in collisional plasmas– suggests a connection between two intrinsically different physical systems (collisionless and collisional).

To further support the idea that the plasma heating is an inhomogeneous process (Osman et al., 2011), Fig. 5 shows the probability density function (P.D.F.) of $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$. The Pi-D P.D.F. is far different from the Gaussian expectation, confirming the presence of intermittency. The Pi-D P.D.F. is also skewed, or slightly elongated in the positive direction of $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$, i.e. in the direction of net increase of the internal energy, as indicated by the positive sign of the $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ average.

Since Eulerian Vlasov simulations are almost noise-free, we can reveal whether the Pi-D term is correlated with kinetic characteristics of the proton VDF. This aspect can be also addressed with PIC simulations, but only when a quite large number of particles per cell are adopted (Camporeale and Burgess, 2011; Franci et al., 2015, 2018); otherwise several of the features shown below are masked due to the presence of numerical noise (Pezzi et al., 2017). Figure 6 reports the joint P.D.F. of the total temperature $T_{p}$ (a), temperature anisotropy $1-T_{p,\perp}/T_{p,\parallel}$ (b) and non-Maxwellianity parameter $\epsilon_{p}$ (c) with $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$. Each joint P.D.F. reveals an asymmetry in the positive direction $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$. A subtle correlation between regions of strong Pi-D and areas where plasma is locally warmer and non-Maxwellian can be seen. However, the Pi-D term is broadly distributed and therefore not strongly correlated with proton temperature heating nor with the plasma non-Maxwellianity, in part due to transport effects. Indeed, strong $\epsilon_{p}$, $T_{p}$ and $1-T_{p,\perp}/T_{p,\parallel}$ are also associated with very weak $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ regions.

The Pi-D term reveals hence regions where kinetic processes may be at work, but the viceversa is not true. This is probably caused by the fact that other terms can be locally important and may break the correlation. With regard to this last comment on correlations between Pi-D and kinetic signatures, it should be borne in mind that, when examining the Pi-D term at the local level or using in-situ measurements(Chasapis et al., 2018), transport terms may also play a role. Some of these terms, that are related to divergence, move energy from one location to another without globally changing the form of energy.

Furthermore, both the magnetic field and the spatial derivatives of the bulk speed enter into the dynamics of the pressure tensor Del Sarto, Pegoraro, and Califano (2016). Consequently, various kinetic and fluid signatures of dissipation, such as Pi-D, current, ${\bm{j}}\cdot{\bm{E}}$, vorticity, temperature anisotropy and non-Maxwellianity, may be coarsely or regionally correlated, but pointwise offset from one another. This type of regional clustering of quantities related to dissipation has been observed in several simulation studies Wan et al. (2012); Servidio et al. (2015); Yang et al. (2018); Camporeale et al. (2018) and is an important facet of intermittency in plasma turbulence.

IV.4 Role of local transport terms

To further support that local effects are important in the energy transfer towards smaller scales, where dissipation is thought to occur, we conclude this section by selecting two vertical cuts in the spatial domain of the numerical simulation that cross strong current sheet.

Figure 7 shows these one-dimensional cuts, taken at $x=41.7d_{p}$ (left column) and $y=122.7d_{p}$ (right column). Each row reports, from top to bottom, (a) magnetic field components $B_{x}$, $B_{y}$ and $B_{z}$; (b) current density ${\bm{j}}^{2}$; (c) proton temperature $T_{p}$; (d) temperature anisotropy $1-T_{p,\perp}/T_{p,\parallel}$; (e) non-Maxwellian measure $\epsilon_{p}$; (f) ${\bm{j}}\cdot{\bm{E}}$; (g) the heat flux divergence $\nabla\cdot{\bm{q}}_{p}$; (h) $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$; and (i) $-P_{p}\theta_{p}$.

In both events, at least one component of the magnetic field changes its sign and the current density has a significant growth, characterized by a double-peaked structure in the left column, as recently observed also in Earth’s magnetosheath (Chasapis et al., 2018). The proton temperature also increases close to the current density peaks. In the left-column panel, plasma is hotter inside the current sheet, while in the right-column panel, plasma heating is slightly displaced from the current sheet Servidio et al. (2015).

Both proxies of non-Maxwellianity reported in Fig. 7, namely $1-T_{p,\perp}/T_{p,\parallel}$ (d) and $\epsilon_{p}$ (e), present a broad structure where both parameters are significantly non-null. The presence of such large-scale structures of non-Maxwellianity also causes the lack of correlation between the non-Maxwellianity parameters and Pi-D. Super-imposed to these large-scale structures, quite strong peaks of $1-T_{p,\perp}/T_{p,\parallel}$, co-located with the temperature increase and associated with a preferential perpendicular heating ($T_{p,\perp}>T_{p,\parallel}$), are also observed.

The ${\bm{j}}\cdot{\bm{E}}$ plot indicates that, although this term is not globally connected with temperature increase, it can be locally significant, close to current sheets and regions where kinetic effects are important and intermittent dissipation may be active (Osman et al., 2011).

The heat-flux divergence panels (g) support the idea that this term can be locally important, although it does not globally change the proton internal energy. Indeed, in the left-column panel, a strong peak close to current sheets is recovered, while in the right-column panel, $\nabla\cdot{\bm{q}}_{p}$ is almost flat.

Both $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ (g) and $-P_{p}\theta_{p}$ (h) report variations across the structures, ranging from positive to negative values. Focusing on the Pi-D term, it is apparent in the left-column panel, that the positive peak of Pi-D is displaced from the current structure and it is not associated with the temperature growth. On the other hand, in the right-column panel, two weak negative $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ peaks are found in correspondence of the current structure, while $P_{p}\theta_{p}$ slightly increases.

Figure 7 helps to realize the complexity of truly understanding how energy is dissipated in a turbulent, weakly-collisional plasma. It is largely established that temperature growth occurs in or around current sheets. The same holds for temperature anisotropy and $\epsilon_{p}$, although structures are wider. The transfer mechanism based on the Pi-D term is somehow correlated with hotter plasma regions, although other terms ($P_{p}\theta_{p}$, heat-flux) may locally play a role. In our opinion, this should motivate extremely accurate in-situ measurements to better appreciate the role of these variables.

V Conclusion

We have examined the production of internal energy in a collisionless plasma employing an Eulerian HVM scheme and the pressure strain formalism Yang et al. (2017a, b, 2018).

To our knowledge, a study of this type has not been reported using an Eulerian approach. Relative to previous reports based on PIC models, the HVM approach has the advantage that the proton velocity distribution is well resolved in three dimensions and for several thermal speeds. Therefore we may have some confidence that the subtle interactions between the pressure tensor and the spatial derivative of fluid velocity entailed in the pressure-strain interaction are well-represented in the computed solution. The HVM method is much less “noisy” than PIC, especially with regard to particle counting noise associated with finite particle number (see e.g., discussion in Pezzi et al. (2017); Haggerty et al. (2017)). This permits HVM to compute accurate structures in velocity space as required for the present study.

Enabled in this way, we have examined characteristic kinetic properties during turbulence evolution, near the time of maximum mean square current density, when nonlinear effects are strongest. Indeed we have observed at this time the emergence of strong coherent current structures in real space, representing intermittency, and strongly non-Maxwellian features in the proton VDFs, as has been previously reported Greco et al. (2012b); Servidio et al. (2015).

One of the very basic demonstrations of the efficacy of the HVM simulation strategy, and of the physical relevance of the pressure-strain relations, is seen in Fig 2. Here we have shown that the time integration of the sum of the two volume-integrated pressure-strain effects –the pressure-dilatation P-$\theta$ and the Pi-D terms– provides a very accurate account of the change of internal energy due to the turbulence cascade. On the other hand, the electromagnetic work on particles, ${\bm{j}}\cdot{\bm{E}}$, does not correlate well with global internal energy increases. This would be anticipated based on the analysis leading to Eqs. (9-11), as one observes that the charge flux of species interacts with the electric field only to change the species flow energy at the expense of magnetic energy.

The present study has also reported that the values of the pressure-strain energy conversion ($-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$) are broadly distributed and skewed towards positive values associated with production of internal energy. Spatial maps furthermore indicate that strong $-{\bm{\Pi}}_{p}:{\bm{D}}_{p}$ values appear in sheet (or core) like patterns indicating intermittent conversion of flow energy into internal energy. Both of these distributions are typical of quantities related to cascade and pathways to dissipation, such as the“local energy transfer” (LET) Yang et al. (2018); Sorriso-Valvo et al. (2018, 2019) or ${\bm{j}}\cdot{\bm{E}}$ in various reference frames Wan et al. (2012, 2015).

It remains to draw some conclusion regarding the idea of “dissipation”. A reasonable definition of dissipation may require the conversion of energy into internal energy. Leaving momentarily aside the role of collisions in energy conversion and introduction of irreversibility, the spatial concentration of Pi-D observed here (and earlier in PIC simulations (Yang et al., 2017b, a)) appears to be precisely what is required for a plasma physics adaptation of the Kolmogorov Refined Similarity Hypothesis. In such a development, the concentration of dissipation is directly related to scaling of velocity increments, including multi-fractal behavior, etc. In this way, it appears likely that the pressure-strain interaction will enter prominently into more formal and more complete theories of plasma turbulence that remain to be developed.

The role of collisions in attenuating velocity-space structures and in securing irreversibility remains to be further explored (Pezzi, Valentini, and Veltri, 2016; Pezzi, 2017; Pezzi et al., 2019). This need is particularly clear pertaining to the relationship of weak collisions to the strength of pressure-strain interactions, since formal (viscous) closure relating the two is not available in the weakly collisional case.

Furthermore, within a collisionless system, the internal energy growth is not uniquely related to the increase of the random particle energy, since the internal energy may also increase for the presence of non-Maxwellian features such as accelerated particle beams. The VDF free-energy contained in such structures(Lesur and Diamond, 2013) can be backward released in form of ordered energy, such as micro-instabilities without increasing the plasma entropy(Hellinger et al., 2017). Collisions would intrinsically inhibit the possibility to transfer back this from particles to field this energy, since –dissipating these structures– they incessantly push plasma toward thermal equilibrium (Pezzi, Valentini, and Veltri, 2016; Pezzi, 2017). This further level of complexity should motivate subsequent studies to provide insights on the strength of the physical ingredient that ultimately dissipates energy in a irreversible way(Pezzi et al., 2019).

Acknowledgements.

This research is supported in part by NASA through the MMS Theory and Modeling Program NNX14AC39G and the Parker Solar Probe project under Princeton subcontract SUB0000165, and HSR grant 80SSC18K1648. This paper has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 776262 (AIDA, www.aida-space.eu). OP and SS are partly supported by the International Space Science Institute (ISSI) in the framework of the International Team 405 entitled “Current Sheets, Turbulence, Structures and Particle Acceleration in the Heliosphere”. Numerical simulations discussed here have been performed on the Marconi cluster at CINECA (Italy), within the projects IsC53_RoC-SWT and IsC63_RoC-SWTB, and on the Newton cluster at the University of Calabria (Italy).