2.1. Highlights of past works on GK turbulence

In classical fluid turbulence (Frisch Reference Frisch1995), the energy cascade, the locality of interactions and the intermittency behaviour are considered standard problems of interest. Although turbulence at kinetic scales inherits all of them, it also adds the phase space mixing problem (includes Landau damping) that affects which route in phase space is selected for the thermalisation of plasma fluctuations. In magnetised plasma, all of these problems can be tackled via the GK approximation (for a review on the formal derivation of the general equations, see Brizard & Hahm Reference Brizard and Hahm2007).

The GK approximation was instrumental in probing turbulence at sub-ion scales ($k\rho _i>1$). GK theory assumes low plasma frequencies compared with the ion cyclotron frequency and small fluctuation levels compared with background quantities to remove the particles’ fast gyro-motion, effectively reducing the relevant phase space to five-dimensions. This approach was adopted by Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006) for the study of KAWs and their turbulent cascade in the dissipative range of the solar wind (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008a, Reference Howes, Dorland, Cowley, Hammett, Quataert, Schekochihin and Tatsunob, Reference Howes, Tenbarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011). Compared with the use of global background profiles in tokamak geometries (see Krommes Reference Krommes2012), the use of local background approximations in a straight-field magnetic geometry, typical for the study of KAW turbulence, simplifies the underlying dynamics.

Following the recipe of classical turbulence, a generalised free energy that is conserved in the absence of collisions was identified for GK turbulence (see Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008, Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). With the idea of a free energy cascade in phase space, the concept of the nonlinear phase mixing for GK was introduced as well (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008, Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009). The nonlinear phase mixing occurs in the direction perpendicular to the magnetic guide field, and it refers in particular to the creation of small-scale structures in velocity space owing to the small-scale structure in position space. This effect results from the nonlinear interaction between the distribution function and the gyroaveraged potential fields. The gyroaverage represents the effect of the fast gyro-motion on the slower dynamics captured by GK theory. In the electrostatic limit, the phase space cascade and the nonlinear phase mixing were studied extensively (Tatsuno et al. Reference Tatsuno, Dorland, Schekochihin, Plunk, Barnes, Cowley and Howes2009, Reference Tatsuno, Barnes, Cowley, Dorland, Howes, Numata, Plunk and Schekochihin2010; Plunk & Tatsuno Reference Plunk and Tatsuno2011; Tatsuno et al. Reference Tatsuno, Plunk, Barnes, Dorland, Howes and Numata2012) for ‘two-dimensional’ GK turbulence (i.e. neglecting parallel dynamics, see Plunk et al. Reference Plunk, Cowley, Schekochihin and Tatsuno2010). For the five-dimensional GK system, while still in the electrostatic limit, the energy balance equation and the energy cascade problem was studied by Navarro et al. (Reference Navarro, Morel, Albrecht-Marc, Carati, Merz, Görler and Jenko2011a, Reference Navarro, Morel, Albrecht-Marc, Carati, Merz, Görler and Jenkob), Nakata, Watanabe & Sugama (Reference Nakata, Watanabe and Sugama2012) and later by Teaca, Navarro & Jenko (Reference Teaca, Navarro and Jenko2014), Cerri et al. (Reference Cerri, Navarro, Jenko and Told2014) and Maeyama et al. (Reference Maeyama, Idomura, Watanabe, Nakata, Yagi, Miyato, Ishizawa and Nunami2015). Measuring the intensity of the energetic exchanges with the increase in separation between scales, the locality of the nonlinear interactions was studied for electrostatic GK turbulence in Teaca et al. (Reference Teaca, Navarro, Jenko, Brunner and Villard2012, Reference Teaca, Navarro and Jenko2014) and for the electromagnetic KAW case in Told et al. (Reference Told, Jenko, Tenbarge, Howes and Hammett2015) and Teaca, Jenko & Told (Reference Teaca, Jenko and Told2017). Although GK turbulence exhibits a strong nonlocal interaction character, Teaca et al. (Reference Teaca, Jenko and Told2017) found that the nonlocal contribution is superimposed on top of a classic asymptotically local contribution that depends only with the separation between scales, rather than substituting the classic local character altogether. This is encouraging when considering modelling the SGS effects. Last, the intermittency problem was looked at in phase space for KAW turbulence by Teaca et al. (Reference Teaca, Navarro, Told, Görler, Plunk, Hatch and Jenko2019), where the deviation from scale invariance was measured directly on the distribution functions.

In relation to the dissipation route for magnetised plasma fluctuations, Told et al. (Reference Told, Jenko, Tenbarge, Howes and Hammett2015) showed via a multi-species GK simulation of KAW turbulence at plasma $\beta =1$ that electrons dissipate most of the free energy at ion scales ($k_{\perp }\rho _i \sim 1$), whereas ions dissipate at small scales ($k_{\perp } \rho _i>1$). Later, Navarro et al. (Reference Navarro, Teaca, Told, Groselj, Crandall and Jenko2016) showed on the same data that the electrons prefer parallel collisions, indicative of parallel linear phase mixing (Hammett, Dorland & Perkins Reference Hammett, Dorland and Perkins1992; Kanekar et al. Reference Kanekar, Schekochihin, Dorland and Loureiro2015), whereas ions enter into a fluid-like cascade in the perpendicular direction. The linear phase mixing problem is tied to the Landau damping problem for GK turbulence (Tenbarge & Howes Reference Tenbarge and Howes2013) and has a non-trivial effect on its structure character (Teaca et al. Reference Teaca, Navarro, Told, Görler, Plunk, Hatch and Jenko2019). The balance between linear phase mixing in the parallel direction and the nonlinear cascade in the perpendicular direction was introduced for a drift-kinetic reduced model in Schekochihin et al. (Reference Schekochihin, Parker, Highcock, Dellar, Dorland and Hammett2016). The four-dimensional drift-kinetic models in question integrate over the perpendicular velocity, while retaining the information for $k\rho _i<1$ scale dynamics (see also Hatch et al. Reference Hatch, Jenko, Bratanov and Navarro2014). The use of these models was helpful in showcasing the linear flux of energy across parallel velocity scales induced by linear phase mixing, and its suppression that leads to the fluidisation of the kinetic turbulent problem (Meyrand et al. Reference Meyrand, Kanekar, Dorland and Schekochihin2019). Finally, depending on the plasma parameters (plasma-$\beta$, in particular), Kawazura, Barnes & Schekochihin (Reference Kawazura, Barnes and Schekochihin2019) showed via hybrid GK simulations that ions can exhibit parallel or (fluid-like) perpendicular dissipation routes in phase space.

Although understanding turbulence is a goal in itself, in tokamak studies, turbulence is seen as a problem that overcomplicates the study of heat and particle transport by energising small-scale fluctuations compared with the scale of the dominant linear instabilities (Görler & Jenko Reference Görler and Jenko2008a,Reference Görler and Jenkob). To model the effect of these small scales on the nonlinear interactions at large scales, LES have been adopted for GK turbulence (Morel et al. Reference Morel, Navarro, Albrecht-Marc, Carati, Merz, Görler and Jenko2011, Reference Morel, Navarro, Albrecht-Marc, Carati, Merz, Görler and Jenko2012), and were refined further in Navarro et al. (Reference Navarro, Teaca, Jenko, Hammett and Happel2014). To put it simply, LES models SGS effects. Although known extensively in the field of turbulence (see Eyink & Sreenivasan Reference Eyink and Sreenivasan2006, and the references within), an SGS analysis for kinetic turbulence was introduced by Eyink (Reference Eyink2018), where the entropy cascade was rigorously defined for a full Vlasov–Maxwell–Landau system and an upper bound scaling computed via functional analysis. Considering that velocity space integrals are performed in addition to position space ones, cancellation effects cannot be overlooked when computing the actual fluxes across scales. To what degree the upper bound estimates overshoot the real levels can only be determined numerically, and it is one of the questions we answer in the current paper for the GK system.

2.2. Plasma parameters and numerical simulation details

Depending on the geometry of the external magnetic guide field and the plasma regime, the GK equations can have an intricate or simple explicit form. Before introducing the GK equations, we start by presenting the main parameters for the plasma considered and list the numerical details used to solve the system in practice.

In this study, we look at a proton–electron plasma that is weakly collisional and strongly magnetised, and which evolves in the presence of a straight magnetic guide field ($B_0\boldsymbol {\hat {z}}$). Proton (referred to as ion) and electron species are included with their real mass ratio of $m_{i}/m_{e}=1836$. The plasma $\beta _{i}\equiv 8{\rm \pi} n_{i}T_{i}/B_{0}^{2}=1$ is chosen to match solar wind conditions at 1 astronomical unit. The plasma background is assumed to exhibit an isotropic thermodynamic equilibrium with a temperature ratio of $T_{i}/T_{e}=1$. The electron collisionality is chosen to be $\nu _{e}=0.06\, \omega _{A0}$ (with $\nu _{i}=\sqrt {m_{e}/m_{i}}\nu _{e}$), and $\omega _{A0}$ being the frequency of the slowest Alfvén wave in the system. This allows for a KAW cascade.

The system is solved numerically with the help of the Eulerian code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000). The data used in this work is from the simulation presented in Told et al. (Reference Told, Jenko, Tenbarge, Howes and Hammett2015), and it is briefly summarised in the following: the evolution of the gyrocentre distribution is tracked on a grid with the resolution $\{N_{x}, N_{y}, N_{z}, N_{v_{\parallel }}, N_{\mu },N_{s}\}=\{768, 768, 96, 48, 15, 2\}$, where ($N_x,N_y$) are the perpendicular, $(N_z)$ parallel, $(N_{v_{\parallel }})$ parallel velocity and $(N_{\mu })$ magnetic moment grid points, respectively. This covers a perpendicular dealiased wavenumber range of $0.2\le k_{\perp }\rho _{i}\leq 51.2$ (or $0.0047\leq k_{\perp }\rho _{e}\leq 1.19)$ in a domain $L_x=L_y=10{\rm \pi} \rho _i$ ($\rho _s= \sqrt {T_{s} m_{s}} c / e B$). In the parallel direction, a $L_z=2{\rm \pi} L_{\parallel }$ domain is used, where $L_{\parallel } \gg \rho _i$ is assumed by the construction of GK theory. A velocity domain up to three thermal velocity units ($v_{T,s}=\sqrt {2T_s/m_s}$) is taken in each direction. The fluctuations in the system are driven to a steady state via a magnetic antenna potential, which is prescribed solely at the largest scale and evolved in time according to a Langevin equation (TenBarge et al. Reference TenBarge, Howes, Dorland and Hammett2014).

2.3. The GK equations

For the system considered previously, we use the $\delta f$-approach. The particle distribution function of each plasma species $s$ is split into a time constant background $F_{s}$ and a perturbed part $\delta f_{s}$, with $\delta f_{s}/F_s \ll 1$. We consider a local approximation for $F_{s}$, which for constant background density $n_s$ and temperature $T_s$ (again, $v_{T,s}=\sqrt {2T_s/m_s}$) has the Maxwellian form,

(2.1)\begin{equation} F_{s}({\boldsymbol{v}})=\frac{n_s}{( v_{T,s}\sqrt{{\rm \pi}})^{3}} \exp\left(-\frac{v_{\|}^{2}+v_{{\perp}}^{2}}{v_{T,s}^{2}}\right). \end{equation}

In the presence of a strong external magnetic field compared to the fluctuating electromagnetic fields, the dynamics of the plasma become strongly anisotropic ($k_{\|}/k_{\perp }\ll 1$). More importantly, particles develop fast cyclotron motions (of gyro-frequency $\varOmega _s={q_sB_0}/{m_s c}$) compared with the rest of the plasma dynamics ($\omega /\varOmega _s\ll 1$). Employing the guiding centre coordinate ($\boldsymbol {R}_s=x\boldsymbol {\hat {x}}+y\boldsymbol {\hat {y}}+z\boldsymbol {\hat {z}}$) transformation

(2.2)\begin{equation} \boldsymbol{R}_s=\boldsymbol{r} + \boldsymbol{v}(\theta)\times\boldsymbol{\hat{z}}/\varOmega_s=\boldsymbol{r} + \boldsymbol{v}_{{\perp}}(\theta)\times\boldsymbol{\hat{z}}/\varOmega_s, \end{equation}

for $\boldsymbol {v}(\theta )=\boldsymbol {v}_{\perp }(\theta )+v_{\|}\boldsymbol {\hat {z}} =v_{\perp } \sin (\theta )\boldsymbol {\hat {x}}+v_{\perp } \cos (\theta )\boldsymbol {\hat {y}}+v_{\|}\boldsymbol {\hat {z}}$, and integrating the dynamics over the gyrophase angle ($\theta$) allows us to reduce the dimension of the phase space by one, obtaining the five-dimensional gyrocentre phase space ($\boldsymbol {R}_s, v_{\|}, v_{\perp }$) of GK theory. We can substitute the perpendicular velocity with the magnetic moment $\mu =m_s v^{2}_{\perp }/2B_0$. Although we do this in practice, some relations are more transparent when utilising $v_{\perp }$.

For this simple case, considering $\delta f_{s}/F_s \sim B_{\perp }/B_0 \sim B_{\|}/B_0 \sim (cE_{\perp }/v_{T,s})/B_0 \sim k_{\|}/k_{\perp } \sim \omega /\varOmega _s \sim \epsilon \ll 1$ as the the GK ordering, expanding all fields in powers of $\epsilon$ and keeping contributions up to the first order, the perturbed distribution function becomesFootnote ¹

(2.3)\begin{equation} \delta f_s(\boldsymbol{r}, \boldsymbol{v},t) ={-} \frac{q_s \phi(\boldsymbol{r},t) }{T_s}F_s(v) + h_s(\boldsymbol{R}_s,v_{\|},v_{{\perp}},t) , \end{equation}

where we see a Boltzmann response contribution and a non-adiabatic part, $h_s(\boldsymbol {R}_s,v_{\|},v_{\perp },t)$, which here is the effective GK distribution function.

The systematic expansion of the Vlasov–Maxwell system gives rise to the GK equations (see Brizard & Hahm (Reference Brizard and Hahm2007) for a general Hamiltonian derivation, or Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006), Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) for a simpler presentation appropriate in our case). For the first-order contribution $h_{s}(x,y,z,v_{\parallel },\mu ,t)$, the GK equations have the form

(2.4)\begin{equation} \frac{\partial h_s}{\partial t}+ \frac{c}{B_0}\{ \langle \chi \rangle_{\boldsymbol{R}_s},h_s\} + v_{{\parallel}}\frac{\partial h_s}{\partial z}=\frac{q_s F_{s}}{T_{s}}\frac{\partial\langle \chi \rangle_{\boldsymbol{R}_s}}{\partial t}+\left(\frac{{\partial} h_s}{\partial t}\right)_c. \end{equation}

Although the electromagnetic potentials are computed at the particle position ($\boldsymbol {r}$), only their gyroaveraged contributions affect the GK dynamics. For clarity, the gyroaveraged GK potential $\langle \chi (\boldsymbol {r}) \rangle _{\boldsymbol {R}_s} = ({1}/{2{\rm \pi} })\int _0^{2{\rm \pi} } \chi (\boldsymbol {R}_s-{\boldsymbol {v}_{\perp }(\theta )\times \boldsymbol {\hat {z}}}/{\varOmega _s} )\,\textrm {d}\theta$ is found via its wave-space representation as

(2.5)\begin{equation} \langle \chi \rangle_{\boldsymbol{R}_s}= \sum_{\boldsymbol{k}} \exp({\textrm{i} \boldsymbol{k} \cdot \boldsymbol{R}}) \left[ \textrm{J}_0(\lambda_s) \left(\hat\phi(\boldsymbol{k}) -\frac{v_{{\parallel}}}{c} \hat{A}_{{\parallel}}{(\boldsymbol{k})} \right)+\frac{2\mu}{q_s}\frac{\textrm{J}_1(\lambda_s)}{\lambda_s} \hat{B}_{{\parallel}}(\boldsymbol{k})\right], \end{equation}

where $\textrm {J}_0(\lambda _s)$ and $\textrm {J}_1(\lambda _s)$ are zero- and first-order Bessel functions, with $\lambda _s=k_{\perp } v_{\perp }/\varOmega _s=k_{\perp }\sqrt { 2\mu B_0/m_s}/\varOmega _s$. The first-order self-consistent electrostatic potential ($\phi$), magnetic potential in the parallel direction ($A_{\parallel }$) and magnetic fluctuation in the parallel direction ($B_{\parallel }$) are obtained in wave space from their respective GK field equations as

(2.6)\begin{gather} \hat{\phi}(\boldsymbol{k},t)=\sum_s 2{\rm \pi} q_s\frac{ B_0}{m_s}\int_{-\infty}^{+\infty}\,\textrm{d}v_{\|}\int_{0}^{+\infty}\,\textrm{d}\mu\, \textrm{J}_0(\lambda_s) \hat{h}_s(\boldsymbol{k}, v_{{\parallel}},\mu,t) \left/ \sum_s \frac{q_s^{2} n_s}{T_s}\right., \end{gather}

(2.7)\begin{gather}\hat{A}_{\|}(\boldsymbol{k},t)=\frac{4{\rm \pi}}{k_{{\perp}}^{2} c} \sum_s 2{\rm \pi} q_s\frac{ B_0}{m_s}\int_{-\infty}^{+\infty}\,\textrm{d}v_{\|}\int_{0}^{+\infty}\,\textrm{d}\mu\, v_{\|} \textrm{J}_0(\lambda_s) \hat{h}_s(\boldsymbol{k}, v_{{\parallel}},\mu,t), \end{gather}

(2.8)\begin{gather}\hat{B}_{\|}(\boldsymbol{k},t)={-}{4{\rm \pi}} \sum_s 2{\rm \pi} \frac{B_0}{m_s}\int_{-\infty}^{+\infty}\,\textrm{d}v_{\|}\int_{0}^{+\infty}\,\textrm{d}\mu\, \mu \frac{2\textrm{J}_1(\lambda_s)}{\lambda_s} \hat{h}_s(\boldsymbol{k}, v_{{\parallel}},\mu,t) . \end{gather}

Considering the values of the Bessel functions $\textrm {J}_0$ and $2\textrm {J}_1(\lambda _s)/\lambda _s$, we see that the gyroaverage operation cannot be ignored for $k_{\perp } \rho _i >1$ scales, where we can imagine the problem as the distribution of a system of electrical charged rings. Conversely, the gyroaverage operation is not that important for $k_{\perp } \rho _i <1$ scales, and drift-kinetic approximations can be obtained in the $k_{\perp } \rho _i \ll 1$ limit, which can still account for gyroaverage effects in a simplified way (Hammett et al. Reference Hammett, Dorland and Perkins1992; Hatch et al. Reference Hatch, Jenko, Bratanov and Navarro2014).

The ${(}{\partial h_s}/{\partial t}{)}_c$ term represents the action of collisions, which are here modelled through the action of a linearised Landau–Boltzmann collision operator (see supplementary material from Navarro et al. Reference Navarro, Teaca, Told, Groselj, Crandall and Jenko2016). Collisions represent the ultimate sink of plasma fluctuations and, in the collisionless limit, they are assumed to occur at very small scales in velocity space. For GK theory, owing to the nonlinear phase mixing, the small scales in the perpendicular velocity and the perpendicular small scales in position space are linked. As a result, dissipation in the perpendicular direction occurs similarly as for a fluid via an effective (hyper) Laplacian term in position space. For GK turbulence, the break from the fluidisation can occur only when the parallel collisions dominate (Navarro et al. Reference Navarro, Teaca, Told, Groselj, Crandall and Jenko2016) and higher-velocity moments in the $v_{\|}$ direction become excited via linear phase mixing.

The nonlinear structure is given in terms of the spatial Poisson bracket (to simplify the notation of gradients, from now on $\boldsymbol {\nabla }\equiv \boldsymbol {\nabla }_{\boldsymbol {R}_s} ={\partial }/{\partial \boldsymbol {R}_s}$),

(2.9)\begin{equation} \{ a, b\}= [\boldsymbol{\nabla} a \times \boldsymbol{\nabla} b] \boldsymbol{\cdot} {\boldsymbol{\hat{z}}}=\frac{\partial a}{\partial x}\frac{\partial b}{\partial y}-\frac{\partial a}{\partial y}\frac{\partial b}{\partial x}, \end{equation}

which possesses all its properties (antisymmetry, bilinearity, etc., see Appendix A for details). To highlight the advective role of the nonlinearity, we can rewrite it as

(2.10)\begin{equation} \frac{c}{B_0}\{ \langle \chi \rangle_{\boldsymbol{R}_s},h_s\}= {\boldsymbol{u}}_s\boldsymbol{\cdot} \boldsymbol{\nabla} h_s=\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{u}}_s h_s) , \end{equation}

where the advective velocity ${\boldsymbol {u}}_s$ is simply the generalised drift velocity for GK,

(2.11)\begin{equation} {\boldsymbol{u}}_s={-}\frac{c}{B_0}[\boldsymbol{\nabla} \langle \chi \rangle_{\boldsymbol{R}_s} \times {\boldsymbol{\hat{z}}} ]. \end{equation}

By definition, it is clear that ${\boldsymbol {u}}_s={\boldsymbol {u}}_s(x,y,z,v_{\|},\mu ,t)$ is zero-divergent ($\boldsymbol {\nabla }\boldsymbol {\cdot } {\boldsymbol {u}}_s=0$) and that it differs slightly for each species due to the gyroaverage, see figure 1. For the analysis of the nonlinear redistribution of free energy, we will utilise the advective velocity form for the nonlinear term. Although key results will also be presented in Poisson bracket form, the advective velocity form allows for a much simpler connection with classical turbulence. As mentioned in § 2.1, the analogue to the energy cascade in classical turbulence is given for GK turbulence by the free energy cascade.

Figure 1. Generalised drift velocity ${\boldsymbol {u}}_s$ for the GK system, at one point in $z$, $v_{\|}$ and $\mu$ ($v_{\|}=-0.31$ and $\mu =0.015$ in thermal velocity units). The difference between the ion and electron species consists in the gyroaverage of $\chi$. The electrons (a), which have a small gyroaverage that can be neglected for the scales depicted, show classical turbulent structures. For the ions (b), we see the phase mixing caused by the gyroaverage.

2.4. The free energy

As presented in Howes et al. (Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006), Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008, Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), the generalised free energy is conserved for GK turbulence in the absence of collisions and external sources. The free energy is defined as,

(2.12)\begin{equation} W=\int\,\textrm{d}^{3}r\left[ \sum_s \int\,\textrm{d}^{3}v \frac{T_s \delta f_s^{2} }{2F_s} +\frac{B^{2}}{8{\rm \pi}}\right],\end{equation}

where we neglect the electric field energy contribution due to free charges, as the scales of interest here are much larger than the Debye length. Considering the quantities that express the GK equation and (2.3), the equivalent definitions are obtained,

(2.13)\begin{align} W&=\int\,\textrm{d}^{3}r\left[ \sum_s\int\,\textrm{d}^{3}v \frac{T_s \langle h_s^{2}\rangle_{\boldsymbol{r}} }{2F_s} - \sum_s\frac{q^{2}_s n_s}{2T_s}\phi^{2} +\frac{{|\boldsymbol{\nabla}_{{\perp}} A_{\|}\,|}^{2}}{8{\rm \pi}} +\frac{B_{\|}^{2}}{8{\rm \pi}}\right] \nonumber\\ &=\sum_{\boldsymbol{k}} \left[ \sum_s \frac{2{\rm \pi} B_0}{m_s}\int \,\textrm{d}v_{\|}\,\textrm{d}\mu \frac{T_s }{2F_s}|\hat{h}_s(\boldsymbol{k},v_{\|},\mu,t)|^{2} \right.\nonumber\\ &\quad \left.-\, \sum_s \frac{q^{2}_s n_s}{2T_s}|\hat{\phi}(\boldsymbol{k})|^{2} +\frac{k^{2}_{{\perp}} |\hat{A}_{\|}(\boldsymbol{k}) |^{2}}{8{\rm \pi}}+\frac{|\hat{B}_{\|}(\boldsymbol{k})|^{2}}{8{\rm \pi}}\right] , \end{align}

with

(2.14)\begin{align} \langle h^{2} \rangle_{\boldsymbol{r}} &= \frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}} h^{2}(\boldsymbol{r}+{\boldsymbol{v}_{{\perp}}(\theta)\times\boldsymbol{\hat{z}}}/{\varOmega_s}, v_{\|}, \boldsymbol{v}_{{\perp}}(\theta) ) \,\textrm{d}\theta \nonumber\\ &=\sum_{\boldsymbol{k}} \exp({\textrm{i} \boldsymbol{k} \cdot \boldsymbol{r}}) \textrm{J}_0\left(\frac{k_{{\perp}} v_{{\perp}}}{\varOmega_s}\right)\hat{h}^{2}(\boldsymbol{k}, v_{\|}, v_{{\perp}}). \end{align}

Considering the contribution of individual terms, with an appropriate selective summation in wave space defined as ${\sum _{\perp }\equiv \sum _{k_z} \sum _{|\boldsymbol {k}_{\perp }| = k_{\perp }}^ {k_{\perp } + \Delta k} }$, we compute the unit band ($\Delta k$) spectra in the perpendicular direction

(2.15)\begin{gather} W_{h_s}(k_{{\perp}}) = \sum_{{\perp}} \frac{2{\rm \pi} B_0}{m_s}\int \,\textrm{d}v_{\|}\,\textrm{d}\mu \frac{T_s }{2F_s}|\hat{h}_s(\boldsymbol{k},v_{\|},\mu,t)|^{2}, \end{gather}

(2.16)\begin{gather}W_{\phi}(k_{{\perp}}) = \sum_{{\perp}} \sum_s \frac{q^{2}_s n_s}{2T_s}|\hat{\phi}(\boldsymbol{k})|^{2}, \end{gather}

(2.17)\begin{gather}W_{B_{{\perp}}}(k_{{\perp}}) = \sum_{{\perp}} \frac{k^{2}_{{\perp}} |\hat{A}_{\|}(\boldsymbol{k}) |^{2}}{8{\rm \pi}}, \end{gather}

(2.18)\begin{gather}W_{B_{\|}}(k_{{\perp}})= \sum_{{\perp}} \frac{|\hat{B}_{\|}(\boldsymbol{k})|^{2}}{8{\rm \pi}} . \end{gather}

The total free energy spectrum can be found simply as the sum,

(2.19)\begin{equation} W(k_{{\perp}})=\sum_s W_{h_s}(k_{{\perp}}) - W_{\phi}(k_{{\perp}}) + W_{B_{{\perp}}}(k_{{\perp}}) + W_{B_{\|}}(k_{{\perp}}). \end{equation}

We plot in figure 2 the spectra for all the contributions to the free energy. We see that the so-called (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008) entropic contributions ($W_{h_s}$) dominate the free energy. The scaling of the magnetic fields is the same as listed in Told et al. (Reference Told, Cookmeyer, Muller, Astfalk and Jenko2016). Notably, for $k_{\perp } \rho _i<1$, $W_{h_i}$ and $W_{\phi }$ have the same energy, as expected in the magnetohydrodynamic (MHD) limit (Howes et al. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2006).

Figure 2. Spectra in the perpendicular direction for the contributions to the free energy, normalised to the global level of free energy in the system.

From the GK equations (2.4), multiplying by $T_s h_s/F_s$ we obtain the balance equation for the $h_s^{2}$ variance,

(2.20)\begin{equation} \frac{T_s}{2F_s}\left[\frac{\partial h^{2}_s}{\partial t}+ \frac{c}{B_0}\{ \langle \chi \rangle_{\boldsymbol{R}_s},h^{2}_s\} + v_{{\parallel}}\frac{\partial h^{2}_s}{\partial z}\right]={q_s h_s}\frac{\partial\langle \chi \rangle_{\boldsymbol{R}_s}}{\partial t}+\frac{T_s}{F_s}h_s\left(\frac{{\partial}h_s}{\partial t}\right)_c. \end{equation}

Integrating over the velocity space, position and summing over all species we can show that we recover the evolution of the free energy (see Appendix B),

(2.21)\begin{equation} \frac{\textrm{d}W}{\textrm{d}t}=\frac{\textrm{d}}{\textrm{d}t}\int\,\textrm{d}^{3}r\left[ \sum_s\left(\int \,\textrm{d}^{3}v \frac{T_s \langle h_s^{2} \rangle_{\boldsymbol{r}}}{2F_s} - \frac{q^{2}_s\phi^{2}n_s}{2T_s}\right)+\frac{B^{2}}{8{\rm \pi}}\right]. \end{equation}

We are interested in the nonlinear contribution to the evolution of free energy for a scale, knowing that for a finite-scale system, globally, nonlinear interactions conserve $h_s^{2}$ (see Appendix A),

(2.22)\begin{align} \left.\frac{\textrm{d}W}{\textrm{d}t}\right|_{\textrm{NL}}=\int \,\textrm{d}^{3}R_s \sum_s\int \,\textrm{d}^{3}v \frac{T_s}{2F_s}\frac{c}{B_0}\{ \langle \chi \rangle_{\boldsymbol{R}_s},h^{2}_s\}=\int \,\textrm{d}^{3}R_s \sum_s\int\,\textrm{d}^{3}v \frac{T_s}{2F_s} \boldsymbol{u}_s\boldsymbol{\cdot} \boldsymbol{\nabla} h_s^{2}=0 . \end{align}

Next, we look at the GK equations and at the nonlinear contribution to the free energy balance for a system coarse-grained in the perpendicular direction in gyrocentre space.

3.1. Definition of coarse graining

The coarse graining of the Vlasov–Maxwell kinetic system was done by Eyink (Reference Eyink2018) using isotropic kernels assumed to be smooth (e.g. infinitely differentiable) and rapidly decaying (e.g.compact) phase space functions. For magnetised plasma turbulence captured by GK theory, the parallel and perpendicular scales are too disjointed in size to justify the use of isotopic filtering kernels. We concentrate here on the perpendicular scales. Accounting that the GK dynamics of interest occur in the gyrocentre space, we define the perpendicular coarse-grid filtering for a $\boldsymbol {R}_{\perp s}$ function as

(3.1)\begin{gather} \bar{a}(\bar{\boldsymbol{R}}_{{\perp} s}) = \int\,\textrm{d} \boldsymbol{R}'_{{\perp} s} G_{\ell}(\boldsymbol{R}'_{{\perp} s}) a(\bar{\boldsymbol{R}}_{{\perp} s}+\boldsymbol{R}'_{{\perp} s}) \end{gather}

(3.2)\begin{gather} = \int\,\textrm{d} \boldsymbol{R}'_{{\perp} s} G_{\ell}( \boldsymbol{R}_{{\perp} s}'-\bar{\boldsymbol{R}}_{{\perp} s}) a(\boldsymbol{R}'_{{\perp} s})=[G_{\ell}\star a](\bar{\boldsymbol{R}}_{{\perp} s}) . \end{gather}

The $\star$ symbol denotes the convolution operation. The filtering functions are considered as $G_{\ell }(\boldsymbol {R}_{\perp s})=\ell _c^{-3}G(\boldsymbol {R}_{\perp s}/\ell _c)$ with the $G$ kernels having a series of desirable properties,

(3.3)\begin{gather} G(\boldsymbol{R}_{{\perp} s})\ge 0 \quad (\text{non-negative}), \end{gather}

(3.4)\begin{gather}\int\,\textrm{d} \boldsymbol{R}'_{{\perp} s} G(\boldsymbol{R}_{{\perp} s}) =1 \quad (\text{normalised}), \end{gather}

(3.5)\begin{gather}\int\,\textrm{d} \boldsymbol{R}'_{{\perp} s}\boldsymbol{R}_{{\perp} s} G(\boldsymbol{R}_{{\perp} s}) ={\boldsymbol{0}}\quad (\text{centred}), \end{gather}

(3.6)\begin{gather}\int\,\textrm{d} \boldsymbol{R}'_{{\perp} s} |\boldsymbol{R}_{{\perp} s}|^{2}G(\boldsymbol{R}_{{\perp} s}) =1\quad (\text{unit variance}) . \end{gather}

A Gaussian kernel,

(3.7)\begin{equation} G(\boldsymbol{R}_{{\perp} s}/\ell_c)=\frac{1}{{\rm \pi}\ell_c^{2}}\exp\left({-\frac{\boldsymbol{R}_{{\perp} s}^{2}}{\ell_c^{2}}}\right) , \end{equation}

represents a good selection for the filtering function, as it obeys the properties (3.3)–(3.6). This has the advantage of a simple wave space representation, $\hat {G}(\boldsymbol {k}_{\perp }/k_c)\sim \textrm {e}^{-(\boldsymbol {k}_{\perp }/k_c)^{2}}$, which reduces the filtering convolution for the wavenumber cut-off $k_c=2{\rm \pi} /\ell _c$ to a simple multiplicative operation. However, in our work we also consider a sharp $k$-filter in wave space (i.e. a Dirichlet kernel in real space). We also consider a general (hyper-Gaussian) kernel,

(3.8)\begin{equation} \hat{G}^{\alpha}(k_{{\perp}}/k_c) \sim \textrm{e}^{-(k_{{\perp}}/k_c)^{\alpha}}, \end{equation}

that is isotropic in the perpendicular direction, knowing that for $\alpha =2$ we recover the Gaussian kernel and for large $\alpha$ we tend towards the sharp $k$-filter with respect to $k_{\perp }$. Figure 3 showcases this for the $W_{h_e}$ spectra, i.e. we filter the $h_e$ before computing the spectra.

Figure 3. Showcasing for $W_{h_e}$, the wavenumber support for different filtering kernels $\hat {G}^{\alpha }$ considered for the same three cut-offs $k_c \rho _i=\{1,4,16\}$, with $\alpha =\{2,4,8\}$ and $\alpha =\infty$ representing the sharp $k$-filter. Not shown, for $\alpha =64$ the hyper-Gaussian filter would overlap near-identically the sharp filter.

3.2. Coarse-grained GK equations

We start from the GK equations given in the advective velocity form and apply the coarse-graining operation ($G_{\ell }\star$) term by term. The overbar notation is moved only on the quantities that undergo coarse graining to obtain,

(3.9)\begin{equation} \frac{\partial \bar{h}_s}{\partial t}+ \boldsymbol{\nabla} \boldsymbol{\cdot} (\overline{{\boldsymbol{u}}_s h_s}) + v_{{\parallel}}\frac{\partial \bar{h}_s}{\partial z}=\frac{q_s F_{s}}{T_{s}}\frac{\partial\bar{\langle \chi \rangle}_{\boldsymbol{R}_s}}{\partial t}+\left(\frac{{\partial}\bar{h}_s}{\partial t}\right)_c. \end{equation}

The field equations are linear in $h_s$ and thus do not pose any complications under the coarse-graining operation. Simply replacing $h_s$ by $\bar {h}_s$ in (2.6)–(2.8) yields the coarse-grained field equations.

Natural for a nonlinear system, the $\overline {{\boldsymbol {u}}_s h_s}$ term, coarse-grained on a grid of resolution $\ell _c$, contains contributions from SGS. In the nonlinear term, to separate the purely coarse-grained contributions from any SGS contributions, we make use of the cumulant

(3.10)\begin{equation} \bar{\tau}_s=\overline{{\boldsymbol{u}}_s {h_s}} -\bar{\boldsymbol{u}}_s \bar{h}_s. \end{equation}

Now $\bar {\tau }_s$ is the term that contains all the SGS contributions to the nonlinear dynamics. An important property of $\bar {\tau }_s$ is that it is Galilean invariant by definition. Indeed, for $\boldsymbol {u}'_s = \boldsymbol {u}_s + {\boldsymbol {U}}$, with $\bar {\boldsymbol {U}}={\boldsymbol {U}}$, we have $\bar {\tau }'_s=\overline {({\boldsymbol {u}}_s +{\boldsymbol {U}}) {h_s}} -(\bar {\boldsymbol {u}}_s + \bar {\boldsymbol {U}}) \bar {h}_s =\overline {{\boldsymbol {u}}_s {h_s}} -\bar {\boldsymbol {u}}_s \bar {h}_s + {\boldsymbol {U}} \bar {h}_s-{\boldsymbol {U}} \bar {h}_s=\bar {\tau }_s$. In a similar way, $\bar {\tau }_s$ is also invariant to a $h'_s = h_s +H$ transformation, for $\bar {H}={H}$.

The nonlinear term simply becomes

(3.11)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (\overline{{\boldsymbol{u}}_s h_s}) =\boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s + \bar{\tau}_s) , \end{equation}

where we now separate the coarse-grid-scale $\bar {\boldsymbol {u}}_s$ and $\bar {h}_s$ from the SGS $\bar {\tau }_s$ contributions.

Last, considering the definition (2.11) for $\boldsymbol {u}_s$, we obtain the equivalent $\bar {\tau }_s$ formula,

(3.12)\begin{equation} \bar{\tau}_s={-}\frac{c}{B_0}{\boldsymbol{\hat{z}}} \times [\overline{ \boldsymbol{\nabla} \langle \chi \rangle_{\boldsymbol{R}_s} h_s} - \boldsymbol{\nabla} \bar{\langle \chi \rangle}_{\boldsymbol{R}_s} \bar{h}_s ], \end{equation}

which gives the SGS contribution to the nonlinear term expressed in term of the Poisson bracket structure as

(3.13)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \bar{\tau}_s=\frac{c}{B_0} [ \{\overline{\langle \chi \rangle_{\boldsymbol{R}_s} ,h_s }\} - \{ \bar{\langle \chi \rangle}_{\boldsymbol{R}_s}, \bar{h}_s\} ]. \end{equation}

3.3. The coarse-grained redistribution of free energy

We look at the evolution of the coarse-grained free energy as result of the nonlinear interactions. This is obtained from (3.11) by multiplying with $T_s\bar {h}_s/F_s$, integrating over the velocity and position space and summing over the species,

(3.14)\begin{equation} \left.\frac{\textrm{d}\bar{W}}{\textrm{d}t}\right|_{\textrm{NL}}= \sum_s \int\,\textrm{d}^{3}R_s\int \,\textrm{d}^{3}v \frac{T_s}{F_s}[\boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s + \bar{\tau}_s)]\bar{h}_s = \sum_s \bar{\varPi}_s (\ell_c) , \end{equation}

where $\bar {\varPi }_s (\ell _c)$ represents the SGS net flux of free energy through the coarse-grained scales $\ell _c$ for the species $s$,

(3.15)\begin{equation} \bar{\varPi}_s (\ell_c)= \int\,\textrm{d}^{3}R_s\int \,\textrm{d}^{3}v \frac{T_s}{F_s}[-\bar{\tau}_s \boldsymbol{\cdot}\boldsymbol{\nabla}\bar{h}_s ]. \end{equation}

As the free energy is a nonlinear invariant, for a finite-scale system (as considered numerically in the current paper), we find the SGS net flux through an infinitely small coarse-grain scale to be equal to zero,

(3.16)\begin{equation} \left.\lim_{\ell_c\rightarrow 0}\frac{\textrm{d}\bar{W}}{\textrm{d}t}\right|_{\textrm{NL}}=\lim_{\ell_c\rightarrow 0} \sum_s \bar{\varPi}_s (\ell_c)=0 . \end{equation}

Furthermore, the contributions to the free energy from each plasma species are independently invariant under the action of the nonlinear terms, i.e. $\lim _{\ell _c\rightarrow 0} \bar {\varPi }_s (\ell _c)=0$. Naturally, for an infinite-scale system given by the $Do \rightarrow \infty$ limit, where $Do\sim (k_{\perp }^{\nu _i} \rho _i)^{5/3}$ is the Dorland number (for the definition convention see Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009) defined here on the ion dissipation scale (i.e. the scale at which the finite collisional dissipation peaks in amplitude), taking the $\ell _c \rightarrow 0$ limit will give a constant flux value once parallel mixing can be neglected. In the current paper $Do\approx 228$ and strong parallel mixing affects the scaling of the electron flux.

With this knowledge, we consider the $\boldsymbol {R}_s$-density of free energy for each plasma species, i.e. $W_s(\boldsymbol {R}_s,t)$, see Appendix B, and look at its coarse-grained variation due to the action of nonlinear interactions,

(3.17)\begin{align} \left.\frac{\partial \bar{W}_s(\boldsymbol{R}_s,t)}{\partial t}\right|_{\textrm{NL}}&=\int \,\textrm{d}^{3}v [\boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s + \bar{\tau}_s)]\bar{h}_s \frac{T_s }{F_s}\nonumber\\ &=\int \,\textrm{d}^{3}v \left[ \boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}} \bar{h}_s^{2}/2+ \bar{\tau} \bar{h}_s)\frac{T_s }{F_s}\right] + \int \,\textrm{d}^{3}{v} \left[ -\bar{\tau}_s\boldsymbol{\cdot}\boldsymbol{\nabla}\bar{h}_s\frac{T_s }{F_s} \right]. \end{align}

The first term on the right-hand side corresponds to a transport in position space of free energy, whereas the second corresponds to the density of the SGS scale flux. We introduce the following definitions,

(3.18)\begin{gather} \bar{\varUpsilon}_s(\boldsymbol{R}_s,t)=\int\,\textrm{d}^{3}v \frac{T_s }{F_s}[ \boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s^{2}/2+ \bar{\tau}_s \bar{h}_s)] , \end{gather}

(3.19)\begin{gather}\bar{\varPi}_s(\boldsymbol{R}_s,t)= \int \,\textrm{d}^{3}v \frac{T_s }{F_s} [ -\bar{\tau}_s \boldsymbol{\cdot}\boldsymbol{\nabla}\bar{h}_s ] . \end{gather}

As these definitions are the two integrals from (3.17), we see that we could move the $\boldsymbol {\nabla } \boldsymbol {\cdot } ( \bar {\tau }_s \bar {h}_s)$ term from (3.18) to (3.19) in an attempt to contain the SGS contributions into a single $\bar {h}_s \boldsymbol {\nabla } \boldsymbol {\cdot } \bar {\tau }_s$ quantity. In fact, the SGS net flux (3.15) is defined up to a divergence term, and the integrant could be written simply as $\bar {h}_s \boldsymbol {\nabla } \boldsymbol {\cdot } \bar {\tau }_s$. However, although the SGS net flux is not changed by adding or subtracting a divergence term, the resulting spatial density would be different. It is hard to see a good reason to add an extra contribution to the flux density $\bar {\varPi }_s(\boldsymbol {R}_s,t)$ that does not contribute at all to the net flux across a scaleFootnote ² . As such, we decide on the definitions given by (3.18)–(3.19) for our current work.

In terms of the Poisson bracket structure, the integrant of (3.19) becomes

(3.20)\begin{equation} \bar{\tau}_s \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{h}_s =\frac{c}{B_0} \left[ \overline{{h}_s \left\{\langle \chi \rangle_{\boldsymbol{R}_s}\right.} ,\left.\overline{h}_s\right\} - \frac{1}{2}\{ \bar{\langle \chi \rangle}_{\boldsymbol{R}_s}, \bar{h}^{2}_s\} \right]. \end{equation}

For clarity,

(3.21)\begin{equation} \overline{{h}_s \left\{\langle \chi \rangle\right.}_{\boldsymbol{R}_s} ,\left.\overline{ h}_s\right\} =\overline{{h}_s\frac{\partial \langle \chi \rangle_{\boldsymbol{R}_s}}{\partial x}}\frac{\partial \overline{ h}_s}{\partial y}-\overline{ {h}_s\frac{\partial \langle \chi \rangle_{\boldsymbol{R}_s}}{\partial y}}\frac{\partial \overline{ h}_s}{\partial x} \end{equation}

and it shows why the Poisson bracket notation becomes cumbersome when dealing with coarse graining. Only terms of the form $\overline { \bullet \{\bullet } ,{ \bullet }\}$, that coarse grain across the Poisson bracket structure, give SGS contributions. From the properties of the Poisson bracket (see Appendix A) we know that the second term integrates spatially to zero ($x,y$ integration suffices). However, this second term is important to ensure the gauge invariance of the SGS flux density. Simple algebra shows that the transformation $h'_s=h_s+H$ and $\langle \chi ' \rangle _{\boldsymbol {R}_s}=\langle \chi \rangle _{\boldsymbol {R}_s}+a$, with $\bar {H} =H$ and $\bar {a}=a$ leaves (3.20) invariant. We also ask for the SGS flux density to be Galilean invariant, meaning that a change in the system of reference cannot change the intensity of turbulence, and see that the definition (3.19) fulfils this requirement. The link with Galilean invariance mentioned for $\bar {\tau }_s$ is given by ${\bar {\boldsymbol {U}}}={\boldsymbol {U}}=-({c}/{B_0}){[}\boldsymbol {\nabla } a \times {\boldsymbol {\hat {z}}}{]}$. In the same spirit, the $h'_s=h_s+H$ invariance shows that by adding or subtracting background density values to $h_s$ (during the $\delta f$ splitting, for example), we cannot change the intensity of turbulence. This also highlights why the $\bar {h}_s \boldsymbol {\nabla } \boldsymbol {\cdot } \bar {\tau }_s$ quantity does not make for a good SGS flux density, whereas (3.19) does.

We normalise the nonlinear results with respect to

(3.22)\begin{equation} \varepsilon_{\textrm{NL}}=\sum_s\left[\int\,\textrm{d}^{3}R_s(\varUpsilon_s(\boldsymbol{R}_s,t)+\varPi_s(\boldsymbol{R}_s,t))^{2}\right]^{1/2}. \end{equation}

As $\bar {\varUpsilon }_s(\boldsymbol {R}_s,t)$ is defined as the divergence of a vector field, we clearly see that it integrates to zero for periodic or appropriate asymptotic boundary conditions (i.e.$\int \,\textrm {d}^{3}R_s \bar {\varUpsilon }_s(\boldsymbol {R}_s,t) = 0$). Here $\bar {\varUpsilon }_s(\boldsymbol {R}_s,t)$ does not contribute to the redistribution of free energy across the cut-off scale. Its role is to transport free energy in position space. The nonlinear transport of free energy can be seen as being due to the coarse-grained advective velocity and due to the SGS interactions, $\bar {\varUpsilon }_s= \bar {\varUpsilon }_{{\boldsymbol {u}},s}+\bar {\varUpsilon }_{{\tau },s}$, with

(3.23)\begin{gather} \bar{\varUpsilon}_{{\boldsymbol{u}},s}=\int\,\textrm{d}^{3}v \frac{T_s }{F_s}[ \boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s^{2}/2)] , \end{gather}

(3.24)\begin{gather}\bar{\varUpsilon}_{{\tau},s} =\int\,\textrm{d}^{3}v \frac{T_s }{F_s}[ \boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\tau}_s \bar{h}_s)]. \end{gather}

The density of the SGS flux is much more interesting to us. Performing the spatial integration, we recover the $\bar {\varPi }_s (\ell _c)$ flux,

(3.25)\begin{equation} \int\,\textrm{d}^{3}R_s \bar{\varPi}_s(\boldsymbol{R}_s,t) = \bar{\varPi}_s(\ell_c), \end{equation}

which we plot in figure 4. As noted, although the integral in (3.14) recovers $\bar {\varPi }_s(\ell _c)$, it does not provide for a good definition for the SGS flux density.

Figure 4. The ion and electron SGS net flux $\bar {\varPi }_s(\ell _c)$ normalised to $\varepsilon _{\textrm {NL}}$, for $\ell _c=2{\rm \pi} /k_c$. The vertical dashed lines represent the same $k_c$ cut-off values as considered by Teaca et al. (Reference Teaca, Jenko and Told2017). We obtain the same behaviour as reported previously by Teaca et al. (Reference Teaca, Jenko and Told2017), where the flux was computed via a triple-scale decomposition. Although the electron flux exhibits a $-2/3$ scaling, the ion flux tends to be constant across the ‘inertial range’ scales. The small bottleneck around $k_{\perp } \rho _i\approx 26$ is due to the relatively abrupt transition towards the range of scales dominated by collisional dissipation.

Scaling laws predicted via functional analysis, like in the work of Eyink (Reference Eyink2018) for the Vlasov–Maxwell system, are computed for absolute values (i.e. $L^{p}$ norms). This prohibits cancellation effects from occurring when integrating any sign indefinite quantity. To make a comparison, we define the maximal (upper bound) values for the spatial transfer and SGS flux. We do so by taking the absolute value before integrating the respective quantities in velocity space:

(3.26)\begin{gather} \bar{\varUpsilon}^{\rm MAX}_s(\boldsymbol{R}_s,t)=\int\,\textrm{d}^{3}v \frac{T_s }{F_s}| \boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s^{2}/2+ \bar{\tau}_s \bar{h}_s)| , \end{gather}

(3.27)\begin{gather}\bar{\varPi}^{\rm MAX}_s(\boldsymbol{R}_s,t)= \int\,\textrm{d}^{3}v \frac{T_s }{F_s} | \bar{\tau}_s \boldsymbol{\cdot}\boldsymbol{\nabla}\bar{h}_s | . \end{gather}

Next, we present a numerical analysis of the SGS flux density and spatial transport of free energy, concentrating on one aspect at a time.

4.1. The free energy transport in position space

We plot the space density of the nonlinear transport of free energy in figure 5 for the ions and in figure 6 for the electrons. In addition, in each figure, we plot the upper bound transport density ($\bar \varUpsilon ^{\rm MAX}_s$) for the two species. Varying the cut-off value in dyadic increments (i.e. $k_c\rho _i=\{1,2,4,8\}$) allows us to observe the change in transport as smaller and smaller structures are accounted. For $\bar {\varUpsilon }_s$, the cut-off scale indicates how a structure of that size perceives the spatial transport of free energy. As the cut-off scales are taken to be smaller and smaller, we see more fine-structure being added to the transport behaviour. In particular, for the electrons, this is seen best from the plots of their upper bound transport ($\bar {\varUpsilon }^{\rm MAX}$). For the transport, while more fine structures are added for small scales, the peaks tend not to change. This is natural, as the advection of large scales by the small scales is negligible in most turbulent systems.

Figure 5. (a–d) Transport of free energy in position space for the ions, $\bar {\varUpsilon }_i(\boldsymbol {R}_s,t)$. The plots depict the typical perpendicular structures at an arbitrary $z$ slice and are normalised to their maximal in-plane value. From left to right, the $k$-filtering cut-offs are $k_c\rho _i=\{1,2,4,8\}$. (e–h) The same plots for $\bar \varUpsilon ^{\rm MAX}_i(\boldsymbol {R}_s,t)$.

Figure 6. (a–d) Transport of free energy in position space for the electrons, $\bar {\varUpsilon }_e(\boldsymbol {R}_s,t)$. The plots depict the typical perpendicular structures at an arbitrary $z$ slice and are normalised to their maximal in-plane value. From left to right, the $k$-filtering cut-offs are $k_c\rho _i=\{1,2,4,8\}$. (e–h) The same plots for $\bar \varUpsilon ^{\rm MAX}_e(\boldsymbol {R}_s,t)$.

Since globally the spatial transport integrates to zero for any coarse-grained cut-off, we look instead in figure 7 at the global variation with scale of $\bar {\varUpsilon }^{\rm MAX}_s$. Defined similarly to (3.26), we also plot in the same figure the upper bound values for the individual contributions

(4.1a,b)\begin{equation} \bar{\varUpsilon}^{\rm MAX}_{{\boldsymbol{u}},s}=\int\,\textrm{d}^{3}v \frac{T_s }{F_s}| \boldsymbol{\nabla} \boldsymbol{\cdot} ( \bar{\boldsymbol{u}}_s \bar{h}_s^{2}/2)|\quad \textrm{and} \quad \bar{\varUpsilon}^{\rm MAX}_{{\tau},s}=\int\,\textrm{d}^{3}v \frac{T_s }{F_s}| \boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\tau}_s \bar{h}_s)|. \end{equation}

As expected, doing so allows us to clearly see that the total transport is mainly due to the advective velocity and not due to the SGS terms. For ions, small-scale contributions add up fast, which we believe is due to the advective velocity and its fine perpendicular velocity structure induced by the gyroaverage, structure that cannot cancel out when taken in absolute value. In fact, past the initial large scales, the density plots for the upper bound transport are indistinguishable from the advective velocity contribution (not shown here). For reference, we plot in figure 8 a $z$-slice in the density of $\bar {\varUpsilon }^{\rm MAX}_{{\tau },s}$ for the $k_c\rho _i=2$ cut-off. Not surprisingly, the $\bar {\varUpsilon }^{\rm MAX}_{{\tau },s}$ structures are closer in shape but not location to those observed for the energy flux density, as we show next.

Figure 7. The global variation with the cut-off scale of the upper bound transport of free energy ($\bar {\varUpsilon }^{\rm MAX}_s$) for ions (a) and electrons (b). The individual contributions from $\bar {\boldsymbol {u}}_s$ and $\bar {\tau }_s$ are shown as well. All the plots are normalised to $\varepsilon _{\textrm {NL}}$.

Figure 8. A $z$ slice in the density of $\bar {\varUpsilon }^{\rm MAX}_{{\tau },i}(\boldsymbol {R}_s,t)$ (a) and $\bar {\varUpsilon }^{\rm MAX}_{{\tau },e}(\boldsymbol {R}_s,t)$ (b) for $k_c\rho _i=2$. The plots are normalised to their maximal in-plane value. Although the scale defined structures are visible as for the flux (see § 4.2), their overall distribution follow the $\bar {\varUpsilon }^{\rm MAX}_{s}$ plots.

4.2. The free energy flux density

We plot the density of the SGS flux of free energy for the ions (figure 9) and for the electrons (figure 10), respectively. The upper bound (maximal) value of the flux density for each species are presented as well. Compared with the transport density, the SGS flux density shows that as the cut-off scales become smaller, the small-scale information replaces the larger-scale information. We do not observe more fine-scale structures being added on top of a larger one, but small scales replacing the larger one. This is one of the best ways to perceive the flux of free energy across a scale (we refine further this argument to account for the filter type in § 4.3).

Figure 9. (a–d) The SGS flux density of free energy for the ions, $\bar {\varPi }_i(\boldsymbol {R}_s,t)$. The plots depict the typical perpendicular structures at an arbitrary $z$ point. From left to right, the $k$-filtering cut-offs are $k_c\rho _i=\{1,2,4,8\}$. (e–h) The same plots for the absolute SGS flux density $\bar {\varPi }^{\rm MAX}_i(\boldsymbol {R}_s,t)$.

Figure 10. (a–d) The SGS flux density of free energy for the electrons, $\bar {\varPi }_e(\boldsymbol {R}_s,t)$. The plots depict the typical perpendicular structures at an arbitrary $z$ point. From left to right, the $k$-filtering cut-offs are $k_c\rho _i=\{1,2,4,8\}$. (e–h) The same plots for the absolute SGS flux density $\bar {\varPi }^{\rm MAX}_e(\boldsymbol {R}_s,t)$.

For the ions, the SGS flux density tends to homogenise for smaller and smaller structures. The electrons show an opposite behaviour, with structures of higher intensity than the background occupying a smaller and smaller volume. These behaviours are clearly seen in figure 11, where we plot the normalised histogram of the SGS flux density values. We see the histogram tails for the electrons becoming more pronounced as the cut-off scales become smaller, whereas the ions’ values tend towards a Gaussian distribution at small scales. This is inline with the intermittency measurements performed on the distribution function in Teaca et al. (Reference Teaca, Navarro, Told, Görler, Plunk, Hatch and Jenko2019).

Figure 11. The normalised histogram of the SGS flux density $\bar {\varPi }_s(\boldsymbol {R}_s,t)$ for (a) ions and (b) electrons. Observe the slight asymmetry in favour of the positive values. For visual reference, the dashed lines depict the corresponding Gaussian distributions.

One of the advantages of measuring SGS flux density is the ability to separate positive $\bar {\varPi }^{(+)}_s(\ell _c)$ and negative $\bar {\varPi }^{(-)}_s(\ell _c)$ valued contributions to the net flux,

(4.2)\begin{equation} \bar{\varPi}^{(+)}_s (\ell_c)= \int\,\textrm{d}^{3}R_s\bar{\varPi}_s(\boldsymbol{R}_s,t), \quad \text{for}\ \bar\varPi_s(\boldsymbol{R}_s,t)>0 , \end{equation}

(4.3)\begin{equation}\bar{\varPi}^{(-)}_s (\ell_c)={-}\int\,\textrm{d}^{3}R_s\bar{\varPi}_s(\boldsymbol{R}_s,t) , \quad \text{for}\ \bar\varPi_s(\boldsymbol{R}_s,t)<0 . \end{equation}

The positive value indicates a transfer towards the small scales, whereas a negative value indicates a backscatter from small scales towards large ones. From figure 11, we clearly see that the positive branch dominates. We plot in figure 12 the positive ($\bar {\varPi }^{(+)}_s$) and negative ($\bar {\varPi }^{(-)}_s$) contributions to the net flux. The difference of the positive and negative contributions give the net flux, i.e. $\bar {\varPi }_s=\bar {\varPi }^{(+)}_s-\bar {\varPi }^{(-)}_s$. Across the entire range of scales, we see how the net flux for the electrons is the result of density cancellations of an order of magnitude higher. For the ions, a drastic cancellation is only observed up to approximately $k_{\perp }\rho _i\sim 2$. We can say that more energy is moved up and down the energy cascade in scale space than the secular-like net flux that is ultimately dissipated into heat. This is important as this diffusion in scale space has an impact on the self-organisation of turbulence. The fact that the net flux through a given scale is smaller in value than typical values of the flux density, shows the benefit of using upper bound calculations in determining the intensity of nonlinear dynamics.

Figure 12. The net SGS flux $\bar {\varPi }_s(\ell _c)$ computed as the difference of positive $\bar {\varPi }^{(+)}_s(\ell _c)$ (direct cascade) and negative $\bar {\varPi }^{(-)}_s(\ell _c)$ (backscatter) contributions for (a) ions and (b) electrons. The positive contributions dominate. For both species, the maximal flux $\bar {\varPi }^{\rm MAX}_s(\ell _c)$ given by (3.27) is plotted as well. All are normalised to $\varepsilon _{\textrm {NL}}$. The net SGS flux and the $\ell _c=2{\rm \pi} /k_c$ cut-off values depicted by vertical short dashed lines are the same as in figure 4.

4.3. Filtering kernel dependence

We want to understand whether the filtering kernel affects our results. In theory, the results should be insensitive to the type of filters used, but, in practice, especially when dealing with finite resolution numerical effects, they matter. We also state that we are less concerned with the representation of the electromagnetic fields and $h_s$ as a result of the filter (we found no visual difference; not shown), and are more concerned with the change of the SGS flux and its density.

Although a sharp filter can be seen as a scale separation, a Gaussian filter is best seen as separating compact structures in real space. With our choice of definition (3.8), we can transition from the Gaussian filter to the sharp one by increasing the value of $\alpha$. From figure 13 we see that the more compact structures of an approximate scale give way to more spread out structures of well-defined scale size. The flux is shown in figure 14, where no change in the scaling is observed once we are past the smallest of wavenumbers. However, from figure 15 we see that the distribution of flux density values has more pronounced tails for a Gaussian filter.

Figure 13. The election's SGS flux density $\bar {\varPi }_e(\boldsymbol {R}_s,t)$ for (a–d) $k_c \rho _i =2$ and (e–h) $k_c \rho _i =4$. Left to right we use $\hat {G}^{\alpha }$ with $\alpha =\{2,4,8,\infty \}$. The tendency to spatially delocalise the structures in position space, while at the same time increase the consistency of their scale size representation, can be observed for larger $\alpha$.

Figure 14. SGS net flux $\bar {\varPi }_s$ for (a) ions and (b) electrons normalised to $\varepsilon _{\textrm {NL}}$, for $\hat {G}^{\alpha }$ with $\alpha =\{2,4,8,\infty \}$.

Figure 15. The normalised histogram of the SGS flux density $\bar {\varPi }_s(\boldsymbol {R}_s,t)$ at $k_c \rho _i=2$ using a Gaussian filter ($\hat {G}^{2}$) and a sharp one ($\hat {G}^{\infty }$) for (a) ions and (b) electrons. The inset pictures show the same for $k_c \rho _i=16$. We see that the Gaussian filter has increased tail contributions and that the positive branch (direct cascade) dominates. For visual reference, the dashed lines depict the corresponding Gaussian distributions.