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 Reference Kolmogorov1941). The solar wind magnetic energy spectrum obeys the Kolmogorov power law in the ‘inertial range’ (e.g. Coleman Reference Coleman1968; Matthaeus & Goldstein Reference Matthaeus and Goldstein1982; Goldstein, Roberts & Matthaeus Reference Goldstein, Roberts and Matthaeus1995), 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. Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009; Alexandrova et al. Reference Alexandrova, Lacombe, Mangeney, Grappin and Maksimovic2012). 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 Reference Ghosh and Goldstein1997; Stawicki, Gary & Li Reference Stawicki, Gary and Li2001; Boldyrev & Perez Reference Boldyrev and Perez2012), 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 Reference Marsch2006) 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 Reference Hollweg and Isenberg2002, and references therein), Landau damping (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009, and references therein), particle energization at current sheets and reconnection sites (e.g. Dmitruk, Matthaeus & Seenu Reference Dmitruk, Matthaeus and Seenu2004; Sundkvist et al. Reference Sundkvist, Retinò, Vaivads and Bale2007; Parashar et al. Reference Parashar, Shay, Cassak and Matthaeus2009, Reference Parashar, Servidio, Shay, Breech and Matthaeus2011; Osman et al. Reference Osman, Matthaeus, Greco and Servidio2011; Wan et al. Reference Wan, Matthaeus, Karimabadi, Roytershteyn, Shay, Wu, Daughton, Loring and Chapman2012; Karimabadi et al. Reference Karimabadi, Roytershteyn, Wan, Matthaeus, Daughton, Wu, Shay, Loring, Borovsky, Leonardis, Chapman and Nakamura2013) and stochastic heating (e.g. McChesney, Stern & Bellan Reference McChesney, Stern and Bellan1987; Chaston et al. Reference Chaston, Bonnell, Carlson, McFadden, Ergun, Strangeway and Lund2004; Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010; Bourouaine & Chandran Reference Bourouaine and Chandran2013; Xia et al. Reference Xia, Perez, Chandran and Quataert2013). 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. Reference TenBarge, Howes and Dorland2013) and stochastic heating (Chandran et al. Reference Chandran, Li, Rogers, Quataert and Germaschewski2010) 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. Reference Schekochihin, Cowley, Rincon and Rosin2010), accretion disks (e.g. Kunz, Schekochihin & Stone Reference Kunz, Schekochihin and Stone2014), astrophysical jets (e.g. Romanova & Lovelace Reference Romanova and Lovelace1992), the interstellar medium, planetary magnetospheres, geologically active moons within planetary magnetospheres or magnetic moons within planetary magnetospheres (e.g. Saur et al. Reference Saur, Politano, Pouquet and Matthaeus2002; Mousis & Gautier Reference Mousis and Gautier2004; Saur Reference Saur2004; Kivelson Reference Kivelson2006; von Papen, Saur & Alexandrova Reference von Papen, Saur and Alexandrova2014).

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 Reference Kingsep, Chukbar and Ian’kov1987; Gordeev, Kingsep & Rudakov Reference Gordeev, Kingsep and Rudakov1994; Cho & Lazarian Reference Cho and Lazarian2004) 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. Reference Parashar, Shay, Cassak and Matthaeus2009) treat protons as particles and electrons as a neutralizing fluid (usually massless and isothermal). Gyro-kinetics (e.g. Brizard & Hahm Reference Brizard and Hahm2007) 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 Reference Gurnett and Bhattacharjee2005), Hall MHD (e.g. Ghosh & Goldstein Reference Ghosh and Goldstein1997), Hall-FLR MHD (e.g. Mjølhus Reference Mjølhus2009) with finite Larmor radius (FLR) corrections, reduced MHD (RMHD) (e.g. Zank & Matthaeus Reference Zank and Matthaeus1992), electron reduced MHD (ERMHD) (e.g. Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), Landau fluid (e.g. Goswami, Passot & Sulem Reference Goswami, Passot and Sulem2005; Hunana et al. Reference Hunana, Goldstein, Passot, Sulem, Laveder and Zank2013), hybrid Eulerian Vlasov (e.g. Valentini, Califano & Veltri Reference Valentini, Califano and Veltri2010) and full PIC (e.g. Karimabadi et al. Reference Karimabadi, Roytershteyn, Wan, Matthaeus, Daughton, Wu, Shay, Loring, Borovsky, Leonardis, Chapman and Nakamura2013).

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 Reference Harten and Clark1995), CLUSTER (Escoubet, Fehringer & Goldstein Reference Escoubet, Fehringer and Goldstein2001), ACE (Stone et al. Reference Stone, Frandsen, Mewaldt, Christian, Margolies, Ormes and Snow1998), Helios (Rosenbauer et al. Reference Rosenbauer, Miggenrieder, Montgomery and Schwenn1976), Ulysses (Wenzel et al. Reference Wenzel, Marsden, Page and Smith1992) and Voyager (Behannon et al. Reference Behannon, Acuna, Burlaga, Lepping, Ness and Neubauer1977), 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. Reference Henri, Cerri, Califano, Pegoraro, Rossi, Faganello, Šebek, Trávníček, Hellinger, Frederiksen, Nordlund, Markidis, Keppens and Lapenta2013; Keppens et al. Reference Keppens, Porth, Galsgaard, Frederiksen, Restante, Lapenta and Parnell2013), 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. Reference Greco, Chuychai, Matthaeus, Servidio and Dmitruk2008; Wicks et al. Reference Wicks, Horbury, Chen and Schekochihin2010, Reference Wicks, Forman, Horbury and Oughton2012; Wu et al. Reference Wu, Perri, Osman, Wan, Matthaeus, Shay, Goldstein, Karimabadi and Chapman2013). Hence, we choose to work with systems with a fairly large Reynolds numberFootnote ² 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.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. Reference Svidzinski, Li, Rose, Albright and Bowers2009) and whistler waves (e.g. Saito et al. Reference Saito, Gary, Li and Narita2008; Gary, Chang & Wang Reference Gary, Chang and Wang2012; Chang, Gary & Wang Reference Chang, Gary and Wang2013, Reference Chang, Gary and Wang2014). 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. Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsuno2008; Sahraoui et al. Reference Sahraoui, Goldstein, Robert and Khotyaintsev2009; Shaikh & Zank Reference Shaikh and Zank2009; Podesta, Borovsky & Gary Reference Podesta, Borovsky and Gary2010; Salem et al. Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012; Chen et al. Reference Chen, Boldyrev, Xia and Perez2013).

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 Reference Gary, Saito and Li2008; Saito et al. Reference Saito, Gary, Li and Narita2008; Saito, Gary & Narita Reference Saito, Gary and Narita2010; Saito & Gary Reference Saito and Gary2012), fully 3D whistler turbulence (e.g. Gary et al. Reference Gary, Chang and Wang2012; Chang, Gary & Wang Reference Chang, Gary and Wang2011, Reference Chang, Gary and Wang2013, Reference Chang, Gary and Wang2014, Reference Chang, Gary and Wang2015) and 2.5D magnetosonic turbulence (e.g. Svidzinski et al. Reference Svidzinski, Li, Rose, Albright and Bowers2009). 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 Reference Markovskii and Vasquez2010; Markovskii, Vasquez & Chandran Reference Markovskii, Vasquez and Chandran2010) and Alfvénic and kinetic Alfvénic turbulence in full 3D (e.g. Vasquez, Markovskii & Chandran Reference Vasquez, Markovskii and Chandran2014) geometries. Tronko et al. (Reference Tronko, Nazarenko and Galtier2013) and Howes (Reference Howes2015) 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. Reference Verscharen, Marsch, Motschmann and Müller2012; Comişel et al. Reference Comişel, Verscharen, Narita and Motschmann2013; Comişel, Narita & Motschmann Reference Comişel, Narita and Motschmann2014) 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 Reference Gary1986, Reference Gary2005), 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 (Reference Gary and Nishimura2004). 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 Reference Gary1986). 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. Reference Chang, Gary and Wang2011, Reference Chang, Gary and Wang2015; Gary et al. Reference Gary, Chang and Wang2012) 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 Reference Gary1986) 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 Reference Parashar, Vasquez and Markovskii2014). 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. (Reference Parashar, Vasquez and Markovskii2014). 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. Reference Greco, Chuychai, Matthaeus, Servidio and Dmitruk2008; Parashar et al. Reference Parashar, Servidio, Shay, Breech and Matthaeus2011), local intermittency measure (LIM) (Farge et al. Reference Farge, Guezennec, Ho, Meneveau and Spinks1990; Farge Reference Farge1992; Bruno et al. Reference Bruno, Bavassano, Bianchini, Pietropaolo, Villante, Carbone, Veltri and Wilson1999), phase coherence index (e.g. Hada et al. Reference Hada, Koga and Yamamoto2003; Koga et al. Reference Koga, Chian, Miranda and Rempel2007) and partial variance of increments (PVI) (e.g. Greco et al. Reference Greco, Chuychai, Matthaeus, Servidio and Dmitruk2008), and structure functions (e.g. Bruno & Carbone Reference Bruno and Carbone2013). 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 Reference Frisch1996) 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. Reference Greco, Chuychai, Matthaeus, Servidio and Dmitruk2008, Reference Greco, Matthaeus, Servidio, Chuychai and Dmitruk2009; Salem et al. Reference Salem, Mangeney, Bale and Veltri2009; Wan et al. Reference Wan, Oughton, Servidio and Matthaeus2010). 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. Reference Greco, Matthaeus, Servidio, Chuychai and Dmitruk2009) or to the ion inertial scale (as in Wu et al. Reference Wu, Perri, Osman, Wan, Matthaeus, Shay, Goldstein, Karimabadi and Chapman2013). 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 Reference Bruno and Carbone2013; Salem et al. Reference Salem, Mangeney, Bale and Veltri2009; Chandran, Schekochihin & Mallet Reference Chandran, Schekochihin and Mallet2015). 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. Reference Salem, Mangeney, Bale and Veltri2009) and dissipation range fluctuations (e.g. Osman et al. Reference Osman, Kiyani, Chapman and Hnat2014).

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 Reference De Wit2004); this problem is even worse for odd moments that involve substantial cancellations (e.g. Podesta et al. Reference Podesta, Forman, Smith, Elton, Malécot and Gagne2009). 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. Reference Haugbølle, Frederiksen and Nordlund2013; Henri et al. Reference Henri, Cerri, Califano, Pegoraro, Rossi, Faganello, Šebek, Trávníček, Hellinger, Frederiksen, Nordlund, Markidis, Keppens and Lapenta2013).

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. Reference Podesta, Forman, Smith, Elton, Malécot and Gagne2009). 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 Reference Higdon1984; Goldreich & Sridhar Reference Goldreich and Sridhar1995; Montgomery & Matthaeus Reference Montgomery and Matthaeus1995; Oughton, Dmitruk & Matthaeus Reference Oughton, Dmitruk and Matthaeus2006), etc. This property has also been observed to hold at kinetic scales (e.g. Chen et al. Reference Chen, Horbury, Schekochihin, Wicks, Alexandrova and Mitchell2010; Chang et al. Reference Chang, Gary and Wang2013). 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 Reference Parashar2011). 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 Reference Parashar2011).

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. (Reference Chang, Gary and Wang2014), 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 Reference Gary and Smith2009). 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. Reference Saito, Gary, Li and Narita2008) and observations (e.g. Gary & Smith Reference Gary and Smith2009; Salem et al. Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012). The differences between these quantities are useful not only in identifying the fluctuating constituents of observed short-wavelength turbulence (e.g. Salem et al. Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012) 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. Reference Salem, Howes, Sundkvist, Bale, Chaston, Chen and Mozer2012), 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. Reference Bale, Kellogg, Mozer, Horbury and Reme2005). 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.

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.