Understanding the level of Turbulence by Asymmetric Distributions: a motivation for measurements in Space Plasmas

In this article, on the basis of the Langevin equation applied to velocity fluctuations, we numerically model the Partial Variance of Increments, which is a useful tool to measure time and spatial correlations in space plasmas. We consider a Coupled Map Lattice model to relate the spatial scale of fluctuations, $k$, to some macro parameters of the systems, as the Reynolds number, $R_\lambda$, the $\kappa$ parameter of Kappa distributions, and a skewness parameter, $\delta$. To do so, we compute the Velocity Probability Density Function (PDF) for each spatial scale and different values of Reynolds number in the simulations. We fit the PDF with a Skew-Kappa distribution, and we obtain a numerical relationship between the level of turbulence of the plasma and the skewness of obtained distributions; namely $\langle \delta \rangle \sim R_\lambda^{-1/2}$. We expect the results exposed in this paper to be useful as a tool to characterize the turbulence in the context of space plasma and other environments.


INTRODUCTION
Understanding why systems are turbulent and how these features are manifested usually leads to non-trivial concerns. Even more, if in this discussion we consider that plasma systems, in particular, present several of these characteristics, the complexity of studying plasma turbulence increases considerably. Therefore, the study of plasma turbulence starting from known facts of turbulence in other systems is customary. Currently, we know that plasmas and neutral fluids have many similarities when it comes to turbulence (Matthaeus & Velli 2011;Matthaeus 2021). To date, there are reports about the turbulent cascade in variables as velocity, electric and magnetic field (Bruno & Carbone 2013). There is also the record of the high space-time correlations (Matthaeus et al. 2005(Matthaeus et al. , 2016. And in this last topic there are many reports about the Partial Variance of Increments (PVI) of the fields in space plasmas. With this definition, time series can be constructed associating them to a characteristic time or spatial scale. Some authors studied these time series of velocity and magnetic field in space plasmas and they could observe some non-extensive statistical features (Pollock et al. 2018;Chasapis et al. 2018;Chhiber et al. 2020). This non-Maxwellianity of the statistics is mainly reflected in the distributions of the PVI, which present heavy tails out of the Maxwellian distributions.
In other words, distributions with a kurtosis greater than 3 (Greco et al. 2017).
Regarding non-extensive statistics, Tsallis entropy provide a framework to describe mathematically these distributions out of Maxwellian statistics (Tsallis 1988). This is possible through a new parameter, that enables to obtain fits with a good agreement with distributions that behave as power law at large energies. The family of Tsallis distributions have many forms to account for these properties. In the context of theoretical non-extensive statistics. Tsallis entropy typically works with q-canonical distributions (Wilk & W lodarczyk 2000;Picoli Jr et al. 2009). In plasma physics is usual to see Kappa distributions (Hasegawa et al. 1985;Yoon et al. 2006;Kim et al. 2017), represented by the so-called κ parameter, with κ = 1/(q − 1) (Pierrard & Lazar 2010). But in essence they are quite similar. However, recently, a theoretical framework has also been constructed from thermodynamic physics that manages to explain κ as a non-extensive parameter, as well as the temperature of a system (Livadiotis & McComas 2021). The main conclusion from the latter reference is that, in the same way that we relate the temperature to the average kinetic energy of the system, the kappa parameter is the result of the average correlations between the elements of a system.
In the literature we can find different models and approaches that address the issue of electromagnetic fluctuations in the context of turbulence in space plasmas. Some of these are based on particle (Biancalani et al. (Beck 1994;Hilgers & Beck 1997). Under this context, Beck (2000) obtained an analytic expression for the Probability Density Function (PDF). Among different properties of the obtained PDF, Beck (2000) showed how some turbulent systems naturally develop asymmetric distributions. In fact, that author argued that the asymmetry of the PDF is related with the level of turbulence of the flow, and obtained an analytical expression of the asymmetry as a function of the Reynolds number.
In the nature of space plasmas we can also find asymmetries in the velocity distributions. An example of that is the electron Velocity Distribution Function (eVDF) in the solar wind, such as Pilipp et al. (1987);Nieves-Chinchilla & Viñas (2008) have shown. This is a very relevant phenomenon that can introduce non-trivial changes the properties of the system. For example, Zenteno-Quinteros et al.
(2021); Zenteno-Quinteros & Moya (2022) have studied the dispersion properties of space plasma for similar conditions in which the solar wind is found. In their works they analyze the whistler mode and characterize the heat flux instability, as expected of a system that has a small asymmetry in its particle distribution. To do this, the authors compute the dispersion relation, with a Beck-like distribution, or as Zenteno-Quinteros et al. (2021) called, a Skew-Kappa distribution, namely: where A κ,δ is a normalization quantity, and w, is related to the thermal speed of the distribution. In addition, κ and δ are dimensionless parameters related with the high energy tails and skewness of the eVDF, respectively.
Under these context, recently Gallo-Méndez & Moya (2022) performed a Langevin approach of a CML, and contextualized its use as a tool to understand the fluctuations in a turbulent system with no skewness. The authors focused on the relation between the spatial scale and the κ parameter in Kappa distributions, and obtained an empirical scaling relation between the level of turbulence and κ; namely κ ∼ R λ k −5/3 , where k represents the spatial scale, and R λ the Reynolds number of the system, respectively. In the present article we expand the work by Gallo-Méndez & Moya (2022) to the case of turbulent systems exhibiting skewness, emphasizing on the relation between the level turbulence, quantified by the Reynolds number, and the skewness parameter, δ, in a velocity distributions as in (1). For such a purpose, we generate the steady state velocity distribution at each spatial scale, and we show that the generated distributions can be well modeled by a Skew-Kappa distribution. We organized this paper as follows: First, in section 2, we explain the basis of the model used in Gallo-Méndez & Moya (2022) but now using the Ulam Map as a chaotic driver force.
Then, in section 3, we expose the methodology for the treatment of the data obtained. Furthermore, we show how we modeled the relation between the spatial scale, k, the Reynolds number, R λ , the κ parameter, and δ parameter. Finally, in section 4, we summarize our findings and present the main conclusions.

MODEL: COUPLED MAP LATTICE
Following the model proposed by Beck & Roepstorff (1987) and Beck (1994), we consider a CML approach to model the eddies in a turbulent system as a Langevin type form; i.e., the force equation for each element is the following: where γ is the viscosity of the media, and ζ(t) is the noise forcing which serves as a driver. Here it is necessary to highlight that, even though there is no direct mention about the magnetization of the plasma in Eq.
(2), the effect of electromagnetic forces are indeed included through the different parameters and terms of the equation. Following Chun et al. (2018Chun et al. ( , 2019 the magnetic field can be included as part of the friction term γ. Similarly, the effect of the electric field can be introduced through the stochastic term ζ, as in Beck (1996). Furthermore, the model can also be used to study the nature of small scale magnetic fluctuations in a turbulent plasma. For example, in Carbone et al. (2021Carbone et al. ( , 2022 the authors propose a stochastic model based on the Langevin equation to describe the nature of kinetic scales magnetic fluctuations in the interplanetary medium, based on the principles of Brownian motion. Analogously to this study, although their model does take directly into account wave-particle interactions, their influence is indeed included through the parameters and statistical properties of the Langevin equation they analyze. In summary, Langevin-like models may be also used to describe the statistical behavior of plasma fluctuations, using a stochastic approach that takes into account the interactions among large an small spatial scales. Thus, this approach is a valid and useful way to study the behavior of plasmas, even if it does not explicitly include the effects of a background magnetic field or wave-particle interactions.
In Eq.
(2), the dynamic variable v is given by and denotes the difference of the longitudinal component of the velocity field u at two points separated by a characteristic distance r. Here, v is proportional to the PVI (Greco et al. 2017) for space increments, where r determines the observed spatial scale. In order to include the interaction between scales, the CML assumes that there is a fully developed turbulence state, and this state is characterized by the existence of eddies in different spatial scales. Here, we labeled these scales by an index k. For k = 1, the model contemplates eddies in the largest spatial scale, ℓ 0 /2, which are moving through a medium surrounded by other smaller eddies with size ℓ 0 /2 k in the inertial range.
We model the first scale by the Langevin force equation in (2), with a chaotic forcing given by Here, we interpret the characteristic time τ as the time average between "collisions" or interactions in the largest scale, which could be produced by an arbitrary driving force in the system, and x n is the normalized amplitude of the fluctuations or noise. In order to study the formation of heavy tailed and asymmetric VDFs, we consider that x n follows the chaotic Ulam map, x n+1 = 1−2x 2 n , which keep the fluctuations in the range of -1 to 1, with a probability distribution given by f As shown by Beck Beck (1994, 1996, the use of the chaotic Ulam map is necessary to model the chaotic forces that produces skewed VDFs in the flow. Since the smallest timescale is τ , it is convenient to discretize equation (2) using (4). Therefore, integrating the Langevin equation in a small time interval ∆t = τ , we obtain where we define v n as v(τ n) and λ = e −γτ . Further, to model the interactions between scales we consider the conservation of momentum between scales k − 1 and k. The basic assumption is that the momentum loss at level k − 1 serves as a driving force at level k. Moreover, daughter eddies at level k gets only a random fraction ξ k−1 n of the momentum loss of the mother eddy at level k − 1, and other part is dissipated. Thus, the momentum balance for the eddies at scale k is given by where m k and m k−1 are the masses of eddies at scales k and k − 1, respectively, λ k = e −kγτ is a constant damping at level k. Clearly, here we approach the breaking factor linearly in k. In addition, we considered ξ (k−1) n as a uniform random distribution with values between 0 and 1. This parameter can be understood as a random fraction of transferred momentum from the scale k − 1 to the next scale k. In other words, the term m k−1 (1 − λ k−1 )v (k−1) n represent the lost of momentum by friction in the k − 1 scale that serves as a driving force for the next scale, k. Moreover, as the friction at the scale k is produced by several smaller eddies, only a fraction of energy is transferred to one eddy in the next scale. It is important to mention that other choices for ξ (k−1) n may be considered. However, as argued by Beck (1994), and numerically shown by Gallo-Méndez & Moya (2022), those details do not significantly affect the properties of the VDFs in the stationary state, which is the focus of our study.

The equations are given by
Here it must be noted that the factor c k represent the inverse of β k = m k /m k−1 , a characteristic parameter of the β-model (Benzi et al. 1984), which serves as a coupling factor between spatial scales, through the ratio between the masses of the eddies at the scales k − 1 and k, respectively.
Here, such relation comes from equation (6). In addition, the value of c k is also related to the width of the PDF at the scale k, which is related with the thermal speed of the system in the case of velocity fluctuations. As mentioned by Beck (1994), the relation between the c k values and formation and properties of heavy tails in the PDF is weak. Therefore, for the sake of simplicity and to focus on the description of the tails of the PDFs, it is customary to select c k = c, constant for all scales. Hence, in order to compute in a simulation, we consider c = 2 for all scales. This means that, on average, a mother eddy preferentially interacts with eddies what are half its mass. Moreover, v (k) n can be interpreted as velocity of the center of mass of an eddy in the k-scale at time nτ . Finally, we notice that the system is determined by two characteristic times, which are condensed in the parameter γτ .
This factor is actually the ratio between these characteristic time scales. First we have τ , the time scale of fluctuations in the large scale ∼ ℓ 0 , given by τ ∼ ℓ 0 /w, where w is the thermal velocity of the fluctuations. On the other hand, in the Langevin Equation (2), γ is the dynamic viscosity given by γ ∼ ν 0 /ℓ 2 0 , with ν 0 the kinematic viscosity. Then, we have γτ ∼ ν 0 /wℓ 0 ∼ R −1 λ . Therefore, γτ factor is inversely proportional to the Reynolds number of the fluid, R λ (Beck 1994(Beck , 1996. Then, as γτ is the only free parameter of the model, in the CML, the dynamics of the system is determined by the level of turbulence of the media, measured by R λ . The CML equations (8) represent a straightforward model to study the consequences of turbulence as a one dimensional system at different scales.

ANALYSIS AND RESULTS
To numerically solve the CML equations (8) we consider an ensemble of N = 10 6 eddies in each of the k scale, from k = 1 up to k max = 50, taking into account that the coupling factor between scales is constant through all scales (c k = 2). In addition, to test our model with different levels of turbulence (measured by the Reynolds number) we solve the equations using five different values for γτ , between 10 −4 and 10 −2 , or, similarly, using different values of Reynolds number, between 10 2 and 10 4 .
As initial condition we select a delta function distribution for velocities centered at v (k) 0 = 0, and evolve the CML equations in time until the system reaches a steady state with constant energy E.
Namely, (1), where A κ,δ , w, κ and δ are fit parameters. We had special attention to the relation between these quantities respect to the scale k for different values of R λ .
Skew-Kappa distributions correspond to a generalization of the Maxwellian distribution mainly in two ways. The first is that it is widely used to model out-of-equilibrium systems in which the VDF exhibits large kurtosis due to powerlaw supra-thermal tails for higher energies (or velocity). Under such context, Eq. (1), depending on the value of the κ parameter the VDF will model the distribution with a quasi-thermal core (for v ≤ √ κw), and power-law tails for larger velocities. The extent and energy of the tails increases with decreasing κ, and the VDF collapses to a Maxwellian in the limit κ → ∞. And the second, Skew-Kappa distributions, as its name says, have a skewness parameter, δ, which is related to the turbulent nature of the system (Beck 2000;Zenteno-Quinteros et al. 2021).
Furthermore, it should be noted that this fit is only valid for small values of skewness parameter, i.e.
δ ≪ 1, and in a limited range of v such that v/w ≤ |2/δ + δ/8|.  (bottom), depending on the spatial scale k, for each value of Reynolds number, R λ = 10 2.0 , 10 2.5 , 10 3.0 and 10 3.5 . For the first two variables a power law function with exponents −α and β is overploted to the data data, respectively. The last variable, δ, is compared with its average value in each case.
In order to characterize this behaviour, we performed a systematic analysis of the relation between δ, κ, R λ and the scale k. Top panels of Figure 2 show the adjusted κ as a function of k for different values of R λ = (γτ ) −1 (from left to right, R λ = 10 2.0 , 10 2.5 , 10 3.0 and 10 3.5 , respectively). It should be noted that with this selection of Reynolds number values we were simulating the behavior of a plasma system in the context of space plasmas, as the solar wind, the plasma sheet, among others.
Due to numerical issues, we could not evaluate the dynamics of Reynolds numbers above 10 4 . This limitation was due to the numerical resolution required to accurately simulate the behavior of the plasma system. However, regarding the range of Reynolds numbers that we were able to consider, we can find that the effective Reynolds number for the solar wind was determined to be 260.000 ± 20.000 at 1 AU from the Sun as shown by Weygand et al. (2007) In addition, regarding κ parameter, we observe a clear decaying power law behavior κ ∼ k −α .
For each value of R λ the best power law is also included as a superposed black line in each panel.
Remarkably, for all considered values of Reynolds number, α ≈ 1.74. Furthermore, from the figure it is also possible to notice how this power law behaviour also depends on the value of R λ , as, in general, κ increases with increasing the Reynolds number. However, from Figure 2 we can recognize two regimes where the purposed power law model is no longer representative: first, for low values of k (or high spatial scales), given that the points are expected to be dispersed since the system tends to the Maxwellian state. As Maxwellian is the limit when κ → ∞, when the system is more Maxwellian than Kappa distributed, the algorithm found high values for κ (larger than 10 2 ), and therefore there is low precision in the convergence of the fit method. Second, for high values of k (small spatial scales), where the points saturate to a fixed value of 3/2. Since the κ-like fit we use has a singularity at this value, where the tails of the PDF seem too large to be described using a Kappa distributions.
In such cases, other models may provide a better description of the PDFs.
It is important to mention that the validity range in k also depends on the Reynolds number, as shown on the top of the panel. In short, the slope seem to be constant for all values of R λ , but we can see how the curve moves to the right for increasing Reynolds number, and also the power law will be presented as long as 3/2 < κ ≲ 10 2 . For all other values of κ, the trend is broken. We attribute this behavior of the Kappa parameter to the correlation gain as we advance in the k index.
Let us recall that scale k contains part of the information of all k − 1 previous ones. In this sense, we can say that smaller scales are more correlated. Therefore, the higher the level of correlation, the smaller the value of κ. This is consistent with the arguments given by Livadiotis & McComas (2021). There, the authors discuss the relationship between κ and the level of correlation of a system.
They justify a non-additive entropy since a correlation term in the total entropy which depends on the entropy of two subsystems and κ. In order to include more information about our system, we defined δ κ = δ/(κ − 3/2), which is indeed a combination of two fit parameters, such as δ and κ. This definition is not a coincidence since δ κ is naturally generated from Eq. (1). In this way, we can interpret δ κ as an effective skewness parameter. As we can notice, on the bottom of Figure 2, we plot δ as a function of the scale level, k.
As expected, due to the arguments exposed by Beck (2000), we should not see a dependence between skewness parameter and the scale, but only on the Reynolds number, namely, δ ∼ R −1/2 λ . In fact, we can see, from the plot of δ parameter, a non-strong-dependence on the scale. Therefore, it is expected that, if κ decays as a power law, for large k, then δ κ will be a power law with a positive slope, just as we see in the Figure mentioned.
As the shape of the data suggest, we fit a power law for the effective skewness, namely, δ κ ∼ k β .
Due to δ κ decays as the inverse of κ, and δ is apparently constant in a certain range, we should expect α ≈ β in that range of k. Although with some error this is fulfilled. As we can see in the figure, δ κ is well fitted by a power law in the same range as κ, but there is a small discrepancy between the mean value of α and β, given by ⟨β⟩ − ⟨α⟩ = 0.29. Nevertheless, this difference is easily attributable to the sensitivity of the δ parameter when we measure it. Indeed, we can see the above in more detail in Figure 3a and 3b, where we plot the values of slopes α and β for each value of R λ , and we compare them with our reference value 5/3, which is apparently a constant value as studied in Gallo-Méndez & Moya (2022). Finally, in Figure 3c, we plot the values of skewness parameter, δ, for each value of Reynolds number, R λ . In order to compare these values with the theory, we overplot the curve R −1/2 λ . We notice that there is a correspondence between the data and the theoretical curve, except for the last point at R λ = 10 4.0 , which is likely to have primarily numerical error. Since in that case the γτ number in the equations (8) is quite small, it is expected to generate numerical issues in the convergence of the system.
Furthermore, considering the relation between the scale k and the κ-parameter, our numerical results indicate that is the relation between κ and the spatial scale. Therefore, the equation (10) is a functional relation which is valid for a certain range of spatial scales, according to the CML. And not only that, since as we see on the bottom of the Figure 3, we conclude that, indeed, the skewness parameter relates with the level of turbulence, namely, This last result could be a very useful tool at the time to measure the level of turbulence in a certain system. As the same way in Beck (2000); Rizzo et al. (2004), when we have two measuring tools (i = {1, 2}), at two different points, separated by a characteristic length ℓ, which measure a variable u i , we could construct the difference v = u 1 − u 2 (as in Eq. (3)), framed within a PVI quantity.
In this way we can obtain the PDF of PVI, such as the authors show in their works. There, it is clear that when the distribution is built, it has a slight asymmetry, which can well fitted by a Skew-Kappa distributions. The latter allows an indirect measurement of the level of turbulence through the quantification of the skewness parameter δ.  (2022)). On the other hand we observe that the skewness, quantified by δ parameter in a Skew-Kappa distribution, is also strongly related with the level of turbulence but not with the scale k.
The numerical results of our study show the relation between the spatial scale, the Reynolds number, and the κ and δ parameters. We observe how κ exhibits a power law behavior given by Eq. (10). In that regard, we can conclude that the nature of the force used in the model does not change the shape of the dependence between κ and the spatial scale. In this way, when the Reynolds number grows, the κ parameter also increases. This fact is expected since in a system with large enough Reynolds number, strong turbulence should inhibit the internal correlations leading to non Maxwellian PDFs.
Secondly, we compare the values of skewness, parameterized by δ, with respect to the spatial scale k. It is possible to notice that the asymmetry is weakly dependent on the scale, in the same way as Beck (2000) suggest and is concluded in Gallo-Méndez & Moya (2022). This is because the theory suppose that δ does not depend on the κ parameter i.e. not depends on the scale either. Therefore we are able to found a match between the simulation approach with theoretical works.
While our research primarily focuses on other aspects of the system dynamics, it is also worth mentioning that our model exhibits the characteristics of the inertial range. Specifically, our model reaches equilibrium at larger spatial scales initially and gradually achieves a stationary behavior at progressively smaller scales, with clear scaling laws between δ κ , κ, and k. However, as all considered scales follow the same equations, the model is not able to exhibit a spectral brake between the inertial and kinetic ranges as observed in the solar wind. To do so we should have to introduce new rules to Equations (7) and (8). Moreover, other characteristic of turbulent flows (as intermittency) are also not considered. We acknowledge that this may happen but this it is not the primary focus of our work, and we have not thoroughly explored these phenomena in our article. In summary all these issues are relevant and worth to study, but investigating the potential intermittent behavior of plasma fluctuations, the properties of the inertial range, and the possible inclusion of a spectral break, are all beyond the scope of our current study.
That being said, as far as κ parameter is concerned, the results are qualitatively consistent with the observational study by Pollock et al. (2018) for the other definition of PVI, which is construct with a characteristic length scale, the results seem to be quite similar to those reported by Beck (2000), who analyzed data from a jet experiment of turbulent flows.
To complement those studies, we propose the systematization of the quantification of skewness parameter in PDFs of spatial PVI, evaluating the fit for several fluids with different Reynolds numbers.
Not only that, but also we suggest that the methodology exposed in this article could be included in future observational works in order to extend the study of PVI in space plasma. And though this has been explored in studies as Chasapis et al. (2018), we emphasize the importance of making the evaluation on the PDF without absolute value by making the difference between both measurements, and thus obtain an asymmetric distribution. In summary, on the basis of a simple numerical model of a turbulent system we have found that the PDFs of velocity differences fluctuations follow Skew-Kappa distributions, and that there is a robust relation between the scaling of κ, the skewness parameter,δ, the spatial scale level, k, and the Reynolds number of the flow. We expect these results to be useful to characterize turbulent systems in different contexts, specially in space plasma environments, where the Reynolds number is not always easy to obtain (Parashar et al. 2019), and our numerical predictions to be tested by spacecraft observations or in experimental setups.