Kinetic plasma turbulence: recent insights and open questions from 3D3V simulations

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

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

DATA SETS
In the following, we consider three recent kinetic simulations in a six-dimensional phase space ("3D3V") using: (i) CAMELIA, a hybrid particle-in-cell (PIC) code with massless electrons (Franci et al., 2018a), (ii) HVM, an Eulerian hybrid-Vlasov code with finite electron-inertia effects (Valentini et al., 2007), and (iii) OSIRIS, a fully kinetic PIC code (Fonseca et al., 2002(Fonseca et al., , 2013. Unless otherwise specified, parallel ( ) and perpendicular (⊥) directions are defined with respect to the global mean magnetic field B 0 = B 0 e z . Franci et al. (2018b) employed the CAMELIA code to investigate freely decaying, Alfvénic fluctuations in a cubic box (L = L ⊥ = 128d i with 512 3 grid points and 2048 particles per cell (ppc)) for β i = β e = 0.5, where β s = 8πn 0 T s /B 2 0 is the species beta. Cerri et al. (2017b) instead adopted the HVM code to study freely decaying compressive fluctuations in an elongated box (L = 2L ⊥ ≃ 63d i with 384 2 × 64 grid points in real space, and 51 3 points in a velocity space bounded by |v/v th,i | ≤ 5) for β i = β e = 1 and with a reduced ionelectron mass ratio of m i /m e = 100 (viz. including d e -effects in a generalized Ohm's law). Spectral filters were applied at runtime, determining a cutoff in the turbulent spectrum at k ⊥ d i > 20 and at k z d i > 2. Finally, Grošelj et al. (2019) use the OSIRIS code to investigate continuously driven Alfvénic fluctuations in a β i ≈ β e ≈ 0.5 plasma with m i /m e = 100. An elongated box was used (L = 2.56L ⊥ ≃ 48.3d i with 928 2 × 1920 grid points and 150 ppc per species). An example of δ B ⊥ = δB ⊥ /δB (rms) ⊥ in a two-dimensional cut perpendicular to B 0 is given in Figure 1A, along with a schematic representation of these simulations in the (k ⊥ , k ) plane ( Figure 1B).
In the following, the analysis of freely decaying simulations (viz., CAMELIA and HVM) is performed at the peak of the turbulent activity (cf., e.g., Servidio et al., 2015), while for the continuously driven OSIRIS run we consider the turbulence at the end of the simulation when the kinetic range spectra appear converged. Following Franci et al. (2018b) and Grošelj et al. (2019), PIC data have been filtered before performing the analysis to remove spectral regions dominated by particle noise. The OSIRIS data have been filtered for k ⊥ d i > 30 or k z d i > 12 and downsampled to a grid 464 2 × 640. Note that OSIRIS simulations require to resolve the Debye scale, while the physical scales of interest are well represented at a lower resolution. A short-time average over t ce = 2 ( ce being the electron cyclotron frequency) was also performed to further reduce electron-scale noise . The CAMELIA data have been filtered for k ⊥ d i > 10 or k z d i > 2. We also considered alternative filtering approaches confirming that our results are not very sensitive to such particular choice.
In Figure 1C the two-dimensional Fourier spectra, E(k ⊥ , k z ), are shown. The wavenumber region (k ⊥ , k z ) occupied by the turbulent fluctuations already highlights the anisotropic nature of the cascade, with energy preferentially flowing to high k ⊥ . However, note that the 2D Fourier spectrum may exhibit a weaker anisotropy than the one typical of turbulent eddies, which are elongated along the local field direction (see, e.g., Cho and Vishniac, 2000). We perform a local analysis of anisotropy in section 4.1.
In Figures 1D,E, the reduced 1D spectra, E(k ⊥ ) (upper panels), and their local slope (lower panels) are reported. To remove the effects of different energy injection conditions, the k ⊥ -spectra have been normalized so that they overlap in the subion range, at k ⊥ d i ≃ 5. According to the spectral anisotropy in Figure 1C, the k z -spectra have been consequently matched at k z d i ≃ 0.5. For our choice of low-pass filter (see section 2), CAMELIA spectra artificially flatten beyond k ⊥ d i 7 due to PIC noise, and therefore we do not show CAMELIA data in this range in Figure 1D. Overall, the spectral slopes are consistent with each other, although the spectra obtained from the three simulations do not quite assume a universal shape. Close to the box scale, the spectral exponents are likely affected by the turbulence injection details. It is also possible that some of the sub-ion scale spectral exponents are not fully converged in terms of the box size (which was different for each simulation) and of the limited extent of sub-ion range itself. 3D3V simulations with a significantly larger sub-ion range are required to clarify this point. To some degree, differences at sub-ion scales could also be physical. In particular, the HVM simulation includes electron inertia effects in Ohm's law, while the OSIRIS results include the full electron kinetics, such as electron Landau damping and finite electron Larmor radius corrections. It is interesting to notice that OSIRIS spectra become steeper than the hybrid counterparts beyond k ⊥ d i 3, for our particular choice of the mass ratio (m i /m e = 100). This feature has been usually explained in terms of electron Landau damping (Grošelj et al., 2017;Chen et al., 2019), which is not included in the hybrid-kinetic model.
In Figure 1F, we report the comparison of spectral ratios, C 1 δB 2 z /δB 2 ⊥ (top), C 2 δn 2 /δB 2 ⊥ (middle), and C 3 δn 2 /δB 2 z (bottom). The ratios are normalized to the β-dependent kinetic Alfvén wave (KAW) eigenvalue from asymptotic linear theory and k ≪ k ⊥ ), namely C 1 = (2 + β)/β, C 2 = (2 + β)β/4, and C 3 = β 2 /4, where β = β i + β e (see, e.g., Boldyrev et al., 2013, for details). In the normalized units, asymptotic KAW theory predicts a value of unity for all three ratios. This is essentially the result of KAWs developing a certain degree of magnetic compressibility at sub-ion scales, which sets the relation between δB ⊥ and δB , and requiring that compressive magnetic fluctuations are pressure balanced, which in turn provides a relation between δB and δn (see, e.g., Schekochihin et al., 2009;Boldyrev et al., 2013). As found in previous studies (e.g., Salem et al., 2012;TenBarge et al., 2012;Chen et al., 2013;Cerri et al., 2017b;Franci et al., 2018b;Grošelj et al., 2018), the spectral field ratios are overall consistent with KAW-like turbulence at sub-ion scales. This is not completely surprising, as both in CAMELIA and OSIRIS simulations, Alfvénic fluctuations are injected. On the other hand, compressible magnetic fluctuations (i.e., including δB ) are injected in HVM run, and yet KAW-like fluctuations still develop. It was also proposed that KAWs may, quite generally, emerge as a result of Alfvén waves interacting with large-scale inhomogeneities (Pucci et al., 2016). Thus, the KAW-like spectral properties at sub-ion scales appear to be a relatively robust feature, independent of the details of the turbulent fluctuations injected at the MHD scales (cf. Cerri et al., 2017a). While the results are overall consistent, some differences are also seen, most notably in the high-k ⊥ range (k ⊥ d i 10), which could be presumably attributed to various numerical artifacts. However, some deviations could also relate to differences between the hybrid-kinetic and fully kinetic model [for instance, some dispersion relation properties not being exactly the same (e.g., Told et al., 2016)].
So are the sub-ion-scale field polarizations indeed KAW-like? As discussed above, recent observations and kinetic simulations are consistent with such idea, although linear wave predictions are not necessarily satisfied precisely (e.g., Chen et al., 2013;Kiyani et al., 2013;Cerri et al., 2017b;Franci et al., 2018b). Chen et al. (2013) report an average value of 0.75 for the normalized ratio C 2 δn 2 /δB 2 ⊥ , whereas (asymptotic) KAW theory predicts a value of unity. That latter may be due to different reasons, among which we remark the following two: (i) sub-ion-range turbulence is not made of purely KAW-like fluctuations, and/or (ii) the asymptotic conditions that are used in the derivation of linear theory predictions are not met exactly because of the limited subion range of scales and/or because of the inherently non-linear dynamics of turbulence. These two explanations are not mutually exclusive, of course. Indeed, sub-ion-scale turbulence can in principle include contributions from wave-like fluctuations of other nature. This may include fluctuations consistent with whistler (e.g., Gary and Smith, 2009), ion-cyclotron (e.g., Omidi et al., 2014;Zhao et al., 2018), or ion Bernstein waves (e.g., Podesta, 2012;Del Sarto et al., 2017;Grošelj et al., 2017), to name a few. On the other hand, the spectral ratios could also deviate from linear KAW predictions as a result of non-linear dynamics.
For example, Boldyrev et al. (2013) propose that, specifically the (normalized) C 2 δn 2 /δB 2 ⊥ ratio may fall somewhat below the KAW prediction due to a (yet to be investigated) non-linear effect, analogous to the residual-energy phenomenon in MHD turbulence.

MULTI-POINT STRUCTURE FUNCTIONS
Beyond energy spectra, fluctuations across different scales may be investigated in more detail via structure functions, i.e., the moments of local field increments (e.g., Frisch, 1995;Biskamp, 2008). Two-point structure functions, S (2) m (m being the order), are most common. However, these cannot quantitatively produce the correct scaling for fluctuations with power spectra steeper than ∼ k −3 , assuming a clean power-law spectrum (Falcon et al., 2007;Cho and Lazarian, 2009). Therefore, structure functions using more than two points are generally required at kinetic scales. Essentially, higher-order increments yield a scale decomposition that is more effective in filtering out the large-scale fluctuations below k ≈ π/ℓ in spectral space, where ℓ is the increment scale. We also mention that if the signal is a polynomial of degree N − 2, its corresponding 2ndorder, N -point structure function vanishes (Cho, 2019). This makes multi-point structure functions more suitable for the analysis of relatively smooth signals with steep spectra (Schneider et al., 2004). A detailed review of N -point increments, as well as their physical interpretation can be found in Cho (2019). Here, we consider for some field f (x) the conditional, five-point structure functions: . . x is a space average, and ϑ B loc is the angle between the increment vector ℓ and the local mean magnetic field B loc . The term "conditional" implies that S m are defined as conditional averages of | f (x, ℓ)| m , using only those points in the statistical sample that fall within a given (narrow) range for ℓ and ϑ B loc . We also considered three-point structure functions (see Figure 2A) and, for a limited number of cases, seven-point structure functions (not shown). Comparison between the three-point, five-point and seven-point structure functions shows not only qualitative similarities among the three cases, but an apparent quantitative convergence with increasing number of points. We chose to illustrate the results in Figure 2 in terms of five-point structure functions in order to provide better constraints for the theoretical predictions. Similar to two-point increments, where the local mean field is often defined as B loc (x, ℓ) = [B(x) + B(x + ℓ)]/2 (e.g., Cho and Vishniac, 2000;Mallet et al., 2016), we obtain B loc by averaging over the points used for the increment. For five-point increments, a reasonable definition is B loc (x, ℓ) = [B(x + 2ℓ) + 4B(x + ℓ) + 6B(x) + 4B(x−ℓ) + B(x−2ℓ)]/16. It is straightforward to check that such mean field definition filters out fluctuations around the scale of the increment ∼ ℓ, while preserving the contribution from scales larger than ℓ.
In what follows, we consider field-perpendicular, S m (ℓ ⊥ ) ≡ S (5) m (ℓ ⊥ , 90 • − ϑ ≤ ϑ B loc ≤ 90 • ), and field-parallel, S m (ℓ ) ≡ S (5) m (ℓ , ϑ B loc ≤ ϑ), five-point structure functions of the magnetic field and density fluctuations, where ϑ represents a finite angular tolerance used in practice to determine the local perpendicular and parallel directions. We reduce ϑ until the scalings appear converged. The field increments, from which we obtain the conditional structure functions, are evaluated at every grid point. In each grid point and at every scale, increments are sampled along random directions. The numbers of these random directions per grid point have been tested to provide a statistically significant (i.e., converged) sample. The sample that is used in the following is such that any structure function S m (ℓ, ϑ B loc ) counts at least 1.5 × 10 5 points per scale ℓ, in any given band for ϑ B loc .

Spectral Anisotropy
A delicate point concerns the sub-ion-range spectral anisotropy, k ∼ k ⊥ α (cf., e.g., Schekochihin et al., 2009;Boldyrev and Perez, 2012;Cerri et al., 2018;Landi et al., 2019). As is known from MHD, electron-MHD (EMHD), and kinetic-reduced-MHD (KRMHD) turbulence (Cho and Vishniac, 2000;Cho and Lazarian, 2009;Meyrand et al., 2019), the true anisotropy is often revealed only when measured with respect to the local, scale-dependent mean magnetic-field direction. Somewhat contradicting estimates, obtained with different methods, for the sub-ion-scale anisotropy have been presented in recent works. Here, we revisit this issue using the above-mentioned implementation of five-point structure functions, consistently applied to all data.
In Figure 2A we show the perpendicular and parallel secondorder structure function scalings, and in Figure 2B we show the inferred anisotropy, ℓ (ℓ ⊥ ). The characteristic parallel length ℓ (ℓ ⊥ ) at a given perpendicular scale ℓ ⊥ is obtained by finding the value of ℓ , at which the amplitudes of S 2 (ℓ ) and S 2 (ℓ ⊥ ) match. To illustrate the sensitivity to the local mean field direction, we show in Figure 2A the convergence with respect to the angular tolerance ϑ. The parallel scalings appear converged at ϑ ≃ 3 • for CAMELIA data and at around ϑ ≃ 1.5 • for HVM, whereas the OSIRIS results are somewhat less sensitive to ϑ (converging already for ϑ ≃ 6 • ). This difference may occur because OSIRIS simulation exhibits the weakest anisotropy (in absolute values). Physically, ϑ should be approximately no larger than ∼ ℓ ⊥ /ℓ of the small-scale turbulent eddies. Thus, smaller ϑ are needed if a stronger anisotropy develops at the energy-containing scales.
All quantities seem to converge to a scaling close to ℓ ∼ ℓ ⊥ 2/3 (although δB ⊥ fluctuations in HVM exhibit a scaling closer to 1/3 over the range of scales across ℓ ⊥ ∼ d i (= ρ i )). It is worth noticing, however, that this is not the end of the story, as the scaling is not quite 2/3 and additional effects such as B-field curvature may slightly change the anisotropy. Indeed, the field increments are taken along a straight line. If the local magnetic field lines are significantly curved over the extent of the increment stencil (= 4ℓ for five-point increments), the field increments will mix contributions from different field lines, in which case the anisotropy may be somewhat underestimated. It is worth mentioning that a scaling ℓ ∼ ℓ ⊥ 2/3 was proposed in 2 , of δB ⊥ (top row) and δB (bottom row) vs. ℓ ⊥ (continuous lines) and ℓ (dashed lines). Here, and ⊥ are defined with respect to the local field direction (i.e., δB = δB z ) with different angular tolerance, ϑ (colored lines; see text for definition). S A 1/ℓ ⊥ scaling is given for reference. Boldyrev and Perez (2012), based on a filling-factor correction for the fluctuation energy. Assuming the energy is concentrated in intermittent, two-dimensional structures as in Boldyrev and Perez (2012), the filling factor should scale as k −1 ⊥ ∼ l ⊥ . The filling factor may be approximately estimated from the inverse scaling of the excess kurtosis ; see section 4.2). Our results shown in Figure 2C are indeed roughly consistent with an excess kurtosis scaling ∼ ℓ −1 ⊥ , although this approximate scaling is overall better satisfied for δB and δn than for δB ⊥ . Finally, we mention that an alternative anisotropy estimate, based on a spectral band-pass filter (Cho and Lazarian, 2009), gives a somewhat stronger anisotropy than the five-point structure functions (not shown). On the other hand, qualitatively similar results are still obtained for all data. Thus, all simulations analyzed exhibit a similar sub-ion-scale anisotropy according to the particular diagnostics employed. Therefore, the differences that were previously reported in the literature could be mainly related to the different methods employed.

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

CONCLUDING REMARKS
So, what is the nature of sub-ion-scale fluctuations? From our independently performed 3D3V (hybrid and fully) kinetic simulations, a picture consistent with KAW turbulence phenomenology emerges. Moreover, our results imply a scaledependent anisotropy, together with intermittent statistics of magnetic and density fluctuations at sub-ion scales. Thus, we conclude that within the range of parameters explored here, the statistical properties of ion-scale plasma turbulence (at β ∼ 1) definitely show a certain degree of similarity, regardless of the precise details of the large-scale energy injection. On the other hand, slight differences can also be identified, some of which may be also model-dependent.
A number of key aspects will have to await the nextgeneration of 3D3V kinetic simulations. Ideally, future numerical experiments should aim to resolve both larger (MHD) scales, as well as a broader range between the ion and the electron scales by adopting significantly higher (if not realistic) mass ratios. These two aspects indeed appear to be both required in order to achieve (i) a possible saturation of the kurtosis at ion scales and (ii) a relevant sub-ion range of scales before electron-scale effects significantly come into play. Moreover, different aspects other than the spectral and statistical properties of the turbulent fluctuations will need to be considered in characterizing kineticrange turbulence, as for instance, the dissipation mechanisms of turbulent fluctuations under different plasma conditions and the consequent energy partition among different species (e.g., Matthaeus et al., 2016;Parashar et al., 2018;Arzamasskiy et al., 2019;Kawazura et al., 2019;Zhdankin et al., 2019).
While certain progress was definitely achieved in recent years, many other plasma regimes and setups may need to be explored, and the process(es) underlying a possible universality of kinetic-range plasma turbulence (e.g., magnetic reconnection) need to be fully worked out. Moreover, a few relevant discrepancies between numerical simulations, theories and in-situ observations appear. These "anomalies" definitely call for an explanation by the space physics community. In this context, advances cannot be achieved without investing in next-generation multi-spacecraft missions. Multi-point in situ measurements of turbulent fluctuations from a large number of spacecrafts are indeed fundamental in order to disentangle the non-linear spatio-temporal character of plasma turbulence (see, e.g., Klein et al., 2019;Matthaeus et al., 2019;TenBarge et al., 2019). This includes answering fundamental questions about, for instance, (i) the distribution of turbulent energy in space and time, (ii) the three-dimensional anisotropic structure of energy transfer across scales, (iii) the highorder statistics of the fluctuations, and (iv) the validity of Taylor's hypothesis over a broad range of time and spatial scales. Alongside observations, advances in computational capabilities are required to perform more realistic numerical simulations as discussed above, and compare these with spacecraft measurements. Finally, following the same spirit promoted by the "Turbulent Dissipation Challenge" , we would like to end this Perspective by stressing that our community could benefit from comparisons such as the one performed here, involving various codes, models and diagnostics.
Note added: Arzamasskiy et al. (2019) recently reported a scale-independent anisotropy at ion scales (i.e., ℓ ∼ ℓ ⊥ ) based on a set of 3D driven hybrid-kinetic turbulence simulations. Using our structure function diagnostic applied to their data, we were able to qualitatively (and quantitatively) reproduce their result. A more detailed investigation along these lines is currently ongoing, but beyond the scope of this Perspective and will be presented elsewhere.

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

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