Turbulent dissipation challenge: a community-driven effort

Abstract

Many naturally occurring and man-made plasmas are collisionless and turbulent. It is not yet well understood how the energy in fields and fluid motions is transferred into the thermal degrees of freedom of constituent particles in such systems. The debate at present primarily concerns proton heating. Multiple possible heating mechanisms have been proposed over the past few decades, including cyclotron damping, Landau damping, heating at intermittent structures and stochastic heating. Recently, a community-driven effort was proposed (Parashar & Salem, 2013, arXiv:1303.0204) to bring the community together and understand the relative contributions of these processes under given conditions. In this paper, we propose the first step of this challenge: a set of problems and diagnostics for benchmarking and comparing different types of 2.5D simulations. These comparisons will provide insights into the strengths and limitations of different types of numerical simulations and will help guide subsequent stages of the challenge.

1. Motivation

Near-Earth space is a dynamic system where plasmas from the Sun and Earth’s magnetosphere interact. The dynamics of this system affect not only our spacecraft and astronauts but also wireless communications and, in some extreme cases, electrical grids on Earth. It is very important to have good predictive capabilities for this system. An integral part of the system is the medium that connects the Sun and the Earth: the solar wind.

The solar wind is observed to be turbulent. By this we mean that the energy in the electromagnetic field and plasma fluctuations at the largest scales is ‘cascaded’ to smaller scales via nonlinear interactions, creating a broadband spectrum of incoherent fluctuations in fields and flows. The distribution of energy over a range of scales usually follows a power law (e.g. Kolmogorov 1941). The solar wind magnetic energy spectrum obeys the Kolmogorov power law in the ‘inertial range’ (e.g. Coleman 1968; Matthaeus & Goldstein 1982; Goldstein, Roberts & Matthaeus 1995), which covers scales larger than a break scale ${\it\lambda}_{b}$ that is comparable to the proton inertial length $d_{i}$ and/or proton gyro-radius ${\it\rho}_{i}$ , where $d_{i}=c/{\it\omega}_{pi}$ , $c$ is the speed of light and ${\it\omega}_{pi}$ is the proton plasma frequency. (Near Earth, $d_{i}$ is comparable to ${\it\rho}_{i}$ .) Near observed frequencies of about 0.5 Hz, magnetic fluctuation spectra in the solar wind usually show a break from the well-known $k^{-5/3}$ form of the inertial range to steeper spectra of the form $k^{-{\it\alpha}}$ with ${\it\alpha}>2$ (e.g. Sahraoui et al. 2009; Alexandrova et al. 2012). Under the Taylor hypothesis that the observed frequency of solar wind fluctuations should scale as $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{sw}/2pi$ , where $\boldsymbol{v}_{sw}$ is the solar wind flow velocity, this spectral break corresponds to $kd_{i}\simeq 1$ . Although it is generally agreed that inertial range turbulence is well described by magnetohydrodynamic (MHD) models, it is also generally agreed that the steeper spectra at wavelengths shorter than the spectral break require a kinetic, that is, velocity-space, description for their complete representation. This short-wavelength, relatively steep spectral regime is often labelled the ‘dissipation range’, suggesting that it plays an important role in turbulent dissipation and ion heating. But the steepening of this regime may also be attributed to fluctuation dispersion introduced by velocity-space effects and electron dynamics (e.g. Ghosh & Goldstein 1997; Stawicki, Gary & Li 2001; Boldyrev & Perez 2012), so that we label this the ‘kinetic range’.

The relative roles of dissipation and dispersion are poorly understood and is an important topic of study in plasma turbulence research. Dissipation (that is, the transfer of energy from the field fluctuations to the plasma particles) by fluctuations at wavelengths near the spectral break is primarily collisionless and must explain features observed in the solar wind, such as the strong anisotropies of the ion velocity distributions (e.g. Marsch 2006) which often correspond to temperatures such that $T_{\bot }>T_{\Vert }$ . By ‘temperature’ in this context, we mean the second moment of the velocity distribution function as measured by spacecraft observations. The anisotropy is introduced by the collisionless dynamics and eventually thermalized by collisions. A number of damping mechanisms have been proposed, including ion–cyclotron damping (Hollweg & Isenberg 2002, and references therein), Landau damping (Schekochihin et al. 2009, and references therein), particle energization at current sheets and reconnection sites (e.g. Dmitruk, Matthaeus & Seenu 2004; Sundkvist et al. 2007; Parashar et al. 2009, 2011; Osman et al. 2011; Wan et al. 2012; Karimabadi et al. 2013) and stochastic heating (e.g. McChesney, Stern & Bellan 1987; Chaston et al. 2004; Chandran et al. 2010; Bourouaine & Chandran 2013; Xia et al. 2013). If the power in electromagnetic fluctuations at scales ${\sim}{\it\rho}_{i}$ becomes sufficiently large, the ion’s orbit becomes chaotic and it undergoes ‘stochastic’ heating in the perpendicular direction. Ion–cyclotron damping requires the presence of waves with frequencies comparable to the proton cyclotron frequency ${\it\Omega}_{ci}$ , while particle energization at current sheets requires the presence of coherent structures. On the other hand, models of Landau damping (TenBarge et al. 2013) and stochastic heating (Chandran et al. 2010) have been applied to both linear waves and coherent structures.

There are problems associated with each of these mechanisms that make it difficult to propose any of them as the dominant mechanism for heating solar wind protons. Cyclotron heating requires sufficient energy at $k_{\Vert }d_{i}\sim 1$ but plasma turbulence with a mean magnetic field (like the solar wind) becomes highly anisotropic with most of the energy residing in high $k_{\bot }$ . There is still no consensus on how to provide energy to high enough $k_{\Vert }$ for this mechanism to be valid. Landau damping inherently produces high $T_{\Vert }$ and hence would provide parallel proton heating instead of the observed perpendicular heating. Intermittent structures have been shown to have a strong correlation with highly anisotropic heating but the cause–effect relation is not completely understood and the exact details of how the heating takes place are debated.

While we focus on the solar wind as an example, there are many astrophysical plasmas for which these processes are important. A few examples are: the early Universe, the intergalactic medium (e.g. Schekochihin et al. 2010), accretion disks (e.g. Kunz, Schekochihin & Stone 2014), astrophysical jets (e.g. Romanova & Lovelace 1992), the interstellar medium, planetary magnetospheres, geologically active moons within planetary magnetospheres or magnetic moons within planetary magnetospheres (e.g. Saur et al. 2002; Mousis & Gautier 2004; Saur 2004; Kivelson 2006; von Papen, Saur & Alexandrova 2014).

We anticipate that most of these processes should be active in the solar wind at any given moment. The fundamental problem of interest is to quantify their contributions under various solar wind conditions. Determining the relative contributions of these mechanisms to turbulent dissipation in the solar wind is the central goal of the ‘turbulent dissipation challenge’. We envision this challenge as consisting of several steps. The overall goal of the ‘turbulent dissipation challenge’ is stated as follows:

What is the minimum physical model required to produce accurate intermittency, fluctuation dissipation and particle heating for collisionless plasma turbulence?

With this in mind, we propose a suggestion for the first steps that should be taken to answer this question: a comparison between different numerical codes simulating the same initial conditions and plasma parameters. The goal of this first step is to gain insights into the strengths and limitations of different simulation techniques, which will help us to design the future stages of this challenge. The eventual goal of the challenge is to run large three-dimensional simulations, but the high cost of such work limits the number of attempts that can be made, and so the challenge will initially focus on the more achievable goal of comparing 2.5D (two-dimensional (2D) grid and three-dimensional (3D) vectors) simulations with the aim of performing only the most relevant simulations in three dimensions.

We describe our overall goals in more detail in § 2, and the numerical set-up for the 2.5D challenge problems in § 3. We outline a set of common diagnostics in § 4. We suggest additional tasks related to observations in § 4.3, and in § 5 we outline future directions for this challenge.

2. Definition of goals

An underlying issue that has hindered progress towards a more complete understanding of dissipation in the solar wind is a lack of comparison between results from various studies. Most of the numerical studies are done using models that are vastly different in their underlying assumptions and numerical schemes. For example, electron magnetohydrodynamics (EMHD) (e.g. Kingsep, Chukbar & Ian’kov 1987; Gordeev, Kingsep & Rudakov 1994; Cho & Lazarian 2004) is the fluid description of the electron dynamics that treats protons as an immobile neutralizing background. Hybrid particle-in-cell (hybrid PIC) simulations (e.g. Parashar et al. 2009) treat protons as particles and electrons as a neutralizing fluid (usually massless and isothermal). Gyro-kinetics (e.g. Brizard & Hahm 2007) is a Vlasov–Maxwell system in which the gyro-motion of the particles has been averaged out from the system. Other examples of possible models include magnetohydrodynamics (MHD) (e.g. Gurnett & Bhattacharjee 2005), Hall MHD (e.g. Ghosh & Goldstein 1997), Hall-FLR MHD (e.g. Mjølhus 2009) with finite Larmor radius (FLR) corrections, reduced MHD (RMHD) (e.g. Zank & Matthaeus 1992), electron reduced MHD (ERMHD) (e.g. Schekochihin et al. 2009), Landau fluid (e.g. Goswami, Passot & Sulem 2005; Hunana et al. 2013), hybrid Eulerian Vlasov (e.g. Valentini, Califano & Veltri 2010) and full PIC (e.g. Karimabadi et al. 2013).

Complicating matters further, studies performed using these different models do not have the same initial conditions or even the same parameters. Hence it is not possible to compare the findings of these studies.

Similar problems concerning the variety of methods used and difficulty of benchmarking these methods are present with spacecraft data analysis. Data can be chosen from different instruments on many spacecraft, e.g. WIND (Harten & Clark 1995), CLUSTER (Escoubet, Fehringer & Goldstein 2001), ACE (Stone et al. 1998), Helios (Rosenbauer et al. 1976), Ulysses (Wenzel et al. 1992) and Voyager (Behannon et al. 1977), to name a few. Even if the data are from the same instruments and spacecraft, then the intervals chosen could be vastly different, ranging from fast wind to slow wind to all-inclusive. Moreover, the analysis techniques have not been methodically benchmarked against one another or simulations.

The above-mentioned problems call for a systematic comparison study where different types of simulation models are used to study the same initial conditions under as similar physical and numerical parameter regimes as possible. Artificial spacecraft data from these simulations can then be used to benchmark different spacecraft data analysis techniques with simulations as well as with each other.

The first step of the endeavour is to develop a better understanding of the strengths and limitations of different techniques for simulating intermittent structures and wave physics. Some work has been done in this direction (e.g. Henri et al. 2013; Keppens et al. 2013), but a comparison of statistical features of turbulence, dynamics of kinetic Alfvénic fluctuations as well as plasma heating have to be addressed. Hence we define the first step of this endeavour to be the following:

With given initial conditions and physical parameters, how do different simulation models of turbulence compare in capturing the physics of (i) intermittent structures and (ii) wave physics, in a turbulent setting?

In the sections to follow, we describe a plan for carrying out such comparisons. We focus on two specific problems with fixed plasma parameters and initial conditions.

1. (a) The first problem will be designed to generate strong turbulence with intermittent structures. It will have 3D vectors defined on a 2D computational grid that is almost perpendicular to the mean magnetic field to emphasize the role of intermittent structures.

2. (b) The second problem will be designed to develop a turbulent spectrum of waves. It will have 3D vectors on a 2D grid that is in the plane containing the mean magnetic field to emphasize wave processes like such as Landau and cyclotron damping.

Towards the end of the first stage, depending on the availability of time and resources, preliminary 3D simulations with minimal departures from the above-mentioned initial conditions could be performed. This will provide important clues about the importance of out-of-plane couplings missing in the 2.5D description.

By performing these simulations using different simulation models and performing the same set of diagnostics, we will be able to quantitatively compare the results from different models and gain insights into their relative strengths and weaknesses. Artificial spacecraft data will be provided to observers in order to establish a system that will be used to address physics questions at a later stage in the challenge.

3. Problem description

The conditions for the first simulations will be chosen to represent the conditions of ambient fast solar wind at 1 AU. Fast wind fills the majority of the heliosphere and is the cause of geo-effective magnetic sub-storms in the Earth’s magnetosphere. The fast wind also generally contains larger-amplitude but more homogeneous fluctuations than the slow wind, making it ideal for studies of plasma turbulence. Hence we will set ${\it\beta}_{i}\equiv 8{\rm\pi}k_{B}n_{0}T_{i}/\boldsymbol{B}_{0}^{2}=0.6$ along with $T_{e}=T_{i}$ , where $n_{0}$ is the initial proton density, and $T_{i}$ and $T_{e}$ are proton and electron temperatures. We will set $m_{e}/m_{i}=0.01$ and ${\it\omega}_{pe}/{\it\omega}_{ce}=1.5$ , where $m_{e}$ and $m_{i}$ are electron and proton masses, ${\it\omega}_{pe}$ is the electron plasma frequency and ${\it\omega}_{ce}$ is the electron cyclotron frequency. These latter two conditions imply that $V_{A}/c=1/15$ where $V_{A}=B_{0}/\sqrt{4{\rm\pi}n_{0}m_{i}}$ is the Alfvén speed. These parameters give us a Debye length of ${\it\lambda}_{D}=0.05d_{i}$ . It should be noted that these artificial values are used to reduce the computational costs for fully kinetic plasma codes. At 1 AU, ${\it\omega}_{pe}/{\it\omega}_{ce}\sim 113$ and $V_{A}/c\sim 2\times 10^{-4}$ (assuming typical fast wind values of $B_{0}\sim 11~\text{nT}$ and $n\sim 15~\text{cm}^{-3}$ ). These realistic numbers along with the realistic mass ratio $m_{e}/m_{i}=1/1836$ would make the computational cost of the proposed systems exorbitantly high. Our suggested parameters correspond to an electron thermal speed of order $c/2$ , implying that relativistic electron effects may play a role in our simulations. In the case of relativistic codes, varying ${\it\omega}_{pe}/{\it\omega}_{ce}$ does not affect the spectral features significantly. However, this ratio will be important in comparing results between relativistic and non-relativistic codes. This difference could provide guidance on the interpretation of discrepancies between the two codes.

To further reduce computational costs, we will perform the simulations in 2.5D, i.e. the simulation dynamics will be in a plane but with 3D components for all the vectors. By working in 2.5D, we will be able to simulate a much larger range of scales than would be possible in 3D simulations. Observations show that the anisotropy and intermittency grow with the width of the inertial range (e.g. Greco et al. 2008; Wicks et al. 2010, 2012; Wu et al. 2013). Hence, we choose to work with systems with a fairly large Reynolds number ² to compare the codes, which is possible in 2.5D (at least for the simulations described in 3.1). We understand that full 3D simulations are required for a complete description of dynamics, but we postpone 3D simulations until a later stage of the challenge. After the code comparison, large 3D simulations can be designed using the insights gathered from the first stage of the challenge.

3.1. Case 1: plane nearly perpendicular to $\boldsymbol{B}_{0}$

These simulations will study the intermittent structures that emerge from the nonlinear development of a Kelvin–Helmholtz (KH) instability (e.g. Chandrasekhar 1961; Miura & Pritchett 1982; Henri et al. 2013; Karimabadi et al. 2013). The KH instability gives rise to large-scale vortices and current sheets. As the vortices roll up, the current sheets get thinner and give rise to secondary tearing instabilities. This generates a turbulent ‘soup’ of current sheets ranging in scales from proton to electron scales. The ease of setting up KH and the broad range of turbulent current layers generated by it make it an ideal candidate for studying how well different codes capture intermittent physics. In the linear limit, it has been used to test PIC codes by comparing the linear damping rates to the ones captured in PIC simulations (Haugbølle, Frederiksen & Nordlund 2013).

We will follow the set-up used in Karimabadi et al. (2013) with a slight modification for this test. The initial density $n_{0}$ and magnetic field $\boldsymbol{B}_{0}$ will be uniform. The field $\boldsymbol{B}_{0}$ will be predominantly out of the simulation ( $xy$ ) plane, with an inclination such that $\boldsymbol{B}=B_{0}[\boldsymbol{e}_{y}\sin {\it\theta}+\boldsymbol{e}_{z}\cos {\it\theta}]$ with ${\it\theta}=2.86^{\circ }$ . To allow for use of periodic boundary conditions, we will use a double shear layer of flow. The shears will be present at $x_{1}=0.25L_{x}$ and $x_{2}=0.75L_{x}$ , with the shear layer flow defined as $\boldsymbol{U}=U_{0}(\text{tanh}(x-x_{1})/L_{V}-\text{tanh}(x-x_{2})/L_{V}-1)\boldsymbol{e}_{y}$ , where $U_{0}=10V_{A}^{\ast }$ , $V_{A}^{\ast }=B_{0}\sin ({\it\theta})/\sqrt{4{\rm\pi}n_{0}m_{i}}$ and $L_{V}=4d_{i}$ . The system size will be $L_{x}=125d_{i}$ , $L_{y}=125d_{i}$ , $N_{x}=2048$ and $N_{y}=2048$ . Also a perturbation of the form ${\it\delta}\boldsymbol{U}={\it\delta}U_{0}\sin (2k_{0}y)\exp (-(x-x_{i})^{2}/L_{V}^{2})\boldsymbol{e}_{y}$ , where $x_{i}=x_{1},x_{2}$ ,  $k_{0}=2{\rm\pi}/L_{y}$ and ${\it\delta}U_{0}=0.15U_{0}$ , will be added in both shear layers to expedite the growth of the instability. For PIC codes, particles will be loaded so as to create a net charge in the transition layer consistent with Gauss’s law and the electric field associated with the cross-field flow, $\boldsymbol{E}=-\boldsymbol{U}\times \boldsymbol{B}_{0}/c$ , which has non-zero divergence.

Fluid codes need to use some kind of resistivity and/or viscosity. In MHD codes, this can be set such that the Kolmogorov scale appears at the wavenumber that is equivalent to proton scales in the fully kinetic simulations. Hence viscosity should be set up to be significant at $k\sim 125(2{\rm\pi}/L_{box})$ .

3.2. Case 2: 2.5D plane containing $\boldsymbol{B}_{0}$

These simulations will address the temporal development of wave turbulence in a collisionless plasma. The specific technique will use several different simulation models to compute the evolution of initial field fluctuation spectra that are relatively long-wavelength and relatively isotropic in wavevector distribution. Recent observations and simulations have demonstrated that the fluctuations of kinetic range turbulence in the solar wind may consist of various normal modes including kinetic Alfvén waves (KAWs), magnetosonic waves (e.g. Svidzinski et al. 2009) and whistler waves (e.g. Saito et al. 2008; Gary, Chang & Wang 2012; Chang, Gary & Wang 2013, 2014). It is generally agreed that, in solar wind at 1 AU, turbulence near and at wavelengths shorter than $kd_{i}\sim 1$ consist primarily of KAWs (e.g. Howes et al. 2008; Sahraoui et al. 2009; Shaikh & Zank 2009; Podesta, Borovsky & Gary 2010; Salem et al. 2012; Chen et al. 2013).

The procedure to be followed in this second set of simulations will follow that used in several previous simulation studies of the kinetic range of solar wind turbulence. The computations are initialized with a narrowband spectrum of relatively isotropic, long-wavelength normal modes; nonlinear processes lead to a forward cascade that develops a broadband, anisotropic spectrum of turbulence. This technique has been used with PIC simulations to study 2.5D whistler turbulence (e.g. Gary, Saito & Li 2008; Saito et al. 2008; Saito, Gary & Narita 2010; Saito & Gary 2012), fully 3D whistler turbulence (e.g. Gary et al. 2012; Chang, Gary & Wang 2011, 2013, 2014, 2015) and 2.5D magnetosonic turbulence (e.g. Svidzinski et al. 2009). The same technique has also been used with hybrid simulations to examine the forward cascade of magnetosonic turbulence in 2.5D (e.g. Markovskii & Vasquez 2010; Markovskii, Vasquez & Chandran 2010) and Alfvénic and kinetic Alfvénic turbulence in full 3D (e.g. Vasquez, Markovskii & Chandran 2014) geometries. Tronko et al. (2013) and Howes (2015) have argued that, for incompressible MHD, the nonlinear cascade is suppressed in 2.5D geometry with the mean field in the plane of simulation. The implications of their conclusions are even stronger in the case of kinetic Alfvénic fluctuations and reduced models like RMHD and ERMHD. It is however unclear how the conclusions would carry over to compressible and/or kinetic limits. Moreover, computations in the same 2.5D geometry using hybrid PIC simulations of Alfvénic turbulence (e.g. Verscharen et al. 2012; Comişel et al. 2013; Comişel, Narita & Motschmann 2014) have shown that the cascade proceeds to kinetic scales anisotropically. Since the turbulent dissipation challenge addresses dissipative processes at proton scales, we will use an initial condition consisting of kinetic Alfvén waves in the second set of simulations.

In designing these simulations, we are guided by the following considerations.

1. (a) As the challenge is focused on the proton inertial scales, the initial condition should be a spectrum of KAWs around $kd_{i}\sim 1$ .

2. (b) The spectral range simulated should be the same over all the simulations, with in-plane or out-of-plane $B_{0}$ : approximately a decade above $d_{i}$ and at least a decade below  $d_{i}$ .

3. (c) The problem should run in reasonable amount of time on modern computing clusters. With $k_{\bot }\gg k_{\Vert }$ , we can reduce the perpendicular extent of the box to save computational time.

Given the above constraints, we set $L_{x}=125d_{i}$ , $L_{y}=31.25d_{i}$ , $N_{x}=2048$ and $N_{y}=256$ . This will give us ${\rm\Delta}x=0.061d_{i}$ corresponding to about 16 grid points across $d_{i}$ and three grid points across $2d_{e}$ . This will also correspond to the spectral range extending from $k_{min}\equiv 2{\rm\pi}/L_{x}=0.05d_{i}^{-1}$ to $k_{max}\equiv {\rm\pi}/{\rm\Delta}x=51.5d_{i}^{-1}$ .

Figure 1 shows the dispersion curves for various angles of propagation of KAWs in a Maxwellian proton–electron plasma, calculated using linear Vlasov theory (e.g. Gary 1986, 2005), with ${\it\beta}=0.6$ , $T_{e}=T_{i}$ , $m_{e}/m_{i}=0.01$ and $V_{A}/c=1/15$ . Almost all of these curves are dispersive at $kd_{i}\sim 1$ and have very low and comparable damping rates. We choose the window $0.4\leqslant kd_{i}\leqslant 1.85$ , which will give us not only dispersive waves and similar damping rates but also ${\sim}28$ grid points across the smallest wavelength in the initial spectrum. This means that the shorter wavelengths generated by the cascade will have sufficient resolution across them.

Figure 1. Dispersion curves for KAWs propagating at various angles in a proton–electron Maxwellian plasma with ${\it\beta}=0.6$ , $T_{e}=T_{i}$ , $m_{e}/m_{i}=0.01$ and $V_{A}/c=1/15$ . Solid lines represent real frequency and dashed lines represent the damping.

The next constraint on the choice of $|k|$ values is that the modes need to be periodic across the box in the $x$ and $y$ directions. With this in mind, we choose the eight numbers for $k_{x}\equiv k_{\Vert }$ and $k_{y}\equiv k_{\bot }$ shown in table 1. The table shows values in units of $kd_{i}$ as well as number of wavelengths in the box in parentheses. The complete eigenvalues for these modes are given in tables 2–5.

Table 1. The parallel ( $k_{x}$ ) and perpendicular ( $k_{y}$ ) wavenumbers in units of $kd_{i}$ as well as the number of wavelengths in the box (in parentheses). Corresponding $|k|$ and ${\it\theta}_{kB}$ (the angle between $\boldsymbol{k}$ and $\boldsymbol{B}_{0}$ ) are also shown.

Table 2. Frequencies and damping rates for the suggested KAW modes.

Table 3. Electric field vectors for the KAW modes. Here, and in the following tables, subscript $z$ denotes the direction parallel to the mean field $\boldsymbol{B}_{0}$ , the $(y,z)$ plane contains the wavevector $\boldsymbol{k}$ and the $x$ direction is perpendicular to $y,z$ . Subscript $k$ denotes the direction parallel to the $\boldsymbol{k}$ vector.

Table 4. Magnetic field for the selected KAW wavevectors.

Table 5. KAW velocity vectors for ions as well as electrons. In each group, the first line is for ions and the second line is for electrons.

Ideally each kinetic model would, as a first step, set up a single KAW and compare its properties with the predictions of linear theory. However, this should be treated as a basic test for the computer codes, as the chances of gaining new physical insights will be low with this test. A non-trivial set-up is required to compare different computer models. Hence, for the comparison of different computer models, the initial condition will be a spectrum of KAWs with the above-mentioned $|k|$ and ${\it\theta}_{kB}$ created using the prescription in § 3 of Gary & Nishimura (2004). The initial field fluctuations are written as

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}\boldsymbol{B}(x,t=0)=\mathop{\sum }_{{\it\alpha}=x,y,z}\hat{\boldsymbol{e}}_{{\it\alpha}}({\it\delta}B_{{\it\alpha}})_{0}\sin (k_{x}x+k_{y}y+{\it\phi}_{B{\it\alpha}}), & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}\boldsymbol{E}(x,t=0)=\mathop{\sum }_{{\it\alpha}=x,y,z}\hat{\boldsymbol{e}}_{{\it\alpha}}({\it\delta}E_{{\it\alpha}})_{0}\sin (k_{x}x+k_{y}y+{\it\phi}_{E{\it\alpha}}), & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}\boldsymbol{v}_{j}(x,t=0)=\mathop{\sum }_{{\it\alpha}=x,y,z}\hat{\boldsymbol{e}}_{{\it\alpha}}({\it\delta}v_{j{\it\alpha}})_{0}\sin (k_{x}x+k_{y}y+{\it\phi}_{vj{\it\alpha}}), & \displaystyle\end{eqnarray}$$

where the mode amplitudes are provided by the linear Vlasov code (Gary 1986). The phases are chosen randomly. We expect the kinetic codes to quickly adjust the phases to the real phases of the fluctuations. The total amplitude of the fluctuations will be such that $|{\it\delta}B|^{2}/B_{0}^{2}\sim 0.1$ , with the fluctuation amplitudes chosen to have the same value for every mode. PIC simulations of initially flat narrowband spectra of whistler fluctuations (Chang et al. 2011, 2015; Gary et al. 2012) show that the fluctuations evolve to broadband turbulent-like spectra, which decrease with increasing $k_{\bot }$ . These results suggest that the forward cascade determines the late-time properties of the spectra, and that those spectral properties are relatively independent of the initial distributions of fluctuation amplitudes. The wavevectors output from the linear Vlasov code (Gary 1986) for the given parameters are listed in tables 1–5.

4. Common diagnostics

To facilitate a quantitative comparison between different simulation models, a common set of diagnostics will be performed on the simulations. In § 4.1 we list the suggested diagnostics for the simulations described in § 3.1, and in § 4.2 we list the suggested diagnostics for the simulations described in § 3.2. Although the diagnostics are listed separately for the two problems, we suggest that these be performed on all the simulations if possible. Some problem-independent diagnostics are as follows:

1. (a) We will plot the change in the thermal energy of protons as a fraction of initial free energy, defined as $E_{tot}-E_{B_{0}}$ , where $E_{tot}=E_{B}+E_{i}+E_{e}+E_{E}$ is the total energy, which includes contributions from electromagnetic fluctuations as well as plasma and $E_{B_{0}}$ is the energy in the mean magnetic field, available at the beginning of the simulation. Where possible, anisotropy as defined by $T_{\bot }/T_{\Vert }$ , with $\bot$ and $\Vert$ defined with respect to the mean field, will also be plotted. As an example, figure 2 shows the change in thermal energy as well as the temperature anisotropy for three hybrid simulations of an Orszag–Tang vortex (OTV) with different equation of state for the electrons (for details see Parashar, Vasquez & Markovskii 2014). This way of plotting makes the results independent of the units used in the simulation code and hence makes a comparison of proton heating in different codes simple.

2. (b) Power spectra for magnetic field should also be compared from all the models. The spectra will be decomposed into $E_{B}(k_{\bot })$ and $E_{B}(k_{\Vert })$ and the slopes of these spectra will be computed.

We now describe the measures for intermittency and wave physics.

Figure 2. Change in thermal energy as fraction of initially available free energy as well as the anisotropy of proton heating for three hybrid simulations of OTV with: cold isothermal (CI), warm isothermal (WI) and warm adiabatic (WA) electrons. The system is turbulent after $t\sim 70{\it\omega}_{ci}^{_{1}}$ . From Parashar et al. (2014). Copyright American Institute of Physics.

4.1. Intermittent structures

There are multiple different measures used to quantify intermittency of turbulence, e.g. kurtosis of derivatives, scale-dependent kurtosis, filtered kurtosis, probability distribution functions (PDFs) of increments (e.g. Greco et al. 2008; Parashar et al. 2011), local intermittency measure (LIM) (Farge et al. 1990; Farge 1992; Bruno et al. 1999), phase coherence index (e.g. Hada et al. 2003; Koga et al. 2007) and partial variance of increments (PVI) (e.g. Greco et al. 2008), and structure functions (e.g. Bruno & Carbone 2013). In particular we will use scale-dependent kurtosis, PDFs of increments and structure functions to quantify intermittency in the simulations described in § 3.1. We now describe these in more detail.

4.1.1. PDFs of increments

The PDFs of turbulent quantities are typically Gaussian (e.g. Frisch 1996) but the PDFs of increments of a turbulent quantity are not Gaussian. By taking the increments, say for the magnetic field, $\boldsymbol{B}(s+{\it\delta}s)-\boldsymbol{B}(s)$ , the gradients (and hence the intermittent structures) are highlighted in the quantity of interest. The strength of the gradients highlighted depends on the lag ${\it\delta}s$ . Hence the increment for a smaller lag ${\it\delta}s$ represents steeper gradients and hence most intermittent structures. When ${\it\delta}s$ becomes comparable to the correlation length of the system, the PDFs typically revert back to approximate Gaussianity.

It has been shown that the non-Gaussian tails on the PDFs of increments correspond to the number of intermittent structures (e.g. Greco et al. 2008, 2009; Salem et al. 2009; Wan et al. 2010). We show an example of this from an MHD simulation in figure 3. Hence by comparing PDFs of a fixed increment ( ${\it\delta}s\sim 1d_{i}$ ), from multiple models, the number of intermittent structures resolved in each model can be compared.

Figure 3. (a,b) The out-of-plane current from a $2048^{2}$ MHD simulation, showing the PDFs for this current. The PDFs clearly show super-Gaussian tails. (c) The currents corresponding to values consistent with the super-Gaussian tails in panel (b). It is clear that the strongest currents are associated with the super-Gaussian tails.

4.1.2. Scale-dependent kurtosis

The $n$ th-order moment of the PDFs of a variable ${\it\phi}$ , also called the $n$ th-order structure function, is defined as

(4.1) $$\begin{eqnarray}S_{{\it\phi}}^{n}(\ell )=\langle |{\it\phi}(x+\ell )-{\it\phi}(x)|^{n}\rangle =\langle {\it\delta}{\it\phi}_{\ell }^{n}\rangle .\end{eqnarray}$$

The kurtosis, the fourth-order moment or fourth structure function, measures the Flatness of the PDFs of a given function ${\it\phi}$ . It is defined as

(4.2) $$\begin{eqnarray}{\it\chi}_{{\it\phi}}=\langle {\it\delta}{\it\phi}_{\ell }^{4}\rangle /\langle {\it\delta}{\it\phi}_{\ell }^{2}\rangle ^{2},\end{eqnarray}$$

where ${\it\phi}$ is the quantity of interest. The lag $\ell$ is properly interpreted as a spatial lag, which should be normalized either to the correlation scale (as in Greco et al. 2009) or to the ion inertial scale (as in Wu et al. 2013). Such normalizations facilitate the comparison across codes, with different system sizes, and with spacecraft data in the super-Alfvénic solar wind. With single-spacecraft data, the lag would be defined in terms of a time lag, multiplied by a solar wind speed. In 2.5D or 3D simulations, one could use all simulation data with separations in one Cartesian direction corresponding to  $\ell$ .

Departures from Gaussianity are measured by ${\it\chi}$ . For a Gaussian function, ${\it\chi}$ is equal to 3, and it is generally greater than 3 for increments of turbulent quantities. Greater non-Gaussianity is expected with smaller $\ell$ . A higher value of kurtosis reveals the presence of larger concentrations of coherent structures. Measuring ${\it\chi}_{{\it\phi}}$ as a function of $\ell$ for different variables is a straightforward diagnostic to implement on the simulation data. Since ${\it\chi}_{{\it\phi}}$ is a dimensionless ratio, its value is independent of the units and/or normalizations used in the simulation codes.

4.1.3. Structure functions and their scaling exponents

Within the inertial range, the structure functions defined in (4.1) scale as power laws,

(4.3) $$\begin{eqnarray}S_{{\it\phi}}^{n}(\ell )=\langle {\it\delta}{\it\phi}_{\ell }^{n}\rangle \propto \ell ^{{\it\zeta}_{n}}.\end{eqnarray}$$

The scaling exponents ${\it\zeta}_{n}$ are measures of intermittency in hydrodynamic and MHD turbulence. Estimating the kurtosis (defined above) requires calculating structure functions up to the fourth order. When it is numerically feasible, it would be worthwhile to go up to sixth order and measure the first six scaling exponents for $v$ , $B$ and density over the scale range $2d_{i}<\ell <20d_{i}$ in order to provide a more complete description of the intermittent fluctuations. Measuring these scaling exponents will also make it possible to tie into a larger literature on intermittency (e.g. Bruno & Carbone 2013; Salem et al. 2009; Chandran, Schekochihin & Mallet 2015). Structure functions are also widely used in the analysis of solar wind fluctuations from spacecraft data, for both inertial range (e.g. Salem et al. 2009) and dissipation range fluctuations (e.g. Osman et al. 2014).

While higher-order statistics are needed for a full description and comparison of turbulence properties, one must also realize that there are inherent difficulties. The main problem is perhaps that, for any model (or even experimental data), higher-order moments become more difficult to estimate (e.g. De Wit 2004); this problem is even worse for odd moments that involve substantial cancellations (e.g. Podesta et al. 2009). In addition, there will be difficulties when comparing across different codes, and examining these will be a greater focus of this challenge. It is already known that differences are encountered across code types (e.g. Haugbølle et al. 2013; Henri et al. 2013).

4.1.4. Choice of second- and higher-order statistics

In view of the above difficulties, meaningful comparisons across code types should begin with examination of the first few moments of the PDFs. Second-order moments (spectra, correlation functions) are already a point of emphasis in §§ 4.2.1 and 4. Third-order moments are non-positive-definite and therefore very difficult (e.g. Podesta et al. 2009). Therefore the higher-order moment of greatest importance in the comparisons will be initially the fourth-order one, i.e. the scale-dependent kurtosis described above. Until such comparisons are established, it makes little sense to focus on multifractal scalings, which can be addressed once the fourth-order comparisons are understood.

4.2. Wave physics

Here we list suggested diagnostics for wave simulations described in § 3.2.

4.2.1. Spectral anisotropy

Turbulent plasmas with mean magnetic field typically show a distribution of energy in high $k_{\bot }$ modes (e.g. Higdon 1984; Goldreich & Sridhar 1995; Montgomery & Matthaeus 1995; Oughton, Dmitruk & Matthaeus 2006), etc. This property has also been observed to hold at kinetic scales (e.g. Chen et al. 2010; Chang et al. 2013). It will be instructive to compare the spectral anisotropy in the $k_{\bot }{-}k_{\Vert }$ plane produced by different models.

4.2.2. Dispersion analysis

An obvious test for the presence of waves is to look for the appropriate dispersion in the energy spectrum. A $k{-}{\it\omega}$ dispersion analysis of the simulation data will show excess of undamped energy along the dispersion curves of the normal modes of the system. As an example, figure 4 shows the spectrum of the magnetic field as a function of $k_{\Vert },{\it\omega}$ and $k_{\bot },{\it\omega}$ from a 2.5D hybrid simulation with ${\it\beta}=0.04$ , large system size, mean field in the plane of simulation and driven at $|k|=2,3$ (from Parashar 2011). The two-fluid dispersion curves for parallel and perpendicular propagation have been over-plotted. Outside the driving wavenumbers, there is clearly enhanced energy along the dispersion curves, indicating the presence of waves.

Figure 4. The $k,{\it\omega}$ spectra of magnetic energy from a 2.5D hybrid PIC simulation, ${\it\beta}=0.04$ driven at large scales of the system. The normal modes of the system are excited and we see enhanced energy along the dispersion curves of Alfvén, slow and magnetosonic modes (from Parashar 2011).

4.2.3. Damping rate

One of the most important features to check is if the damping rate of the waves is appropriately captured by the simulation codes. A comparison of numerically calculated damping rates to the average damping rate calculated from linear Vlasov theory is an important consistency check. It should be noted that, given the noisy nature of PIC algorithms, the damping rate comparisons are not expected to be extremely accurate. We will follow the procedure of Chang et al. (2014), who showed that the integrated damping rate of whistler turbulence is consistent with the predictions of linear kinetic dispersion theory of whistler turbulence.

4.2.4. Compressibility

Several distinct compressibilities can be defined for a plasma (e.g. Gary & Smith 2009). The plasma compressibility of the $j$ th plasma species is defined as

(4.4) $$\begin{eqnarray}C_{j}(\boldsymbol{k})\equiv \frac{|{\it\delta}n_{j}(\boldsymbol{k})|^{2}}{n_{0}^{2}}\frac{B_{0}^{2}}{|{\it\delta}\boldsymbol{B}(\boldsymbol{k})|^{2}}.\end{eqnarray}$$

For a quasi-neutral plasma consisting of only electrons and protons, $C_{e}\sim C_{p}$ , hence it is sufficient to study the compressibility of only a single species. Similarly, the compressibility of the magnetic field can also be computed:

(4.5) $$\begin{eqnarray}C_{\Vert }(\boldsymbol{k})\equiv \frac{|{\it\delta}B_{\Vert }(\boldsymbol{k})|^{2}}{|{\it\delta}\boldsymbol{B}(\boldsymbol{k})|^{2}}.\end{eqnarray}$$

The compressibility of different modes (Alfvén–kinetic Alfvén and magnetosonic–whistler) can be calculated from linear dispersion theory and compared with simulations (e.g. Saito et al. 2008) and observations (e.g. Gary & Smith 2009; Salem et al. 2012). The differences between these quantities are useful not only in identifying the fluctuating constituents of observed short-wavelength turbulence (e.g. Salem et al. 2012) but also as a means of testing a code’s capability to accurately capture the wave physics of a turbulent simulation.

4.2.5. The electric to magnetic field ratio

The electric to magnetic field ratio, $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ (Salem et al. 2012), provides a complementary test for wave-mode identification purposes. The $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ ratio is constant in the inertial range as the electric and magnetic field fluctuations are well correlated (Bale et al. 2005). Indeed, the electric field $\boldsymbol{E}$ in this range is essentially equivalent to $\boldsymbol{v}\times \boldsymbol{B}$ . At scales smaller than the ion gyro-scale, $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ increases as a function of $k$ owing to the dispersive nature of the modes. In linear theory, $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ depends on the propagation angle of the mode, in both the inertial and dissipation range.

As for the compressibility, a direct comparison of the $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ expected from linear dispersion theory with the $|{\it\delta}\boldsymbol{E}|/|{\it\delta}\boldsymbol{B}|$ captured by the simulations will be a useful tool to understand the nature of the fluctuations.

4.3. Comparison with observations

At this initial stage, we cannot expect to see close similarities between the simulations and observations. However, it will still be instructive to make direct comparisons. Many of the diagnostics described above cannot be performed on spacecraft data because spacecraft make 1D cuts through the solar wind. Thus most spacecraft observations are unable to provide $k,{\it\omega}$ diagrams or assumption-free multi-dimensional spectra. Instead we suggest that the modellers take artificial 1D cuts across simulations and use these to compare with the results of solar wind observations. These ‘artificial observations’ can also be provided to observers for more direct comparisons with solar wind observations.

Spacecraft are also unable to observe all of the properties of the plasma that can be calculated from a simulation, for example magnetic and electric fields are not typically measured at scales equivalent to $|k|\sim 100d_{i}^{-1}$ (although some spacecraft are able to do so). We therefore recommend a downsampling or reduction in the parameters that are directly compared to observations, for example, current is not usually calculated from solar wind data since 3D magnetic field gradients are impossible to calculate from single spacecraft, instead magnetic and electric fields and particle velocity, density and temperature are the best data products to use.

The derived parameters that are most suitable for direct comparison between simulations and observations are magnetic field power spectra, compressibilities and magnetic helicity. Structure functions, and therefore PVI, kurtosis, the scaling of higher-order structure functions and the PDF of increments, can all be calculated from both types of data. Correlations between different quantities such as reduced energy and cross-helicity calculated from artificial observations can be compared to real observations also.

This will help set up the communication between simulation modellers and the observers for future ‘critical simulations’.

5. Envisioned roadmap

The simulations outlined in this paper will be only the first step in a multi-step process that will lead to a better understanding of proton scale dissipative processes. Although the exact path taken by this community effort will be decided based on the outcome of this step, we here outline a roadmap for the challenge.

5.1. Code comparison

Comparing the results, specifically the statistical properties, from multiple different models in a turbulent setting will enable us to look critically at the fundamental assumptions, strengths as well as limitations, of multiple simulation models. A few sketchy comparisons exist (e.g. Parashar 2011) but a thorough comparison is required.

5.2. Designing critical 3D simulations

Through a comparison of the 2.5D simulations and a few preliminary 3D simulations, a partial outcome of this effort will be important insights about the relative importance of 3D couplings that are missing in the 2.5D picture. Based on these insights, appropriate 3D simulations can be designed.

For the next step of designing the simulations, significant input from observers will be required. Highlighting a few critical solar wind intervals, parametrizing them and listing important unanswered questions will be required before the final, critical simulations can be designed and performed.

5.3. Diagnostics for reconnection 3D

An important unsolved problem is the role of reconnection in solar wind turbulence. We anticipate that future stages in this challenge will address this question, but a significant parallel effort will be required to design appropriate analysis tools applicable to 3D simulations. For example, identifying a reconnection site is relatively easy and very well documented (e.g. Servidio et al. 2009, 2010) but the nature of reconnection sites changes significantly in 3D (e.g. Daughton et al. 2011). Hence, novel diagnostic techniques will need to be developed.

5.4. Important physics questions

Once the critical simulations have been designed and opened to the community for analysis, effort will be devoted towards understanding the nature of kinetic processes in such systems. Here we provide a partial list of questions that the community could focus on.

1. (a) What fraction of the total dissipation power in the solar wind comes from low-frequency turbulence, and what fraction comes from high-frequency waves through, for example, cyclotron heating?

   An answer to this question might require studies beyond the scope of this multi-year project, as the nature of fluctuations at kinetic scales will depend on many large-scale parameters and the nature of the cascade under given conditions.

2. (b) What kinetic mechanisms damp low-frequency turbulence at small scales?

   This is the main emphasis of our present effort. Carefully designed large and fully 3D simulations are our best bet to address this question.

3. (c) How is the turbulent heating power divided between protons and electrons, and between parallel and perpendicular heating?

4. (d) How do proton/electron/alpha heating rates in numerical simulations compare to rates inferred on the basis of linear wave damping or nonlinear mechanisms such as reconnection or stochastic heating?

5. (e) What are the quantitative measures of intermittency of fluctuations and dissipation in the solar wind and in different types of numerical simulations?

6. (f) Can observations rule some mechanisms out? This is where the proposed comparisons with observations become critical. Cuts from the critical 3D simulations will be taken and provided to observers as ‘artificial spacecraft data’. Analysis performed on these artificial data sets by multiple different groups not only will help benchmark observational techniques, but also will help benchmark physics contained in simulations with reality. Hence, it is one of the critical parts of longer-term goals of the challenge.

6. Conclusion

A quantitative improvement in our understanding of the dissipative processes active in collisionless plasmas such as the solar wind will benefit from a better understanding of the strengths and limitations of different types of numerical simulations. To aid in the development of this understanding, we have proposed two 2.5D ‘challenge problems’ that can be studied using different simulation techniques. The first problem is to study the evolution of KH instability in a plane that is nearly perpendicular to the mean magnetic field. The second problem is to study kinetic Alfvén waves in the plane containing $\boldsymbol{B}_{0}$ . We also carefully described the diagnostics that should be performed on all simulations. To investigate the generation of, and dissipation at, intermittent structures, we require the measurement of the PDF of increments, kurtosis and PVI of the simulations over a range of scales. For the wave dissipation task we have determined that dispersion relations can be tested in simulations, but when comparing to observations other techniques must be used, and we encourage those analyzing simulations to make mock spacecraft trajectories and also measure the properties of the particle distributions, the magnetic fluctuation spectra, helicity and compressibility.

Initially, all of these simulations will be performed in 2.5D owing to a desire to have large Reynolds numbers but maintain high accuracy with the various methods employed by the modelling community. Once these initial studies are complete, we aim to encourage well-targeted, large, 3D simulations to access the truly 3D nature of turbulence and dissipation. These simulations will ideally be the largest possible simulations with initial conditions designed to mimic real solar wind conditions and based upon the results found from studies designed as described above. The data obtained from these critical simulations will be made open to the public in order for the wider research community to make quantitative comparisons between different dissipative processes.

The ‘turbulent dissipation challenge’ aims to help coordinate community efforts, so that we can collectively achieve greater progress on these important and difficult questions. We hope that the plasma modelling and solar wind observation research communities will take up this challenge and use the conditions and measurement techniques described here as a place to start a coherent and well-benchmarked campaign of studies to make considerable advances in our understanding of the dissipation of turbulence and collisionless plasma heating.