Gyrokinetic turbulence: between idealized estimates and a detailed analysis of nonlinear energy transfers

Using large resolution numerical simulations of GK turbulence, spanning an interval ranging from the end of the fluid scales to the electron gyroradius, we study the energy transfers in the perpendicular direction for a proton-electron plasma in a slab magnetic geometry. In addition, to aid our understanding of the nonlinear cascade, we use an idealized test representation for the energy transfers between two scales, mimicking the dynamics of turbulence in an infinite inertial range. For GK turbulence, a detailed analysis of nonlinear energy transfers that account for the separation of energy exchanging scales is performed. We show that locality functions associated with the energy cascade across dyadic (i.e. multiple of two) separated scales achieve an asymptotic state, recovering clear values for the locality exponents. We relate these exponents to the energy exchange between two scales, diagnostics that are less computationally intensive than the locality functions. It is the first time asymptotic locality is shown to exist for GK turbulence and the contributions made by highly non-local interactions, previously reported in the literature, are explained as very local transfers of energy that occur between wavenumbers within the same dyadic signal. The results presented here and their implications are discussed from the perspective of previous findings reported in the literature and the idea of universality of GK turbulence.


I. INTRODUCTION
Plasma turbulence is ubiquitous, being found in astrophysical [1] and laboratory [2] settings. Although laboratory experiments [3,4] offer a controlled environment, it is often solar wind measurements [5] that allow for plasma turbulence to be probed across a wide range of scales [6]. While the dynamics at large scales can be captured by fluid approximations [7,8], the physics of weakly collisional plasma turbulence requires a kinetic description [9] for scales comparable to the proton gyroradius and smaller. Typically, the kinetic scales are associated with the dissipation range of turbulence. In this interval, which is important for the plasma heating problem [10,11], kinetic effects such as cyclotron [12,13] and Landau damping [14,15] need to be considered. However, a kinetic Alfvén wave [16,17] and entropy cascades [18,19] develop in the same kinetic scales interval, through a process similar to the one found in fluid turbulence.
As with neutral fluids, turbulence in a plasma develops due to the existence of couplings between system scales, albeit with a series of complications resulting from the interactions with the electromagnetic field. At kinetic scales these complications are increased even further, as the underlying interactions manifest themselves in a position-velocity phase space. For a kinetic plasma, the energy redistribution due to nonlinear interactions and the particle-wave resonance (e.g. cyclotron and Landau resonance) cannot be studied separately, as these phenomena influence each other in a turbulence setting [14,20]. For example, in magnetised plasmas an anisotropy develops in the position and velocity space, which influences the balance between the linear phase mixing [21][22][23] (that includes linear Landau damping) and nonlinear phase mixing [18,[24][25][26][27] that occurs in the perpendicular direction and is caused by the same nonlinear term responsible for the generation of the turbulence cascade. Perpendicular plasma structures, generated through nonlinear interactions, can be damped in the parallel direction through Landau damping, if a balance between the damping rates and the nonlinear turnover time can be achieved. While this is a problem that is receiving more interest within the community [23,[28][29][30][31] and is far from being solved, analysing the effective perpendicular cascade in position space can offer an important insight into the properties of kinetic turbulence.
Although in kinetic plasma turbulence the electromagnetic fluctuations take part in the particlewave interactions and mediate the nonlinear couplings, the balance between linear and nonlinear phase mixing can be seen as a species dependent process. It is easy to see that particle-wave resonance conditions depend on the characteristics of the particles. However, the fact that the nonlinear interactions conserve free energy for each plasma species independently, makes the energy cascade a species dependent process as well. Since the properties of the nonlinear energy redistribution in kinetic turbulence mirror the ones found in classical fluid turbulence, it is best to investigate the energy cascade for each plasma species independently. This allows us to maximise the applicability of the lessons learned from classical fluid turbulence to our current study of gyrokinetic (GK) turbulence.

A. A small overview of past studies on energy transfers and locality in turbulence
The unsolved problem of turbulence has been posed and analysed extensively in the framework of fluid dynamics [32]. In plasmas, the incompressible magnetohydrodynamic (MHD) limit represents the simplest mathematical representation for the turbulence problem. MHD turbulence, considered either for a simple electrically conducting fluid or seen as capturing large scale plasma effects, shares a lot of its properties with its electrically neutral fluid counterpart. While its Alfvénic nature, i.e. the existence of Alfvén waves that affect the nonlinear interaction time and the sweeping/straining motions of turbulent structures [8], does lead to a series of particularities not seen in neutral fluids [33][34][35][36][37][38], MHD turbulence can still be seen as the nonlinear coupling of scales which leads to an intermittent energy cascade.
In turbulence, the concept of the energy cascade represents the phenomenological interpretation of the effective energy transfer that occurs between any two scales, as a result of the nonlinear interactions [39,40]. For three-dimensional turbulence, the direct energy cascade (i.e. from large to small scales) is assumed to be local by Kolmogorov scaling estimates. However, only in the early 90's did numerical simulations allow the diagnostics that measure the transfers between two scales to be computed directly from the nonlinear terms, and the direct and local character of the cascade to be shown explicitly. These diagnostics were introduced for neutral fluid turbulence [41][42][43][44], then ported to drift-wave plasma turbulence [45,46] and latter extensively used for MHD turbulence [47][48][49][50][51][52][53][54]. Recently, they were analysed in the context of gyrokinetic plasma turbulence [55][56][57].
In addition to the characterisation of the energy cascade as local, the problem of locality of nonlinear interactions was studied directly in fluid [39,[58][59][60][61][62][63][64][65], MHD [66][67][68] and gyrokinetic [69][70][71] turbulence. While the locality of the cascade and the locality of the nonlinear interactions are related, the two concepts possess different characteristics, as we will explore in the present article for gyrokinetic turbulence. Next, we try to clarify the ideas and definitions related to locality in turbulence.

B. Clarifying the meaning of locality in turbulence
While not difficult as a generic concept, the various terms used to refer to scale locality in turbulence can become confusing, especially when a detailed analysis is attempted. We clarify what we mean by these terms, which will be used in the current article.
The locality of the energy exchanging scales refers to the separation between the energy giving and the energy receiving scales that take part in an energy transfer. This can be evaluated regardless of the information of the mediator scale (i.e. the scale of the advecting field, the third scale involved in a coupling), or by integrating first over all possible mediators. In the latter case, we recover the locality of the energy cascade problem (referred simply as the locality of the cascade).
By comparison, the locality of the nonlinear interactions accounts explicitly for the mediator scale. It relates the intensity of an energetic coupling with the maximal separation that exists between the energy receiving scale and the other two scales. For example, a strong energy transfer between two close scales that is mediated by a much larger scale will contribute to the local character of the cascade, while at the same time it will enhance the nonlocal character of the nonlinear interactions. The locality of the nonlinear interactions is measured by locality functions [39,70] and is traditionally referred to as the locality problem in turbulence (or simply as the locality problem). While in both situations we relate the intensity of energy exchanges with the scale separation, the locality of the cascade can be seen as being included in the locality problem.
As nonlinear interactions become scale invariant in the inertial range, i.e. the range where all other interactions are subdominant, the locality problem is also expected to become scale invariant. This implies that the intensity of the energy transfers decreases in the same manner with the increase in separation, regardless of the value of the energy receiving scale (i.e. the reference scale from which we measure the separation in all cases). We now say we recover asymptotic locality, since once turbulence develops an inertial range the locality problem does not change further. In the inertial range, expressing the decreases in intensity of the energy transfers as a power law of the scale separation yields an exponent, which we call asymptotic locality exponent.
The asymptotic locality exponent is a characteristic of turbulence. If by varying the parameters that define the system (e.g. plasma parameters) we obtain the same asymptotic locality exponent, then the nonlinear interactions are invariant in regard to these parameters. This invariance of the nonlinear interactions in the inertial range will ensure that any information introduced at large scales (e.g. via large scale forces or boundaries) will be destroyed (decorrelated) through the cascade process. At the end of the cascade the same information is recovered, which leads to the small scales to be universal. Traditionally this is referred to as the universality of turbulence problem, since universal small scales lead in classical fluids to universal dissipative mechanisms. However, we can allow the small scales to be directly affected by some linear mechanism, while at the same time having a universal character for the nonlinear interactions. Defining the universality of turbulence as the universality of nonlinear interactions is appropriate, as it is the latter we need to have to be able to develop unique sub-grid scale models. Regardless of the accepted definition of universality, a unique asymptotic locality exponent is seen as a necessary condition for the existence of universality in turbulence.

C. Structure of the article
In the current paper, we study the effective energy transfer in the perpendicular direction for a magnetised plasma described by a gyrokinetic (GK) formalism. We are interested in describing the energy cascade and the locality problem for GK turbulence. This analysis relates to the fundamental question of universality of turbulence and in particular the universality of plasma turbulence. Starting from the study of the energy redistribution and the scale locality problem, we show that the general nonlocal nature of GK turbulence, captured via locality functions, contains a subset of interactions that are deemed local, are scale invariant (i.e. a sign of asymptotic locality) and possess an exponent that can be recovered directly from the interaction of solely energy exchanging scales (i.e. the energy cascade).
In section II, we present the GK equations for a slab magnetic geometry and list the parameters of the nonlinear simulation employed throughout this work. In addition, we derive the free energy balance equation for a scale and define a norm for the intensity of the energy transfers in the system. In section III, we discuss succinctly the magnetic geometry effects on the scale representation, present the interaction conditions between three scales and the resulting implications on their separation and introduce the waveband decomposition for the GK system. In section IV, after building the nonlinear energy transfers and the scale flux diagnostics, listing their properties and interconnections, we proceed to present the results for the numerical simulation employed. This section presents diagnostics which can be seen up to a point as classical tools for the analysis of turbulence, i.e. tools used in the past for the analysis of the equivalent problem in fluids.
In an attempt to understand better the connection between energy transfers and the locality of interactions, we introduce in section V an idealised test problem. In section VI, employing the lessons learned from the test problem and applying a series of detailed considerations, we conduct a further analysis of the large resolution GK simulation data. We find that the locality exponents for the energy cascade exhibit an asymptotic behavior, denoting the possibility of a universal character for the energy cascade in GK turbulence. Furthermore, we show that these exponents can be obtained directly from the energy transfers between two scales rather than the computationally intensive locality functions (modified in section VI to account only for the locality of the energy cascade). These implications and their link with past results presented in the literature are discussed last, in section VII.

II. THE GYROKINETIC FRAMEWORK
The gyrokinetic formalism represents a rigorous limit [72] of kinetic theory for strongly magnetised plasmas for which gyrotropy is assumed (i.e. invariance under the gyration of particles of charge q σ and mass m σ around the magnetic guide field of intensity B 0 , for each plasma species σ). We reproduce below the main ideas behind GK for the simpler context of astrophysical plasmas [73].
At the kinetic level, the distribution functions expressed at the particle position f σ (r, v, t) represent the dynamical quantities of interest. The role of the self-consistent electromagnetic fields, obtained from velocity moments of the particles' distributions, is to mediate the linear and nonlinear interactions between structures in the distribution functions. The most common approach for kinetic turbulence is to assume small fluctuations around equilibrium background distribution functions F σ , considered usually to be Maxwellians. In the GK formalism, the dynamics are contained by the perturbed gyro-center distribution functions h σ (x, y, z, µ, v , t), where (x, y, z) are magnetic coordinates, with the z direction coinciding with the direction of the guiding magnetic field lines. Aligning the representation of the system with the direction of homogeneity induced by gyrotropy allows us to remove the gyration motion (of gyroradius ρ σ = √ T σ m σ c/eB) from the kinetic system and effectively reduce the phase space to just five-dimensions (i.e. x, y, z, µ, v ). This reduction is important for the numerical implementations [74], as it substantially saves computational resources by considering only two velocity directions. The velocity along the magnetic field line is v . The magnetic moment (µ = m σ v 2 ⊥ /2B 0 ) is an adiabatic invariant for the GK system and contains implicitly the perpendicular velocity v ⊥ information. Typically, the velocity directions are expressed in thermal velocity units (v th σ = 2T σ /m σ ) and the equilibrium (background) density n σ and temperature T σ are known.
Accounting for a Boltzmann response factor due to the process of restoring plasma electroneutrality and writing explicitly only terms up to the first order in the small parameter introduced by the GK ordering (notably low frequencies for the plasma fluctuations compared to the ion, here proton, cyclotron frequency Ω σ and small fluctuation levels), we obtain [73] a relation between f σ (r, v, t) and the perturbed gyro-center distribution functions h σ (x, y, z, µ, v , t), For such an expansion to be valid and for the removal of the particles' fast gyro-motions to be done systematically, the GK ordering must hold from a physics perspective. In strongly magnetised plasmas, the GK formalism represents a rigorous kinetic representation of the turbulence problem, being able to capture the kinetic Alfvén wave cascade [75,76], as well as the entropy cascade [18,19] and the linear phase mixing associated with Landau damping. In astrophysical studies, while it neglects cyclotron resonance, gyrokinetics can still be seen as a useful tool as it captures [77] the crucial dynamics of kinetic Alfvén wave (KAW) turbulence in three spatial dimensions [78] while offering a more manageable system to be simulated numerically. The GK equation is just the Vlasov equation rewritten for h σ (x, y, z, µ, v , t) and considering gyrotropy, for which the electromagnetic potentials are determined from electromagnetic sources obtained at the location of the particles' gyroradius (i.e. the electromagnetic sources are effectively rings of electric charge centered on the gyro-centers).

A. The gyrokinetic equations
The nonlinear gyrokinetic equations [79] were derived rigorously by Ref. [80] and presented extensively in [72]. An appropriate review of GK turbulence for newcomers is presented by Ref. [81] and in the simplifying context of astrophysical plasma the GK formalism is presented in [73]. In slab magnetic geometry, for h σ = h σ (x, y, z, v , µ) the gyro-center distribution functions, the gyrokinetic equation for a species σ has the form ,σ is the gyrokinetic (gyro-averaged) potential and φ, A and B are obtained [82,83] from their respective field equations for known sources (n σ , j ,σ , p ⊥,σ ), which in turn are obtained from the velocity moments of h σ at the location of the particles' gyroradius. As this is expressed in a simpler form in k ⊥ space, we list the sources equations for a mode k ⊥ as where J 0 is the Bessel function, I 1 is the the modified Bessel function and for k ⊥ ≡ |k ⊥ | we have [84]. The term ∂hσ ∂t coll refers to the impact made on the time evolution of h σ by a linearised Landau-Boltzmann collision operator (see the supplemental material provided by Ref. [85] for the exact form used here). While not standard, we absorb a minus sign in the definition v σ , which is minus the drift velocity, as a way to achieve a more compact notation in the next sections.

B. Nonlinear simulation data
In this study we use gyrokinetic simulations of magnetised proton-electron plasmas. The nonlinear gyrokinetic system of equations is solved using the Eulerian code GENE [86]. The data used in this work is taken from the simulation presented in Ref. [71], one of the largest GK simulation to date, and it is briefly summarised below. We mention that, while our analysis is limited to the use of this pre-existing large resolution simulation, and an even larger velocity space domain and velocity space resolution would be needed ideally for diagnosing the overall phase space mixing problem, we are confident in the results presented here, results which relate to the energy cascade in the perpendicular spatial direction and the locality of interactions problem.
The physical parameters of the simulation are chosen to be close to the solar wind conditions at 1 AU, with β i = 8πn i T i /B 2 0 = 1 and T i /T e = 1. Proton and electron species are included with their real mass ratio of m i /m e = 1836. The electron collisionality is chosen to be ν e = 0.06 ω A0 (with ν i = m e /m i ν e ), and ω A0 being the frequency of the slowest Alfvén wave in the system. The evolution of the gyro-center distribution is tracked on a grid with the resolution {N x , N y , N z , N v , N µ , N σ } = {768, 768, 96, 48, 15, 2}, where (N x , N y ) are the perpendicular, (N z ) parallel, (N v ) parallel velocity, and (N µ ) magnetic moment (µ = m σ v 2 ⊥ /2B 0 ) grid points, respectively. This covers a perpendicular dealiased wavenumber range of 0.2 ≤ k ⊥ ρ i ≤ 51.2 (or 0.0047 ≤ k ⊥ ρ e ≤ 1.19) in a domain L x = L y = 10πρ i . In the parallel direction, a L z = 2πL domain is used, where L ≫ ρ i is assumed by the construction of gyrokinetic theory. A velocity domain up to three thermal velocity units is taken in each direction. The fluctuations in the system are driven via a magnetic antenna potential (a A ant contribution is added to χ), which is prescribed solely at the largest scale and evolved in time according to a Langevin equation [87]. We mention that this antenna potential is removed fromχ σ before the nonlinear diagnostics are computed.

C. The free energy balance equation for a scale
The free energy for a species (E σ ) represents the quadratic quantity of interest for the study of the dynamics of gyrokinetic turbulence [19]. Free energy is the quantity that is injected into the system by instabilities or external drives, dissipated by collisions, while being redistributed in a conservative fashion by the action of the nonlinear terms (E σ is nonlinear conserved for each species σ independently). The free energy for a species σ is defined using the GK system as where the · · · notation stands for the average over the phase space domain, including the appropriate five-dimensional Jacobian (J 5D ) contributions, · · · = (· · · )J 5D dx dy dz dv dµ J 5D dx dy dz dv dµ .
See Ref. [88] appendix B for a full derivation of the free energy definition for GK, starting from its classical form (entropy contribution + electric energy + magnetic energy). While the link with eq. (9) is conceptually straightforward, the derivation is too tedious to be reproduced here and we consider the free energy definition for GK theory to be granted by eq. (9). Rewriting eq. (2) as we can formally obtain the global free energy balance equation for the species σ by applying the operator T 2F h σ · · · on each of the terms in eq. (11). Since the free energy is a nonlinear invariant, we have T 2F h σ (v σ · ∇h σ ) = 0. A much more useful balance equation is obtained for a hierarchy of scales, naturally provided by a Fourier decomposition, by considering The free energy balance equation for a scale is now simply obtained by applying the filtered operator σ · · · on each of the terms in eq. (11) and has the form The T σ (k) contains the contribution of all nonlinear interactions for a scale identified by k and is known as the transfer spectrum. The L σ (k) term contains the linear parallel dynamics (including linear Landau damping). As long as the filtering condition does not depend on z (e.g. dk/dz = 0, the case considered here), L σ (k) is zero for all k's due to the {z, v } integration. Otherwise, the geometry leads to a (linear) k-redistribution of energy that integrates to zero over all k-scales. Last, D σ (k) is a dissipative term that appears due to the presence of collisions.

D. Choosing a norm for the energy transfers
The transfer spectrum T σ (k) integrates to zero over all scales, a result of the conservation of free energy by the nonlinear interactions. It represents the simplest quantity related to the nonlinear energy exchanges that can be computed numerically and it can be recovered from all other nonlinear diagnostics that account for additional scale decompositions of the nonlinear term, as we will see in later sections. Taking all these facts into account, we consider that T σ (k) can serve as a basis for a useful norm that will allow us to gauge the intensity of various energy transfers. We formally define this norm ε σ to be: In general, this is a definition that requires no particular shape for a T (k) curve, except that it integrates to zero over the entire scale domain considered. For steady state turbulence that exhibits a clear inertial range (i.e. the value of T (k) goes from negative to positive with the increase in k and is zero in the inertial range), ε defined above recovers the scale-invariant value for the energy flux in the inertial range. However, the use of this definition is also appropriate for transient state turbulence, for which an inertial range is not observed and a representative energy flux value cannot be determined for the system. Having units of power, ε can act as an indicator of the intensity of the nonlinear energy exchanges present in the system (here for each species σ) and it can be used to compare turbulence between various states and simulations.

III. THE REPRESENTATION OF PERPENDICULAR SCALES
Since the gyrokinetic formalism is strongly anisotropic and assumes by construction that k ⊥ ≫ k , the main effect of the nonlinear interactions is to mix the perpendicular spatial scales. While we recall that the underlying GK dynamics occur in a five-dimensional phase space and the nonlinear phase mixing does become important at scales k ⊥ ρ σ > 1, the effective energetic interaction between three perpendicular modes (k ⊥ +p ⊥ +q ⊥ = 0) can still be measured and the resulting perpendicular scale interactions can be analysed. To simplify the notations, we use k from now on to refer to the perpendicular wave vector k ⊥ . Considering that a scale ℓ can be defined by the norm of a wave vector (e.g. ℓ ∼ 1/k), we will typically identify a scale by the norm k and, by abuse of language, refer to it as a k-scale, even though the k norm has units of inverse length.
A. The impact of the magnetic geometry on the perpendicular scales Before analysing the interactions between perpendicular scales, we stop to talk briefly about the impact made by the geometry of the magnetic guide field on the scale representation. In the (k, z) space, let us consider the magnetic geometry to be prescribed via a contra-variant metric tensor that varies only along the magnetic field, i.e. η ij = η ij (z). This is one of the simplest scenarios, that of the local approximation of the magnetic flux surfaces typically used in tokamak studies [89]. The perpendicular scale, considered as the wave norm k ≡ |k|, is found as We immediately notice that the perpendicular scales have now a z dependence. This simple fact complicates the scale decomposition typically employed in turbulence studies, which now needs to be done in (k, z) rather than simply in k. While in a slab geometry, i.e. the magnetic configuration of choice for astrophysical studies, η xx (z) = η yy (z) = 1, η xy (z) = 0 and dk/dz = 0, let us visualised in figure 1 the k(z) norm for a sheared box given by η xx (z) = 1, η xy (z) =ŝz, η yy (z) = 1 + (ŝz) 2 which recovers the slab configuration forŝ = 0.
The z dependence for a scale is not just a nuisance that complicates the scale decomposition. The nonlinear interactions occur between a triad of resonant modes (k + p + q = 0). In the slab case the interaction of any three modes leads to an interaction solely between three perpendicular scales. Once magnetic curvature effects are considered the same triad interaction now couples multiple scales together, as the same k x , k y mode contributes to multiple k(z) scales. Furthermore, in the free energy balance equation (eq. 13) we now have L(k) = 0 (without going into details, this can be seen as a variation in terms of z leading to a variation in terms of k). This shows a split between the magnetic coordinates used to describe the GK system (i.e. k x , k y , z) and the natural coordinates used for the description of scales in turbulence (i.e. k); a nontrivial problem that requires further analysis from the plasma scientific community. This is the main reason why the slab configuration, which we resume to from this point on, is preferred as a basis for understanding the basic dynamics of gyrokinetic turbulence, even though the general gyrokinetic theory can account for magnetic geometry effects. However, fully understanding the fundamental characteristics of plasma turbulence in arbitrary magnetic geometry is a highly desirable scientific proposition.
FIG. 1. The iso-surface of k(z) for z ∈ [−π, +π] and for the kx, ky domains centered on zero, respectively. The slab case corresponds toŝ = 0. Note that the largest k value that is fully captured in the same (kx, ky, z) domain, as considered in all the figure's panels, decreases with the increase of the shearŝ, as evident from the z mid-plane for which k(z = 0) is a circle of diminishing radius.

B. Interaction conditions for scales
As mentioned, the nonlinear energetic interactions occur due to all possible triads, i.e. three modes for which their wave vectors obey the triad condition This triad condition imposes a limit on the interaction of scales. Using the fact that |k+p+q| = 0, we have the triangle inequalities: These conditions tell us that for given k and p scales, the third scale q that enters the nonlinear interaction, respecting the triad condition given by eq. (16), needs to obey: An example is visually represented in figure 2 and more detailed pictures of allowed scale interactions can be found in Ref [90]. We see that for k ∼ p the scale q can be at most 2k and at least 0. A further consideration on scale separation is undertaken next.

C. Dyadic separation of scales
Let us first consider q ≤ k, p ≤ k and k as three scales coupled by a nonlinear interaction, where k denotes the smallest scale initially available in the system. The nonlinear interaction in question can potentially generate scales smaller than k. From eq. (19) we see that any new smaller scale is contained in the interval [k, 2k], i.e. at most a factor of 2 smaller. A phenomenological interpretation can be given: the shearing of a structure by a larger advecting flow can generate at most scales half the structure's size.
We will refer to the interval [k, 2k] as a dyadic band [62], and refer to a separation between p = k/2 and k as a dyadic separation. In music, a dyadic band represents an octave. Let us now consider a dyadic (octave) separation of scales (i.e. ...k/4, k/2, k, 2k, 4k...) and assess the implications for the interaction of two scales thus separated. The interaction of k with any p ≤ k/2 is mediated by the scales q ∈ [k/2, 3k/2], according to the triangle inequalities introduced above (section III B). This means that only two dyadic bands, one on each side of k captures the entire interaction process (that is, two dyadic separated scales can be seen as interacting rather locally). The scale q can be less than k/2 and tending towards zero (thus increasing the nonlocal character of the interaction) only when scale p is less than a dyadic separation away from k (i.e. k/2 < p < k). For a given k and k/2 < p < k, we designate their coupling as an intra-dyadic interaction. We also notice that for the separation from k of both q < k and p < k to be equal, it is required that q = p = k/2, i.e. both q and p are dyadic separated from k. In later sections, these definitions will help us in the classification of local and nonlocal interactions.
Far from being a simple thought experiment, physical scales do need to be separated sufficiently (in a log(k) space) to be distinguished from each other (e.g. the octave separation of frequencies for musical notes being a classical example in this sense). Thus, a flow structure with scales contained by a wavenumber band can be better associated with the concept of an eddy in turbulence, i.e. a structure well-localised at the same time in real space and in wave space, opposed to the case of a simple linear wave that is identified by a single wavenumber and which is completely un-localised in real space. While a dyadic band decomposition can ultimately be seen as arbitrary, the interaction between two dyadic separated scales can prove to capture better the phenomenological interaction between two scale structures in turbulence [66] and provide a simpler link with the classical phenomenological interpretation of turbulence. Next, we introduce the waveband decomposition for the scales that we will actually use in the measurement of the energy transfers.

D. Waveband representation
As in previous works for GK turbulence [56,[69][70][71], we define a series of scale intervals s n = [k n−1 , k n ], with boundary wavenumbers given as a geometric progression, for n ∈ N * , λ > 1 and k 0 = 0. These structures are called shells in previous works [47][48][49][50][51][52][53][54] (here having the geometric shape of cylindrical shells in a (k, z) space) or bands [41,58,61,66,91], to account for the fact that they represent bands of equal width in a log(k) space. The distribution functions or electromagnetic potentials are filtered in wave space, obtaining their respective band (shell) filtered contributions. For example, the waveband filtered distribution functions h [n] (k) are found in k space as while the real space contributions are simply obtained as It is important to realise that the filtered signals are well defined in real space, the total information being recovered as the superposition of all filtered contributions, e.g.
and that they are orthogonal to each other, i.e. h [n] (x, y)h [m] (x, y)dxdy = 0 for n = m. We mention that a decomposition using infinitesimally thick bands could be performed, equivalent to the recovery of the wave-norm k-scale splitting prescribed by eq. (12). However, a geometric progression is preferred for turbulence studies, since scaling laws play an important part and, we want to separate physical structures without wasting numerical resources.
For the current GK study we take a total of N = 25 wavebands, with k 1 = 0.275ρ i and λ = 2 1/3 . While a dyadic separation (i.e. λ = 2) is most useful for the analysis of energy transfers between structures with a more robust phenomenological equivalence, i.e. eddies, our choice for the λ factor allows us to perform a finer analysis of the nonlinear interactionsvia the analysis of scale fluxes.
Past studies used 2 1/4 (Ref. [68]) and 2 1/5 (Ref. [69]) as the values for λ. In the current work, which makes use of Ref. [71] large scale GK computations, the λ = 2 1/3 choice was done to reduce the number of bands required to spawn the k interval. This is a choice dictated by computational costs, as the most complex diagnostic requires the calculations of the nonlinear term N 2 times.

IV. NONLINEAR ENERGETIC INTERACTIONS IN GK TURBULENCE
In the current work, instead of using the triad transfers (i.e. the nonlinear energy transfers that occur between three modes which respect the k + p + q = 0 resonant condition) as the basis for the conceptual definition of various energy transfers [69][70][71], we employ an alternative presentation. While the two approaches are in fact equivalent, we consider the following introduction of the transfer of energy between three scales to be easier to grasp.

A. Building the triple-scale transfer
Suppressing the plasma species index and using for each species the implicit form of the advecting velocity, i.e. v = − c B0 e z ×∇χ , the nonlinear term entering the GK equation has the compact expression v · ∇h. The global variation of the species free energy due to the nonlinear term has now the form and is zero due to the conservation of free energy by the nonlinear interactions. Considering the . Performing a similar scale splitting on ∇h and v we obtain the equivalent statement for the conservation of free energy by the nonlinear interactions, At this point we make two remarks. First, while globally the energy transfers integrate (sum) to zero, the individual transfers T can have any value and will form the basis for our triple-scale transfers. And second, we see that the splitting employed is true for any decomposition of our quantities, not just for a waveband decomposition. In fact, it is up to us to provide the proper physical scale decomposition and justify what can be seen as an arbitrary choice. We remind the reader that K, P , Q are integers that identify the waveband intervals and are not themselves a wavenumber, e.g. K is the integer that identifies the band spanning the Employing a waveband scale decomposition prescribed by eq. (23), we define here the triple-scale transfer as which measures the energetic interaction between three waveband prescribed scales. As the role of the velocity v is to advect the spatial gradients of the distribution h, the scales on position Q have the role of mediating the transfers between scales P and K. We furthermore consider that the scale K receives energy if the transfer is positive, a choice consistent with the interpretation used in past studies [69][70][71].
Conceptually, considering wavebands of infinitesimal thickness in the continuous limit of the spectral space, we can obtain S(k|p|q) as the transfer between three wave-norm scales. While we will use this to provide some formal definitions that will be simpler to grasp, we emphasise that we only have access to the power law waveband decomposition, which is the most efficient choice from a computational point of view.

B. Properties of the triple-scale transfer
Assigning specific values to K, P and Q, we list the properties for the triple-scale transfer. For each mediator Q, the amount of energy received by scale K is opposite the amount of energy given by P , This can be shown from eq. (28), using derivation by parts, accounting for the periodic boundaries employed here and considering that ∇ · v [Q] ≡ 0, which results from the definition of v. From the relation (29) we trivially find that S(K|K|Q) = 0, i.e. the energy transferred from one scale to itself is zero. The conservation of energy implies that the sum of all transfers occurring between the same three scales is zero, and is shown to be true by employing the antisymmetry property given by eq. (29). All the properties listed above will be inherited by subsequent diagnostics that are constructed on the triple-scale transfers S(K|P |Q).

C. The definition of energy transfer diagnostics and the link between them
For GK turbulence, the scale-to-scale (shell-to-shell [48]) transfers have been studied before in the literature [55,56]. They represent one of the first types of nonlinear diagnostics to be adopted by the field of plasma turbulence [45] from the field of hydrodynamical (classical) turbulence [42]. The scale-to-scale transfers are defined from the triple-scale transfer or directly from a waveband decomposition of the nonlinear term, as A scale-to-scale transfer has the interpretation of the energy received by a scale K from the scale P , accounting for all possible mediations. Due to the conservation of energy, P(K|P ) = −P(P|K) and P(K|K) = 0 for each species. If P(K|P ) is determined directly from the nonlinear term, only N calculations (i.e. the number of bands) of the nonlinear term are required. By comparison, S(K|P |Q) requires N 2 computations of the same type. From the scale-to-scale transfers or directly from the triple-scale transfers, we recover the nonlinear transfer spectrum, defined here as While we can recover the transfer spectrum from more complex scale decompositions, we mention that knowing the full nonlinear term, required for the incremental integration of the GK equations, is sufficient for the computation of T .
Succinctly, the links between the nonlinear transfer spectrum, the scale-to-scale transfers, the triple-scale transfers and the global conservation of free energy by the nonlinear interactions (the sole property needed for their definition), can be summarised as We will present next the transfer spectrum and the the scale-to-scale transfers at a given instant in time for the GK simulation analysed in this paper (section II B).

D. Transfer spectrum
We start by mentioning that for each species σ we approximate the norm ε in terms of the waveband representation for the transfer spectra, This is a practical approximation and the use of the discrete version given by eq. (36) can be seen as being acceptable as long as T (k) is not highly fluctuating in a K interval. For our GK simulation (section II B) we have ε i = 9.65 × 10 1 [GENE power units] and ε e = 2.16×10 2 [GENE power units]. Looking at the ratio ε e /ε i ≈ 2.24, we see that the electron's energy transfers are the most intense in the GK system analysed. Here, from the transfer spectrum T (K) for ions and electrons (presented in figure 3) we observe that the electrons remove more energy from the large (forced) scales. In the absence of time averages, it is hard to properly distinguish properties of the transfer spectra.

E. Scale-to-scale transfers
The scale-to-scale diagnostic provides a way to visualise the energy cascade. Since the waveband boundaries are taken as a power law, the scale-to-scale transfers normalised to their maximal absolute value provide us with information regarding the direction and locality of the energy cascade. We designate a transfer to be direct if it's positive for K > P and we call it local if it occurs primary between scales with P ∼ K. In section VI, we will elaborate on the local character of the cascade. From figure 4 for the ions and figure 5 for the electrons, we do observe that the patterns for the scale-to-scale transfers correspond indeed to a direct and local energy cascade. Since P(K|P ) is systematically positive for the energy received from larger scales K > P (lower-diagonal in panels a), we can say that we observe a direct energy cascade. For both sets of figures, the pseudo-isometric visualisation (panels b) allows us to better compare the intensity of the transfers at different scales, while plotting the scale-to-scale transfers as a function of P − K (panels c) allows us to better gauge their locality and self-similarity (i.e. the curves for different K's would collapse on each other for perfect self-similar transfers; an expected property for scale-to-scale transfers in the inertial range of classical homogenous turbulence).
We again mention that it's important to differentiate between the locality of the energy cascade, one structure giving energy to a similar size structure, and the locality of interactions captured by Kraichnan's locality functions [33], where the mediators of the energetic interaction between two scales are also considered. While the two are related, as we will show in section VI and discuss in the last section, they are not directly equivalent.
For the system analysed here (β = 1), we notice that the ion's scale-to-scale transfers tend to be self-similar in a range starting with the gyroradius (kρ i ∼ 1) and ending around kρ i ∼ 20. The latter limit is due to the collisional dissipation employed, which for the ions is dominated by the field perpendicular contributions (including the k ⊥ Finite Larmor Radius -FLR effects) [85]. Considering the overall shape of T i (large scale source, small scale sink), the ions match the classical picture of turbulence to a certain degree. By comparison, electrons suffer from strong parallel mixing effects at scales kρ i > 1, which help to remove free energy through parallel collisions. For electrons dominated by Landau damping, the role of the perpendicular cascade is just to mop up the reminder of the free energy and pass it down to ever smaller scales, which are in turn also affected by Landau damping. This leads to the "damped" cascade picture observed in figure 5. As the electron dissipative route is more efficient than the ion one, there should be no surprise that in a steady state the electron free energy channel draws in more free energy from a given source (leading to ε e /ε i ≈ 2.24).
where the last two identities relate the scale flux to the scale-to-scale transfer and the transfer spectrum, respectively. We display in figure 6 the scale flux for the ions and electrons for our GK simulation. While in the context of reduced (z invariant) GK turbulence [92], the departure from scale invariance for an energy flux was shown to be dependent on the perpendicular collisions, the scaling displayed by the electron flux should be analysed from the perspective of Landau damping, as electrons dissipate energy mainly due to parallel velocity collisions and a v , z coupling is established via linear phase mixing. This acts as a reminder that phase space dynamics cannot be ignored, even though such an analysis goes beyond the scope of this paper and is left for future work.
Since definitions listed for a continuous system can sometimes offer a clearer understanding, we provide an equivalent definition for the scale flux across the waveband boundaries k c in terms of the infinitesimal triple-scale transfer S(k|p|q), for any a and b sub-domain limits. Unlike energy transfers that depend on the thickness of the bands, the value of Π(k c ) is identical regardless of it being computed via eq. (37) or from eq. (38).

G. The locality functions
Knowing the scale flux through k c , the infrared (IR) locality function is defined by taking a probe wavenumber boundary k p , so that for k p ≤ k c we have The IR locality function measures the contribution to the flux through k c from couplings with at least one scale wavenumber less than k p . In the second term, the sum over waveband Q starts from p + 1 to avoid double counting. In the limit k p → k c , we recover the flux across the cutoff wavenumber k c , i.e.
for any a, b indices that identify a set of bands. It is customary to normalise the locality functions to the flux trough k c , in which case a value of one is obtained for k p = k c and less than one for k p /k c < 1.
Although the IR functions have a clear interpretation as the ratio of energy contributed to the flux through scale k c coming only from larger and larger scales, it should be remembered that for k p /k c ≪ 1 the transfers can only take place between triads with one wave vector leg much smaller compared to the other two. Therefore, these functions can provide information regarding the overall locality of the nonlinear interaction. The rate with which Π ir (k p |k c )/Π(k c ) decreases in value as a function of k p /k c measures the locality of interactions. Larger exponents denote more local behavior while a zero exponent would qualify turbulence as fully nonlocal, i.e. every scale influencing equally the coupling of any other scale in the system. For our simulation data, the IR locality is presented in figure 7. The 1/12 value for the ion scaling exponent is consistent with the value reported by Ref. [69] in the first study of the locality  Using as reference the 4/3 value found for the asymptotic locality exponent in classical turbulence [59,63], we will associate larger values for the GK exponents to local and smaller values to nonlocal interactions. While this is an arbitrary choice, this simple characterisation will allow us to simplify the presentation of the results.
A similar definition is made for the ultraviolet (UV) locality functions, which for k c ≤ k p is given as It measures the contribution to the flux through k c from couplings of scales with at least one scale wavenumber greater than k p , therefore providing information regarding the locality makeup of a scale k c in relation with smaller and smaller scales. For our simulation the UV locality is presented in figure 8. A strong UV locality character is inferred from the large values of the UV locality exponents. We mention that selecting a cutoff k c ρ i = 0.55 gives the same scaling as for k c ρ i > 1, denoting that gyro-averaging effects don't influence the contributions to the flux across a scale emerging from interactions with progressively smaller scales. For both plasma species, a robust 5/2 value can be inferred for the UV locality exponent. For completeness, as we will use them in later sections, we provide the definitions of the locality 10 -2 10 -1 10 0 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 functions in terms of the infinitesimal triple-scale transfer S(k|p|q), These definitions are analogous to their waveband discrete forms, up to a manipulation of the IR integral limits similar to the one performed in eq. (41) and by employing eq. (39).

H. Reviewing the significance of the results
The results listed for the transfers, fluxes and locality functions have been, in one form or another, presented in the literature. We take the time to comment on their significance, before performing a more detailed analysis on the transfers.
The fact that the energy transfers are dominant between neighboring wavebands denote the local character of the cascade. However, considering that we take wavebands separated by less than a factor 2 (i.e. 2 1/3 ), the observed energy exchanges imply strong intra-dyadic interactions, known to lead to an enhanced nonlocal characteristic for the nonlinear interactions (we will detail this statement in section VI). This assertion is validated by the strong nonlocal character captured by the IR locality functions (i.e. exponents reaching values of 1/12 for the ions and 1/3 for the elections, values much smaller than the 2/3 and 4/3 found in MHD and hydrodynamical turbulence). The locality nature of interactions is important when modeling turbulence [93], but it also has a strong phenomenological significance. A strong IR non-locality is a sign of a very large scale flow shearing small scale structures, rather than advecting them. For example, a zonal flow [94] is expected to enhance the non-locality nature of plasma turbulence. On the other hand, the local UV behavior can be seen as an insensitivity of the perpendicular interactions to the type of collisional operators employed. This is fortunate from the perspective of turbulence modeling, even more so when we consider that both ions and electrons recover the same robust 5/2 value for the UV locality exponent, although their phase space dissipation route is found to be different [85].
In the case of IR locality functions (figure 7) we also observe a change in slopes at k p ρ i ∼ 1, an effect we consider specific to magnetised plasma turbulence, as turbulence above and below the ion gyroradius is expected to have different properties. The IR exponent seems to be much larger (more local behavior) at scales kρ i < 1. The ions strong nonlocal behavior occurs at scales kρ i > 1. As IR non-locality increases, the transfers associated with intra-dyadic interactions become more important in the energy cascade. As these intra-dyadic transfers, in the limit q → 0, can be seen as occurring between neighboring waves rather than compact structures, the GK turbulence cascade may have a stronger wave characteristic [76] compared to strong classical turbulence or even MHD turbulence.
For the remainder of this paper, we seek to understand better the link between dyadic separated couplings and intra-dyadic interactions and their contribution to the locality problem.

V. IDEALISED ENERGY TRANSFERS BETWEEN SCALES
While diagnostics based on S(K|P |Q) (e.g. the locality functions) can offer the most amount of information pertinent to the nonlinear interactions, the triple-scale transfers can be very demanding to compute, especially for the five-dimensional gyrokinetic problem. For a scale decomposition that uses N wavebands, the triple-scale transfers require N 2 calculations of the nonlinear terms. By comparison, computing the scale-to-scale transfers P(K|P ) directly from the nonlinear term requires only N nonlinear terms calculations. Since P(K|P ) measures the energy exchange between two scales, a natural question emerges: can we recover information pertinent to scale locality directly from P(K|P )? And if yes, what is its interpretation?

A. Definition
To address these questions and to build confidence in the interpretation of diagnostics that are applied in practice to turbulence systems that don't exhibit clear inertial ranges, we define a test case that will prove helpful in these matters. We choose: where α ir and α uv are here the two control parameters, together with the implicit choice for the wavebands. As the indices' suffixes imply, we take the α ir and α uv exponents to be related to the infrared and ultraviolet locality exponents. Below, we use the same waveband decomposition as the one given in section III D for the GK problem (i.e. λ = 2 1/3 ) and we use α ir = α uv = 5/2, to start. The resulting P test (K|P ) transfers are presented in figure 9. Not surprisingly, we recover an idealised forward cascade that is scaleinvariant (excepting the start and end of the band interval considered, due to numerical truncation effects). The scale invariance of the transfers is best seen from the figure 9-c) panel, where a perfect collapse of the scale-to-scale transfers as a function of P − K is observed. Moreover, from the same panel we observe that the "tails" of the transfers decreases gradually at ±(P − K).
In figure 10 we plot the absolute value of P test (K|P ) as a function of k P for all K values. We clearly see that for this simple test the tails follow a power law. We mention that from these type of plots, especially when considering our graphical representation, we only seek to identify an overall power law for the tails, rather than analyse individual transfers characteristics. The slope of the tails for K > P and K < P are indicated for reference, recovering the values prescribed for α ir and α uv parameters. While this should not come as a surprise due to our choice for P test (K|P ), we want to recover these exponents by means of locality functions. We emphasise the transfer for K = 14 as a thick line, to help visualise the typical K transfer curve. The slope of the tails for K > P and K < P are indicated for reference, recovering the 5/2 value prescribed for α ir and αuv.

B. Capturing the locality of the test transfers via locality type functions
Taking into account that we can recover from P test (K|P ) the flux through a surface k c (displayed in figure 11) by simply employing the definition given by eq. (37), i.e.
where N P =c+1 N K=c+1 P test (K|P ) = 0 due to energy conservation, we take an additional probe surface that limits the separation between the energy giving and the energy receiving scales. In the absence of a mediator, we cannot recover the IR and UV locality functions previously introduced and use the following definitions instead: We see that in the limit k p → k c , the two locality functions recover the energy flux given by eq. (47). The test locality function captures the prescribed exponents, figure 12. This is also seen from the two additional cases presented in figure 13, where we use the same value α uv = 5/2 but different values for α ir . We see that the values prescribed are recovered by the test locality exponents. We mention that the UV plots are exactly the same as the one listed in figure 12.
A few lessons can be drawn from this simple test. First, this simple test allows us to clearly identify characteristics we associate with asymptotic locality, i.e. for different values of k c all curves collapse on each other and exhibit the same power law. Second, for more nonlocal scalings (small value of the locality exponent), the locality functions fall off the asymptotic scaling at small k p /k c ratios. This fall-off is due to the high non-locality nature not being able to isolate the finite domain fast enough and limits our ability to gauge the correct IR locality behavior close to the largest scale. The slope at high k p /k c values is where asymptotic values should be investigated. We mention that for α ir = α uv the test transfer considered does not sum up to zero. This is a simple particularity of the definition employed, as there is no a-priori requirement for IR and UV locality exponents to be identical. This is acceptable as we just use the different α ir values to test the ability of the modified locality functions to capture the correct IR exponents.
The most important lesson to be drawn is also the most obvious: The locality exponents can be determined directly from the scale-to-scale transfer. However, these are not the exponents related to all interactions, but related to the locality of the energy cascade. We go back to the full GK simulation to show this fact.

VI. DETAILED MEASUREMENTS OF THE ENERGETIC EXCHANGES IN GK TURBULENCE
Compared to the idealised test case that was presented in the previous section, when investigating nonlinear energetic interactions and their scale locality in turbulence, we need to account implicitly or explicitly for the contribution made by the mediator scale. As an example, we start by looking at the ion scale-to-scale transfers for various separations between the energy exchanging scales.

A. The impact made by the separation of scales on the ion scale-to-scale transfers
The scale-to-scale transfers (eq. 31) are obtained by integrating over all mediator scales. Here, rather than doing so, we will characterise the exchanges between two scales p and k as a function of the values taken by the mediators q, and integrate accordingly over q to obtain a conditional form of the scale-to-scale transfers.
The minimal value for a scale q that mediates the interaction of two other scales (p and k) is given by eq. (21), i.e. q ≥ max{k, p} − min{k, p}. Expressing the value of max{k, p} in terms of min{k, p}, as max{k, p} = α min{k, p} with α ≥ 1, we obtain We see that the lower limit of the mediator q (i.e. q min = (α − 1) min{k, p}) and the minimal separation between p and k are linked. For k = p, we in fact have α = 1 and q min = 0. Imposing q ≥ q min = 0 as a selection condition would pick up all possible scale couplings for which max{k, p} ≥ min{k, p} (i.e. any p and k would lead to a q which is larger than q min ). This includes the k = p = q couplings. As another example, choosing max{k, p} = 2 min{k, p}, which means α = 2, leads to q min = min{k, p}. Imposing the selection condition q min = min{k, p}, selects all scale couplings (q ≥ q min ) for which max{k, p} ≥ 2 min{k, p}. Imposing the lower limit of the mediator q sets the the minimal separation possible between the energy exchanging scales. We see that in principle, α can act as a control parameter. For α = 1, which recovers k = p, we have q min = 0. Taking α = 2 leads to q min = min{k, p} and ensures that a dyadic separation (i.e. max{k, p} = 2 min{k, p}) is the minimal separation available between the energy exchanging scales. Going further with our example, for α = 3 we obtain max{k, p} = 3 min{k, p} as the minimal separation between energy exchanging scales and q min ≥ 2 min{k, p}.
Since our scales are prescribed as multiples of λ = 2 1/3 , e.g. k Q = k 1 λ Q−1 , the values that α can take are discrete and considered here as α = λ n + 1 for n ∈ Z. Considering that the Q index identifies the waveband s Q = (λ −1 k Q , k Q ], we write the condition given by eq. (50) for the waveband boundaries (i.e. q Q ≥ λ n min{k K , p P }) in terms of the waveband indices as Q ≥ min{K, P } + 1 + n .
The +1 emerges due to our choice for the waveband index (i.e. (k Q−1 , k Q ]) and the integer n can now be used as the control parameter for our discrete waveband representation. Corresponding to α = 2 we have n = 0 and Q min = min{K, P } + 1, while for the previous α = 3 example we have n = 3 and Q min = min{K, P } + 4. Not only that we can recover a minimal dyadic separation between energy exchanging scales for n = 0, but we can also select a minimal separation that is less than dyadic for n < 0. Accounting for the condition expressed by eq. (51), we define the modified scale-to-scale transfers as and plot them in figure 14 for different values of n, using the ion species as an example. We see that for a separation between energy exchanging scales that is at least dyadic (n = 0) the index separation between the intense peaks and the diagonal is 3 (since 2k Q = k (Q+3) ). As expected, selecting larger values for n limits the selection of energy exchanges between scales that are farther and farther apart. For clarity, we present in figure 15 the same transfers, but only for a selected receiving waveband K = 14. We see that only for n = 0, the scale k K is dyadic separated on each side from k P , i.e. from the black diamond we need to count 3 diamonds (due to λ = 2 1/3 ) to reach the either left or right peaks. Compared to the full scale-to-scale transfers that are dominated by the intra-dyadic exchanges, the dyadic separated transfers have a wider k support and decay to zero more gradually. This is important, since the decay of the transfers for P ≪ K and P ≫ K contains information relevant to the locality of the cascade, as we saw for the ideal case. We also note that for a larger separation between the energy exchanging scales, the transfers decrease in intensity, a fact seen from the diminishing value for the ratio between the maximal value of the transfer and the overall transfer norm ε. The results presented in figures 14 and 15 validates our expectation that limiting the minimal value of the possible moderators q results in an effective limitation of the minimal separation between the exchanging scales p and k. From the n = −6 case (for which q min = 1 4 min{k, p} and the minimal separation is max{k, p} = 5 4 min{k, p}) we can infer an additional important result. We see that we allow for a relatively small separation between the k and p exchanging scales and since q ∈ [ 1 4 min{k, p}, 6 4 min{k, p}], the q ∼ k or q ∼ p values are possible for the mediators. We can consider this as a proxy for the intensity of the energy exchanges for the k ∼ p ∼ q subset of couplings. By comparison, the full case (i.e. q ≥ 0), while allowing for these exchanges to take place, it also accounts for the exchanges between k ∼ p mediated by q ≪ k ∼ p. In figure 15, comparing the intensity of the full transfers with the ones for the n = −6 case, we see (from the maximal value used for normalisation in units of the transfer norm ε) that the latter (i.e. n = −6) are ten times smaller in value while exhibiting the same overall location for the peaks. This shows that the dominant intra-dyadic exchanges are the ones mediated by small values of q (i.e. q ≪ k ∼ p), which are intrinsically nonlocal from the perspective of the locality of interactions.

B. Extracting locality information from the scale-to-scale measurements
We seek to extract information pertinent to the locality problem from the scale-to-scale transfers. For ions and electrons, we look in figure 16 at the absolute value of P(K|P ) as a function of k P for all K values. While compared to the ideal (test) case (section V) we see that the intra-dyadic exchanges lead to an abrupt fall-off of the transfer curves (also seen in figure 15, full panel) and we see a large and small scale pollution, we are able to identify (IR and UV) slopes for the tails. In fact, looking at the P −6 (K|P ) case in figure 17, while noting that the missing points are due to our inability to compute eq. (52) for them, we see that the abrupt fall-off is removed while the tail information remains identical. The tail information, or its scaling to be more exact, is what determines the locality information. We also mention that, while P −6 (K|P ) removes the abrupt fall-off, only once the transfers occur between scales more separated (close to a dyadic separation) do we see them fall on a (IR or UV) power law (this can be inferred from the inspection of the thick black lines in figure 17). This is why the dyadic separation represents a useful reference. While not crucial, as the tails can be recovered from the full scale-to-scale transfers, they help identify the correct scalings.
In principle, a comprehensive analysis of the tails could be performed. Here, we resume to simply identify the IR and UV slopes. For both ions and electrons we identify a 5/2 slope for the transfers towards smaller scales (UV transfers). Regarding the exchanges with the larger scales (IR transfers), we identify a 4/2 scaling for the ions and a 3/2 scaling for the electrons. These scalings can be identified as the locality exponents for the energy cascade, which is contained in the full locality problem, i.e. the locality of interactions. We try next to recover these scalings from a set of modified locality functions.

C. Modified locality functions
We proceed now to impose the q ≥ q min condition for the locality functions, where we take q min = min{k, p}. Starting from the locality functions given by the eqs. (44)(45) in terms of the infinitesimal triple-scale transfer S(k|p|q) and accounting for the q ≥ min{k, p} condition, we define  the modified locality functions as In terms of our waveband decomposition, the equivalent definitions used in actual computations are given as We look at the IR in figure 18 and UV in figure 19 for the modified locality functions considered here. The locality curves show a clear power law and they collapse on each other for different cutoffs (k c ). This is a clear sign that we recover asymptotic locality exponents. Furthermore, using these modified locality functions, we find a 4/2 value for the ion's IR exponent, a 3/2 value 10 -2 10 -1 10 0 10 -9 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 for the electron's IR exponent and a 5/2 value for both UV exponents. These are the same values as the ones found for the locality exponents of the cascade. As such, we see that we can recover these exponents directly from the scaling of the tails of the scale-to-scale transfers. However, we still need to relate these asymptotic locality exponents (for the energy cascade) with the ones in section IV G (the exponents for the locality of interaction). We do so next and make our final conclusions.

VII. DISCUSSION AND CONCLUSIONS
A. Discussing the relation between the IR locality function and its modified form We first mention that in the UV limit, the locality exponents for the nonlinear interactions (obtained via the normal locality functions) and the asymptotic locality exponents found for the energy cascade (either directly from the scale-to-scale transfers or via the use of the modified locality functions) recover the same 5/2 value for both ions and electrons.
Compared to the UV case, the IR locality exponents obtained for the full and modified definitions of the locality functions differ drastically. To understand better the significance behind this difference, we plot in figure 20 the (p, q) integration domain of S(k|p|q) for eq. (44) and its two constituent terms. Since in the S(k|p|q) object the position of the k, p, q scales matters, the two terms in eq. (44) have different physical interpretations.
Due to the definition of IR locality functions, the ordering k ≥ k c ≥ k p is always valid. Considering now solely the q ≥ p (first) term, for which k p ≥ p, we see from the condition listed in eq. (21) (or equivalently in eq. 50) that we always have at least a dyadic separation between the giving and receiving scales (i.e. k ≥ 2p) for k c ≥ 2k p . As explained in section III C, such interactions are mediated only by immediate dyadic structures above and below the energy receiving scales. Thus, these interactions are quite local and as a consequence scale very well with the energy receiving scale k, i.e. the dynamics of these interactions tend to become scale invariant and recover asymptotic locality properties. The separation imposed by k p and k c limits now the separation between the giving and receiving scales and the locality exponents thus measured are directly related to the locality of the energy cascade. The modified IR locality functions stated in eq. (53) is simply this first term, which measures the intensity of the energy cascade as a function of the separation between the energy exchanging scales and leads to the recovery of asymptotic locality exponents. We mentions that for k c = k p the q ≥ p conditions can select intra-dyadic exchanges, however, these exchanges tend to be mediated by q ∼ k and do not increase the nonlocal nature of interactions. Moreover, interactions of the type q = p = k cannot be seen as moving energy from one scale to another.
In the second term, we have p ≥ q and k ≥ k c ≥ k p ≥ q. This limitation on the mediator scale q selects only transfers of energy between scales contained in the same dyadic signal. While the energy exchanges are indeed very local, the largest contributions involve couplings for which the mediator is well separated from the energy exchanging scales q ≪ p ∼ k. This leads to a strong nonlocal character for the nonlinear interactions, which is captured by the second term's exponents. The fact that the mediation of intra-dyadic transfers p ∼ k is done by scales comparable to the forcing range, scales which are thus amplified directly by the linear forcing term, is listed as a warning to the impact made by a force on the development of turbulence [66]. The full IR locality function captures both of these distinct aspects and, depending on their intensity, provides an effective locality exponent. As we have shown, asking for k and p to be at least dyadic separated leads to a null contribution for the second term. This fact emphasises the interpretation that the second term measures the locality of plane waves interacting inside the same dyadic signal and its not a measure of the classical turbulence energy cascade (overall, the energy exchange between a dyadic scale and itself is zero). Indeed, for hydrodynamics and MHD turbulence, the contribution of the second term is sufficiently reduced as to recover globally the asymptotic locality exponents. This is not the case for GK turbulence, where the exchange between neighboring waves seems to remain strong even at the smallest scales.

B. Discussing the asymptotic locality exponents
For GK turbulence, the effective IR locality exponents do show that turbulence is strongly nonlocal. However, at the same time, our detailed analysis shows that the energy exchange between well separated scale structures recovers asymptotic values for the locality exponents. In the past [69], a value of (k p /k c ) ±5/6 was estimated for the IR and UV asymptotic locality exponents. Those values were determined using scalings computed for statistically homogenous two-dimensional GK turbulence [27]. Here, we do not confirm those predictions, as we find, respectively, 4/2 and 3/2 values for the ion and electron IR asymptotic locality exponents. For the UV, both ions and electrons converge on the 5/2 value for the asymptotic locality exponent, which are recovered from the full form of the locality functions, denoting their robustness.
The recovery of asymptotic locality exponents is indeed impressive, as it shows that in spite of all existing complications, GK turbulence possesses a strong classical characteristic for the kinetic Alfvén wave cascade. By filtering the large scales mediators and limiting the contributions of intra-dyadic exchanges we allow the nature of the remainder interactions to surface. We see that embedded in the full GK turbulence problem, there exists an asymptotic turbulence component (just like hydrodynamical or MHD turbulence), which under ideal conditions may yet be realised. Even if only a tendency, the asymptotic locality nature of GK turbulence can still be measured and the least expensive way is through the scaling of the scale-to-scale transfers. In the future, measuring these asymptotic locality exponents for different plasma parameters in simple magnetic configurations will tell us if this classical characteristic is also universal for GK turbulence.
For GK turbulence to recover a general asymptotic behavior, i.e. for its nonlinear dynamics to remain scale invariant with the increase of the interval of excited scales, the contribution of intra-dyadic exchanges (mediated by large scales) should tend towards zero. In this scenario, we would recover the asymptotic locality exponent directly from the locality functions (as the second term, p ≥ q, would become sub-dominant) and the locality of the cascade and the locality of the nonlinear interactions would become the same. If furthermore we find the asymptotic locality exponent to be unique, we could say that GK turbulence has a universal character. However, even if the intra-dyadic exchanges will always dominate the GK system, the values of the asymptotic locality exponent can still be determined from the energy cascade and their uniqueness be assessed for various plasma parameters.
Last, we mention that all theoretical estimates that assume an infinite inertial range in turbulence will automatically assume that the large scales mediated intra-dyadic exchanges are zero (or infinitely small), since the largest scales that couple in a nonlocal way will always be removed by the infinite range limit (also seen as: the large shearing flows are removed and homogeneity is restored for turbulence). This last statement can be interpreted as a warning, not to over-rely on simple scaling laws for turbulence in complex systems.

C. Conclusions
Using a large resolution simulation of GK turbulence in slab magnetic geometry, we have analysed the energy cascade and the locality of interactions. While the interactions can be deemed as being nonlocal (1/3 scaling for electrons and 1/12 for the ions in the IR limit), we have shown that embedded in the full GK problem, there exists a set of couplings that tend to develop an asymptotic turbulence behavior. These couplings possess locality exponents that can be recovered directly from measurements on the energy cascade. The energy cascade was shown to be local in nature (3/2 scaling for electrons and 4/2 for the ions in the IR limit) and, more importantly, it was shown that it recovers asymptotic values for the locality exponents. This may prove to be useful in the development of sub-grid scale models for GK turbulence.
In addition, clarifying the diagnostics that can capture the asymptotic exponents for GK turbulence is important as a wide parameter space needs to be explored to evaluate the universality of turbulence at kinetic scales. Being able to extract the asymptotic locality exponents from the less expensive scale-to-scale transfers can prove to be invaluable, especially when seeking more robust time averaged results.
From the perspective of the five-dimensional dynamics, understanding that the exchange of energy between perpendicular scales occurs in an asymptotically local way is important when seeking to understand the balance between nonlinear phase mixing (that includes the perpendicular cascade) and linear phase mixing (Landau damping) using turbulence scaling arguments. This is left for future work.
Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.