Kinetic Plasma Turbulence: Recent Insights and Open Questions From 3D3V Simulations

Cerri, Silvio Sergio; Grošelj, Daniel; Franci, Luca

doi:10.3389/fspas.2019.00064

PERSPECTIVE article

Front. Astron. Space Sci., 11 October 2019

Sec. Space Physics

Volume 6 - 2019 | https://doi.org/10.3389/fspas.2019.00064

This article is part of the Research Topic Improving the Understanding of Kinetic Processes in Solar Wind and Magnetosphere: From CLUSTER to MMS View all 10 articles

Kinetic Plasma Turbulence: Recent Insights and Open Questions From 3D3V Simulations

$\nSilvio Sergio Cerri$ Silvio Sergio Cerri¹^*

Daniel Grošelj²

Luca Franci³

¹Department of Astrophysical Sciences, Princeton University, Princeton, NJ, United States
²Max-Planck-Institut für Plasmaphysik, Garching, Germany
³School of Physics and Astronomy, Queen Mary University of London, London, United Kingdom

Turbulence and kinetic processes in magnetized space plasmas have been extensively investigated over the past decades via in-situ spacecraft measurements, theoretical models and numerical simulations. In particular, multi-point high-resolution measurements from the Cluster and MMS space missions brought to light an entire new world of processes, taking place at the plasma kinetic scales, and exposed new challenges for their theoretical interpretation. A long-lasting debate concerns the nature of ion and electron scale fluctuations in solar-wind turbulence and their dissipation via collisionless plasma mechanisms. Alongside observations, numerical simulations have always played a central role in providing a test ground for existing theories and models. In this Perspective, we discuss the advances achieved with our 3D3V (reduced and fully) kinetic simulations, as well as the main questions left open (or raised) by these studies. To this end, we combine data from our recent kinetic simulations of both freely decaying and continuously driven fluctuations to assess the similarities and/or differences in the properties of plasma turbulence in the sub-ion range. Finally, we discuss possible future directions in the field and highlight the need to combine different types of numerical and observational approaches to improve the understanding of turbulent space plasmas.

1. Introduction

With the establishment of satellite space missions, the near-Earth environment and the solar wind have provided unique opportunities to explore the physics of weakly collisional, magnetized plasmas (e.g., Bruno and Carbone, 2013; Chen, 2016; Verscharen et al., 2019). In particular, increasingly accurate in-situ measurements of plasma fluctuations and particle distribution functions from Cluster and MMS have uncovered an entire new world of kinetic processes occurring in plasma turbulence (e.g., Alexandrova et al., 2009, 2012, 2013; Sahraoui et al., 2009, 2010; Chen et al., 2010, 2019; Greco et al., 2016; Narita et al., 2016; Chasapis et al., 2017; Chen and Boldyrev, 2017; Huang et al., 2017; Roberts et al., 2017; Servidio et al., 2017). These observations highlight a change in the turbulent cascade at plasma microscales, challenging the community for a consistent theory of kinetic-range turbulence. In fact, several collisionless plasma processes may be simultaneously at play and compete with each other in determining the nature of ion-scale and electron-scale fluctuations (e.g., Stawicki et al., 2001; Galtier and Bhattacharjee, 2003; Howes et al., 2008a; Gary and Smith, 2009; Schekochihin et al., 2009; Boldyrev and Perez, 2012; He et al., 2012; Podesta, 2012; Boldyrev et al., 2013; Matthaeus et al., 2014; Passot and Sulem, 2015, 2019; Kunz et al., 2018; Passot et al., 2018) and, consequently, how free energy cascades in phase space (e.g., Schekochihin et al., 2008; Servidio et al., 2017; Adkins and Schekochihin, 2018; Cerri et al., 2018; Eyink, 2018; Pezzi et al., 2018; Kawazura et al., 2019). Many observations at ion and sub-ion scales, specifically, suggest that turbulent fluctuations exhibit properties mainly typical of kinetic Alfvén waves (KAWs) (Leamon et al., 1998; Sahraoui et al., 2009; Podesta and TenBarge, 2012; Salem et al., 2012; Chen et al., 2013; Kiyani et al., 2013; Chen, 2016; Lacombe et al., 2017). The emergence of KAW-like fluctuations in kinetic turbulence has been also supported by means of a large number of theoretical and numerical works (e.g., Hollweg, 1999; Stawicki et al., 2001; Gary and Nishimura, 2004; Howes et al., 2008a; Gary and Smith, 2009; Sahraoui et al., 2012; TenBarge et al., 2012; Vásconez et al., 2014, 2015; Franci et al., 2015; Cerri et al., 2016; Pucci et al., 2016; Zhao et al., 2016; Valentini et al., 2017; Grošelj et al., 2019). Some of these studies rely on the so-called spectral field ratios, which provide a measure of the wave-like polarization properties of the turbulent fluctuations, as compared to what linear theory predicts (see, e.g., Boldyrev et al., 2013) and section 3.

In the above context, direct numerical simulations play a key role by providing a controlled test ground for different theories, providing information not accessible to observations. Enormous efforts have been recently made to understand 3D kinetic turbulence via numerical experiments (e.g., Howes et al., 2008b; Gary et al., 2012; TenBarge and Howes, 2013; Vasquez et al., 2014; Servidio et al., 2015; Told et al., 2015; Wan et al., 2015, 2016; Bañón Navarro et al., 2016; Comişel et al., 2016; Cerri et al., 2017b, 2018; Hughes et al., 2017a,b; Kobayashi et al., 2017; Franci et al., 2018a,b; Grošelj et al., 2018, 2019; Arzamasskiy et al., 2019; Roytershteyn et al., 2019; Zhdankin et al., 2019). In this Perspective, we combine data from our recent 3D3V studies (Cerri et al., 2017b; Franci et al., 2018b; Grošelj et al., 2019) to investigate whether common turbulence features exist in all three independently performed simulations (section 2), thus indicating a certain “universality” of kinetic-scale turbulence. Moreover, we also highlight possible model-dependent differences between the 3D hybrid-kinetic and fully kinetic simulations. We mention that this approach follows the general idea of adopting different models (and/or implementations) to study turbulent heating and dissipation in collisionless plasmas that was initiated within the “Turbulent Dissipation Challenge” framework (Parashar et al., 2015). Here we extend similar comparative analysis of the spectral properties that have been previously performed in a reduced two-dimensional setup (see Cerri et al., 2017a; Franci et al., 2017; Grošelj et al., 2017) to the more realistic three-dimensional geometry (section 3), and we present a new analysis of our data based on local structure functions (section 4). Finally, we discuss possible implications for sub-ion-scale turbulence and future directions emerging from this study (section 5).

2. Data Sets

In the following, we consider three recent kinetic simulations in a six-dimensional phase space (“3D3V”) using: (i) CAMELIA, a hybrid particle-in-cell (PIC) code with massless electrons (Franci et al., 2018a), (ii) HVM, an Eulerian hybrid-Vlasov code with finite electron-inertia effects (Valentini et al., 2007), and (iii) OSIRIS, a fully kinetic PIC code (Fonseca et al., 2002, 2013). Unless otherwise specified, parallel (||) and perpendicular (⊥) directions are defined with respect to the global mean magnetic field B₀ = B₀e_z. Franci et al. (2018b) employed the CAMELIA code to investigate freely decaying, Alfvénic fluctuations in a cubic box (L_|| = L_⊥ = 128d_i with 512³ grid points and 2048 particles per cell (ppc)) for β_i = β_e = 0.5, where $β_{s} = 8 π n_{0} T_{s} / B_{0}^{2}$ is the species beta. Cerri et al. (2017b) instead adopted the HVM code to study freely decaying compressive fluctuations in an elongated box (L_|| = 2L_⊥ ≃ 63d_i with 384² × 64 grid points in real space, and 51³ points in a velocity space bounded by |v/v_{th, i}| ≤ 5) for β_i = β_e = 1 and with a reduced ion-electron mass ratio of m_i/m_e = 100 (viz. including d_e-effects in a generalized Ohm's law). Spectral filters were applied at runtime, determining a cutoff in the turbulent spectrum at k_⊥d_i > 20 and at k_zd_i > 2. Finally, Grošelj et al. (2019) use the OSIRIS code to investigate continuously driven Alfvénic fluctuations in a β_i ≈ β_e ≈ 0.5 plasma with m_i/m_e = 100. An elongated box was used (L_|| = 2.56L_⊥ ≃ 48.3d_i with 928² × 1920 grid points and 150 ppc per species). An example of $δ {\tilde{B}}_{⊥} = δ B_{⊥} / δ B_{⊥}^{(r m s)}$ in a two-dimensional cut perpendicular to B₀ is given in Figure 1A, along with a schematic representation of these simulations in the (k_⊥, k_||) plane (Figure 1B).

FIGURE 1

Figure 1. (A) $δ {\tilde{B}}_{⊥} = δ B_{⊥} / δ B_{⊥}^{(r m s)}$ in a plane perpendicular to B₀. (B) Nominal wavenumber-space representation of simulations. (C) δB_⊥ spectrum in (k_⊥, k_z) space. White dotted lines mark kd_e = 1. (D) Top panels: reduced spectra vs. k_⊥d_i. Spectra have been shifted (see text). Bottom panels: local spectral exponents. Horizontal lines denote −5/3 (dashed) and −8/3 (dash-dotted) slopes. Vertical dotted line marks k_⊥d_e = 1. (D) Same as (E), but vs. k_zd_i. Here −2 (dashed) and −7/2 (dash-dotted) slopes are marked for reference. (F) Spectral ratios vs. k_⊥, normalized to the asymptotic KAW prediction (dashed horizontal lines; see text for details).

In the following, the analysis of freely decaying simulations (viz., CAMELIA and HVM) is performed at the peak of the turbulent activity (cf., e.g., Servidio et al., 2015), while for the continuously driven OSIRIS run we consider the turbulence at the end of the simulation when the kinetic range spectra appear converged. Following Franci et al. (2018b) and Grošelj et al. (2019), PIC data have been filtered before performing the analysis to remove spectral regions dominated by particle noise. The OSIRIS data have been filtered for k_⊥d_i > 30 or k_zd_i > 12 and downsampled to a grid 464² × 640. Note that OSIRIS simulations require to resolve the Debye scale, while the physical scales of interest are well represented at a lower resolution. A short-time average over ΔtΩ_ce = 2 (Ω_ce being the electron cyclotron frequency) was also performed to further reduce electron-scale noise (Grošelj et al., 2019). The CAMELIA data have been filtered for k_⊥d_i > 10 or k_zd_i > 2. We also considered alternative filtering approaches confirming that our results are not very sensitive to such particular choice.

3. Spectral Slopes and Normalized Field Ratios

Here we review and compare the standard set of spectral properties in our independently performed 3D kinetic simulations, namely the slopes of the turbulence power spectra and the spectral field ratios. Early theoretical predictions for sub-ion-range turbulence (e.g., Cho and Lazarian, 2004; Schekochihin et al., 2009) proposed a spectral scaling $~ k_{⊥}^{- 7 / 3}$ for the magnetic energy spectrum. However, solar wind observations typically exhibit much steeper magnetic spectra, namely $~ k_{⊥}^{- 2.8}$ (e.g., Alexandrova et al., 2009, 2013; Sahraoui et al., 2010; Chen, 2016; Kobayashi et al., 2017; Sorriso-Valvo et al., 2018). Similar spectral exponents were also reported in recent 3D kinetic simulations (Told et al., 2015; Cerri et al., 2017b, 2018; Franci et al., 2018a,b; Grošelj et al., 2018, 2019; Arzamasskiy et al., 2019). Recently, refined predictions were proposed to explain steeper spectra. Those include intermittency corrections (Boldyrev and Perez, 2012), dissipative effects (Howes et al., 2011; Passot and Sulem, 2015), and reconnection-mediated turbulence (Loureiro and Boldyrev, 2017; Mallet et al., 2017a). Further insight into the nature of kinetic-scale turbulence can be obtained from the spectral field ratios, which have been used to detect wave-like polarization properties in solar-wind turbulence and in kinetic simulations (e.g., Sahraoui et al., 2009; Salem et al., 2012; TenBarge et al., 2012; Chen et al., 2013; Kiyani et al., 2013; Cerri et al., 2017b; Franci et al., 2018b; Grošelj et al., 2018).

In Figure 1C the two-dimensional Fourier spectra, $ℰ (k_{⊥}, k_{z})$ , are shown. The wavenumber region (k_⊥, k_z) occupied by the turbulent fluctuations already highlights the anisotropic nature of the cascade, with energy preferentially flowing to high k_⊥. However, note that the 2D Fourier spectrum may exhibit a weaker anisotropy than the one typical of turbulent eddies, which are elongated along the local field direction (see, e.g., Cho and Vishniac, 2000). We perform a local analysis of anisotropy in section 4.1.

In Figures 1D,E, the reduced 1D spectra, $ℰ (k_{⊥})$ (upper panels), and their local slope (lower panels) are reported. To remove the effects of different energy injection conditions, the k_⊥-spectra have been normalized so that they overlap in the sub-ion range, at k_⊥d_i ≃ 5. According to the spectral anisotropy in Figure 1C, the k_z-spectra have been consequently matched at k_zd_i ≃ 0.5. For our choice of low-pass filter (see section 2), CAMELIA spectra artificially flatten beyond k_⊥d_i ≳ 7 due to PIC noise, and therefore we do not show CAMELIA data in this range in Figure 1D. Overall, the spectral slopes are consistent with each other, although the spectra obtained from the three simulations do not quite assume a universal shape. Close to the box scale, the spectral exponents are likely affected by the turbulence injection details. It is also possible that some of the sub-ion scale spectral exponents are not fully converged in terms of the box size (which was different for each simulation) and of the limited extent of sub-ion range itself. 3D3V simulations with a significantly larger sub-ion range are required to clarify this point. To some degree, differences at sub-ion scales could also be physical. In particular, the HVM simulation includes electron inertia effects in Ohm's law, while the OSIRIS results include the full electron kinetics, such as electron Landau damping and finite electron Larmor radius corrections. It is interesting to notice that OSIRIS spectra become steeper than the hybrid counterparts beyond k_⊥d_i ≳ 3, for our particular choice of the mass ratio (m_i/m_e = 100). This feature has been usually explained in terms of electron Landau damping (Grošelj et al., 2017; Chen et al., 2019), which is not included in the hybrid-kinetic model.

In Figure 1F, we report the comparison of spectral ratios, $C_{1} δ B_{z}^{2} / δ B_{⊥}^{2}$ (top), $C_{2} δ n^{2} / δ B_{⊥}^{2}$ (middle), and $C_{3} δ n^{2} / δ B_{z}^{2}$ (bottom). The ratios are normalized to the β-dependent kinetic Alfvén wave (KAW) eigenvalue from asymptotic linear theory ( $ρ_{i}^{- 1} ≪ k_{⊥} ≪ ρ_{e}^{- 1}$ and k_|| ≪ k_⊥), namely C₁ = (2 + β)/β, C₂ = (2 + β)β/4, and $C_{3} = β^{2} / 4$ , where β = β_i + β_e (see, e.g., Boldyrev et al., 2013, for details). In the normalized units, asymptotic KAW theory predicts a value of unity for all three ratios. This is essentially the result of KAWs developing a certain degree of magnetic compressibility at sub-ion scales, which sets the relation between δB_⊥ and δB_||, and requiring that compressive magnetic fluctuations are pressure balanced, which in turn provides a relation between δB_|| and δn (see, e.g., Schekochihin et al., 2009; Boldyrev et al., 2013). As found in previous studies (e.g., Salem et al., 2012; TenBarge et al., 2012; Chen et al., 2013; Cerri et al., 2017b; Franci et al., 2018b; Grošelj et al., 2018), the spectral field ratios are overall consistent with KAW-like turbulence at sub-ion scales. This is not completely surprising, as both in CAMELIA and OSIRIS simulations, Alfvénic fluctuations are injected. On the other hand, compressible magnetic fluctuations (i.e., including δB_||) are injected in HVM run, and yet KAW-like fluctuations still develop. It was also proposed that KAWs may, quite generally, emerge as a result of Alfvén waves interacting with large-scale inhomogeneities (Pucci et al., 2016). Thus, the KAW-like spectral properties at sub-ion scales appear to be a relatively robust feature, independent of the details of the turbulent fluctuations injected at the MHD scales (cf. Cerri et al., 2017a). While the results are overall consistent, some differences are also seen, most notably in the high-k_⊥ range (k_⊥d_i ≳ 10), which could be presumably attributed to various numerical artifacts. However, some deviations could also relate to differences between the hybrid-kinetic and fully kinetic model [for instance, some dispersion relation properties not being exactly the same (e.g., Told et al., 2016)].

So are the sub-ion-scale field polarizations indeed KAW-like? As discussed above, recent observations and kinetic simulations are consistent with such idea, although linear wave predictions are not necessarily satisfied precisely (e.g., Chen et al., 2013; Kiyani et al., 2013; Cerri et al., 2017b; Franci et al., 2018b). Chen et al. (2013) report an average value of 0.75 for the normalized ratio $C_{2} δ n^{2} / δ B_{⊥}^{2}$ , whereas (asymptotic) KAW theory predicts a value of unity. That latter may be due to different reasons, among which we remark the following two: (i) sub-ion-range turbulence is not made of purely KAW-like fluctuations, and/or (ii) the asymptotic conditions that are used in the derivation of linear theory predictions are not met exactly because of the limited sub-ion range of scales and/or because of the inherently non-linear dynamics of turbulence. These two explanations are not mutually exclusive, of course. Indeed, sub-ion-scale turbulence can in principle include contributions from wave-like fluctuations of other nature. This may include fluctuations consistent with whistler (e.g., Gary and Smith, 2009), ion-cyclotron (e.g., Omidi et al., 2014; Zhao et al., 2018), or ion Bernstein waves (e.g., Podesta, 2012; Del Sarto et al., 2017; Grošelj et al., 2017), to name a few. On the other hand, the spectral ratios could also deviate from linear KAW predictions as a result of non-linear dynamics. For example, Boldyrev et al. (2013) propose that, specifically the (normalized) $C_{2} δ n^{2} / δ B_{⊥}^{2}$ ratio may fall somewhat below the KAW prediction due to a (yet to be investigated) non-linear effect, analogous to the residual-energy phenomenon in MHD turbulence.

4. Multi-Point Structure Functions

Beyond energy spectra, fluctuations across different scales may be investigated in more detail via structure functions, i.e., the moments of local field increments (e.g., Frisch, 1995; Biskamp, 2008). Two-point structure functions, $S_{m}^{(2)}$ (m being the order), are most common. However, these cannot quantitatively produce the correct scaling for fluctuations with power spectra steeper than ~ k⁻³, assuming a clean power-law spectrum (Falcon et al., 2007; Cho and Lazarian, 2009). Therefore, structure functions using more than two points are generally required at kinetic scales. Essentially, higher-order increments yield a scale decomposition that is more effective in filtering out the large-scale fluctuations below k ≈ π/ℓ in spectral space, where ℓ is the increment scale. We also mention that if the signal is a polynomial of degree $N - 2$ , its corresponding 2nd-order, $N$ -point structure function vanishes (Cho, 2019). This makes multi-point structure functions more suitable for the analysis of relatively smooth signals with steep spectra (Schneider et al., 2004). A detailed review of $N$ -point increments, as well as their physical interpretation can be found in Cho (2019). Here, we consider for some field f(x) the conditional, five-point structure functions:

\begin{array}{l} S_{m}^{(5)} (ℓ, ϑ_{B_{loc}}) = {〈 | Δ f (x, ℓ) |^{m} | ℓ, ϑ_{B_{loc}} 〉}_{x} & (1) \end{array}

where $Δ f (x, ℓ) = [f (x + 2 ℓ) - 4 f (x + ℓ) + 6 f (x) - 4 f (x - ℓ) + f (x - 2 ℓ)] / \sqrt{35}$ is the (normalized) field increment, 〈… 〉_x is a space average, and ϑ_{B_loc} is the angle between the increment vector ℓ and the local mean magnetic field B_loc. The term “conditional” implies that S_m are defined as conditional averages of |Δf(x, ℓ)|^m, using only those points in the statistical sample that fall within a given (narrow) range for ℓ and ϑ_{B_loc}. We also considered three-point structure functions (see Figure 2A) and, for a limited number of cases, seven-point structure functions (not shown). Comparison between the three-point, five-point and seven-point structure functions shows not only qualitative similarities among the three cases, but an apparent quantitative convergence with increasing number of points. We chose to illustrate the results in Figure 2 in terms of five-point structure functions in order to provide better constraints for the theoretical predictions. Similar to two-point increments, where the local mean field is often defined as B_loc(x, ℓ) = [B(x) + B(x + ℓ)]/2 (e.g., Cho and Vishniac, 2000; Mallet et al., 2016), we obtain B_loc by averaging over the points used for the increment. For five-point increments, a reasonable definition is B_loc(x, ℓ) = [B(x + 2ℓ) + 4B(x + ℓ) + 6B(x) + 4B(x−ℓ) + B(x−2ℓ)]/16. It is straightforward to check that such mean field definition filters out fluctuations around the scale of the increment ~ ℓ, while preserving the contribution from scales larger than ℓ.

FIGURE 2

Figure 2. (A) Five-points, 2nd-order structure functions, $S_{2}^{(5)}$ , of δB_⊥ (top row) and δB_|| (bottom row) vs. ℓ_⊥ (continuous lines) and ℓ_|| (dashed lines). Here, || and ⊥ are defined with respect to the local field direction (i.e., δB_|| ≠ δB_z) with different angular tolerance, Δϑ (colored lines; see text for definition). $S_{2}^{(3)}$ with Δϑ = 1.5° are also shown for reference (gray lines). (B) Anisotropy scaling, ℓ_||(ℓ_⊥), of δB_⊥ (left panel), δB_|| (central panel), and δn (right panel), derived from $S_{2}^{(5)}$ with $Δ ϑ_{⊥} = 1 . 5^{°}$ nominal resolution. Three reference scalings are also shown. (C) Excess kurtosis, $K = S_{4} / {[S_{2}]}^{2} - 3$ vs. ℓ_⊥ for δB_⊥ (left panel), δB_|| (central panel), and δn (right panel). A 1/ℓ_⊥ scaling is given for reference.

In what follows, we consider field-perpendicular, $S_{m} (ℓ_{⊥}) \equiv S_{m}^{(5)} (ℓ_{⊥}, 9 0^{°} - Δ ϑ \leq ϑ_{B_{loc}} \leq 9 0^{°})$ , and field-parallel, $S_{m} (ℓ_{∥}) \equiv S_{m}^{(5)} (ℓ_{∥}, ϑ_{B_{loc}} \leq Δ ϑ)$ , five-point structure functions of the magnetic field and density fluctuations, where Δϑ represents a finite angular tolerance used in practice to determine the local perpendicular and parallel directions. We reduce Δϑ until the scalings appear converged. The field increments, from which we obtain the conditional structure functions, are evaluated at every grid point. In each grid point and at every scale, increments are sampled along random directions. The numbers of these random directions per grid point have been tested to provide a statistically significant (i.e., converged) sample. The sample that is used in the following is such that any structure function S_m(ℓ, ϑ_{B_loc}) counts at least 1.5 × 10⁵ points per scale ℓ, in any given band for ϑ_{B_loc}.

4.1. Spectral Anisotropy

A delicate point concerns the sub-ion-range spectral anisotropy, $k_{||} ~ {k_{⊥}}^{α}$ (cf., e.g., Schekochihin et al., 2009; Boldyrev and Perez, 2012; Cerri et al., 2018; Landi et al., 2019). As is known from MHD, electron-MHD (EMHD), and kinetic-reduced-MHD (KRMHD) turbulence (Cho and Vishniac, 2000; Cho and Lazarian, 2009; Meyrand et al., 2019), the true anisotropy is often revealed only when measured with respect to the local, scale-dependent mean magnetic-field direction. Somewhat contradicting estimates, obtained with different methods, for the sub-ion-scale anisotropy have been presented in recent works. Here, we revisit this issue using the above-mentioned implementation of five-point structure functions, consistently applied to all data.

In Figure 2A we show the perpendicular and parallel second-order structure function scalings, and in Figure 2B we show the inferred anisotropy, ℓ_||(ℓ_⊥). The characteristic parallel length ℓ_||(ℓ_⊥) at a given perpendicular scale ℓ_⊥ is obtained by finding the value of ℓ_||, at which the amplitudes of S₂(ℓ_||) and S₂(ℓ_⊥) match. To illustrate the sensitivity to the local mean field direction, we show in Figure 2A the convergence with respect to the angular tolerance Δϑ. The parallel scalings appear converged at Δϑ ≃ 3° for CAMELIA data and at around Δϑ ≃ 1.5° for HVM, whereas the OSIRIS results are somewhat less sensitive to Δϑ (converging already for Δϑ ≃ 6°). This difference may occur because OSIRIS simulation exhibits the weakest anisotropy (in absolute values). Physically, Δϑ should be approximately no larger than ~ ℓ_⊥/ℓ_|| of the small-scale turbulent eddies. Thus, smaller Δϑ are needed if a stronger anisotropy develops at the energy-containing scales.

All quantities seem to converge to a scaling close to $ℓ_{||} ~ {ℓ_{⊥}}^{2 / 3}$ (although δB_⊥ fluctuations in HVM exhibit a scaling closer to 1/3 over the range of scales across ℓ _⊥ ~ d_i(= ρ_i)). It is worth noticing, however, that this is not the end of the story, as the scaling is not quite 2/3 and additional effects such as B-field curvature may slightly change the anisotropy. Indeed, the field increments are taken along a straight line. If the local magnetic field lines are significantly curved over the extent of the increment stencil (= 4ℓ for five-point increments), the field increments will mix contributions from different field lines, in which case the anisotropy may be somewhat underestimated. It is worth mentioning that a scaling $ℓ_{||} ~ {ℓ_{⊥}}^{2 / 3}$ was proposed in Boldyrev and Perez (2012), based on a filling-factor correction for the fluctuation energy. Assuming the energy is concentrated in intermittent, two-dimensional structures as in Boldyrev and Perez (2012), the filling factor should scale as $k_{⊥}^{- 1} ~ l_{⊥}$ . The filling factor may be approximately estimated from the inverse scaling of the excess kurtosis (Matthaeus et al., 2015; see section 4.2). Our results shown in Figure 2C are indeed roughly consistent with an excess kurtosis scaling $~ ℓ_{⊥}^{- 1}$ , although this approximate scaling is overall better satisfied for δB_∥ and δn than for δB_⊥. Finally, we mention that an alternative anisotropy estimate, based on a spectral band-pass filter (Cho and Lazarian, 2009), gives a somewhat stronger anisotropy than the five-point structure functions (not shown). On the other hand, qualitatively similar results are still obtained for all data. Thus, all simulations analyzed exhibit a similar sub-ion-scale anisotropy according to the particular diagnostics employed. Therefore, the differences that were previously reported in the literature could be mainly related to the different methods employed.

4.2. Intermittency: The “Saturation Problem”

Another relevant feature of kinetic plasma turbulence is the excess kurtosis of the fluctuations, $K (ℓ_{⊥}) = S_{4} (ℓ) / {[S_{2} (ℓ)]}^{2} - 3$ . The increase of K(ℓ_⊥) above zero is a measure of non-Gaussian statistics of the turbulent fluctuations (Frisch, 1995; Matthaeus et al., 2015). As seen in Figure 2C, the excess kurtosis gradually increases above the Gaussian value throughout the sub-ion scale range. Moreover, similar statistical trends are seen for δB_⊥, δB_∥, and δn [note that we take here the component of δB_⊥ parallel to ℓ × B_loc to estimate the flatness of δB_⊥ (see also Kiyani et al., 2013)]. In apparent contrast with our results, a number of observational studies of solar wind turbulence find non-Gaussian, yet nearly scale-independent turbulence statistics at sub-ion scales (Kiyani et al., 2009, 2013; Wu et al., 2013; Chen et al., 2014). Thus, it appears a process operates in the solar wind that saturates the turbulence statistics already near the transition to sub-ion scales (ℓ_⊥ ≲ d_i). What could be the reason for this apparent contradiction? One clear difference is that the solar-wind fluctuations are already heavily non-Gaussian at MHD scales (Salem et al., 2009), whereas our 3D kinetic simulations do not quite share the same feature due to the limited simulation domain. We mention that even large-size 2D kinetic simulations (e.g., Wan et al., 2012; Franci et al., 2015; Leonardis et al., 2016) did not yet achieve K(ℓ_⊥) ≫ 1 in the MHD range (ℓ_⊥ ≫ d_i). In this context, it may be worth pointing out that intermittency in MHD turbulence is commonly associated with the emergence of sheetlike structures (e.g., Chandran et al., 2015; Matthaeus et al., 2015; Mallet and Schekochihin, 2017), which may break apart via the tearing instability (causing the field lines to reconnect), once their perpendicular aspect ratio exceeds a critical threshold (Matthaeus and Lamkin, 1986; Boldyrev and Loureiro, 2017; Mallet et al., 2017b). For sub-ion-scale turbulence, the possible role of magnetic reconnection has been as well highlighted in a number of recent works (e.g., Franci et al., 2016, 2017; Cerri and Califano, 2017; Loureiro and Boldyrev, 2017; Mallet et al., 2017a; Papini et al., 2019). Moreover, a recent observational study (Vech et al., 2018) argued that the spectral break at the tail of the MHD cascade may be controlled by reconnection. Therefore, the phenomenology of the cascade may critically depend on the morphology of the intermittent structures at the transition into the kinetic range (Mallet et al., 2017a). If the structures are indeed sufficiently sheetlike to be tearing unstable, collisionless reconnection might be one possible process that limits the growth of the sub-ion scale kurtosis (see also Biskamp et al., 1990; Chen et al., 2014). However, alternative possibilities such as collisionless damping of the fluctuations cannot be ruled out at this stage.

5. Concluding Remarks

So, what is the nature of sub-ion-scale fluctuations? From our independently performed 3D3V (hybrid and fully) kinetic simulations, a picture consistent with KAW turbulence phenomenology emerges. Moreover, our results imply a scale-dependent anisotropy, together with intermittent statistics of magnetic and density fluctuations at sub-ion scales. Thus, we conclude that within the range of parameters explored here, the statistical properties of ion-scale plasma turbulence (at β ~ 1) definitely show a certain degree of similarity, regardless of the precise details of the large-scale energy injection. On the other hand, slight differences can also be identified, some of which may be also model-dependent.

A number of key aspects will have to await the next-generation of 3D3V kinetic simulations. Ideally, future numerical experiments should aim to resolve both larger (MHD) scales, as well as a broader range between the ion and the electron scales by adopting significantly higher (if not realistic) mass ratios. These two aspects indeed appear to be both required in order to achieve (i) a possible saturation of the kurtosis at ion scales and (ii) a relevant sub-ion range of scales before electron-scale effects significantly come into play. Moreover, different aspects other than the spectral and statistical properties of the turbulent fluctuations will need to be considered in characterizing kinetic-range turbulence, as for instance, the dissipation mechanisms of turbulent fluctuations under different plasma conditions and the consequent energy partition among different species (e.g., Matthaeus et al., 2016; Parashar et al., 2018; Arzamasskiy et al., 2019; Kawazura et al., 2019; Zhdankin et al., 2019).

While certain progress was definitely achieved in recent years, many other plasma regimes and setups may need to be explored, and the process(es) underlying a possible universality of kinetic-range plasma turbulence (e.g., magnetic reconnection) need to be fully worked out. Moreover, a few relevant discrepancies between numerical simulations, theories and in-situ observations appear. These “anomalies” definitely call for an explanation by the space physics community. In this context, advances cannot be achieved without investing in next-generation multi-spacecraft missions. Multi-point in situ measurements of turbulent fluctuations from a large number of spacecrafts are indeed fundamental in order to disentangle the non-linear spatio-temporal character of plasma turbulence (see, e.g., Klein et al., 2019; Matthaeus et al., 2019; TenBarge et al., 2019). This includes answering fundamental questions about, for instance, (i) the distribution of turbulent energy in space and time, (ii) the three-dimensional anisotropic structure of energy transfer across scales, (iii) the high-order statistics of the fluctuations, and (iv) the validity of Taylor's hypothesis over a broad range of time and spatial scales. Alongside observations, advances in computational capabilities are required to perform more realistic numerical simulations as discussed above, and compare these with spacecraft measurements. Finally, following the same spirit promoted by the “Turbulent Dissipation Challenge” (Parashar et al., 2015), we would like to end this Perspective by stressing that our community could benefit from comparisons such as the one performed here, involving various codes, models and diagnostics.

Note added: Arzamasskiy et al. (2019) recently reported a scale-independent anisotropy at ion scales (i.e., ℓ_∥ ~ ℓ_⊥) based on a set of 3D driven hybrid-kinetic turbulence simulations. Using our structure function diagnostic applied to their data, we were able to qualitatively (and quantitatively) reproduce their result. A more detailed investigation along these lines is currently ongoing, but beyond the scope of this Perspective and will be presented elsewhere.

Data Availability Statement

The data used in this study are available from the authors upon reasonable request.

Author Contributions

SC, DG, and LF provided their HVM, OSIRIS, and CAMELIA simulation data, respectively. SC performed the spectral analysis, produced the figures, and wrote the paper, taking into account suggestions from DG and LF. DG implemented and performed the structure function analysis. All authors discussed the results.

Funding

SC was supported by the National Aeronautics and Space Administration under Grant No. NNX16AK09G issued through the Heliophysics Supporting Research Program. LF was supported by the UK Science and Technology Facilities Council (STFC) grant ST/P000622/1.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

SC and DG acknowledge the generous hospitality of the Wolfgang Pauli Institute in Vienna, where the first discussions leading to this work took place. We acknowledge PRACE for awarding us access to Marconi at CINECA, Italy, where the calculations with the HVM code were performed under the grant No. 2017174107. The Cray XC40, Shaheen, at the King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia was utilized for the simulation performed with the OSIRIS code. LF acknowledges PRACE for awarding him access to Cartesius at SURFsara, the Netherlands, through the DECI-13 (Distributed European Computing Initiative) call (project HybTurb3D) where the HPIC simulation was performed, and INAF and CINECA for awarding him access to Marconi within the framework of the MoU New Frontiers in Astrophysics: HPC and New Generation Data Exploration (project INA17_C4A26), where new further analysis of the HPIC data has been performed. The authors would like to acknowledge the OSIRIS Consortium, consisting of UCLA and IST (Lisbon, Portugal) for the use of OSIRIS 3.0 and for providing access to the OSIRIS 3.0 framework. SC acknowledges Dr. C. Cavazzoni and Dr. M. Guarrasi (CINECA, Italy) for their contributions to HVM code parallelization, performance and implementation on Marconi-KNL. The authors also acknowledge useful discussions with Alfred Mallet, Lev Arzamasskiy, Bill Dorland, Matt Kunz, Simone Landi, Emanuele Papini, Frank Jenko, and David Burgess.

References

Adkins, T., and Schekochihin, A. A. (2018). A solvable model of Vlasov-kinetic plasma turbulence in FourierHermite phase space. J. Plasma Phys. 84:905840107. doi: 10.1017/S0022377818000089