Influence of MHD Turbulence on Ion Kappa Distributions in the Earth's Plasma Sheet as a Function of Plasma β Parameter

The possible influence of MHD turbulence on the energy distributions of ions in the Earth's plasma sheet was studied using data taken by the THEMIS satellites. Turbulence levels were traced using eddy diffusion coefficients (D), of which we measured one for each Geocentric Solar Magnetospheric (GSM) coordinates every 12 min. Ion fluxes between 1.75 and 210.5 keV during the same time windows that correspond to mainly suprathermal populations were fitted to Kappa distribution functions, which approximate a Maxwellian distribution when the κ-index (κ) is large. We found that the distribution of the eddy diffusion coefficients is bimodal, independently of both the eddy diffusion component and the plasma beta (β) parameter, which is defined as the ratio between plasma and magnetic pressures. The main peak corresponds to turbulent plasma flows with D > 103 km2 s−1. In such cases, the impact of turbulence on the κ index depends on the value of β and also on the direction of the turbulent transport. For eddy diffusion perpendicular to the neutral sheet, the values of κ decrease as Dzz increases for β < 2; while for higher values of β, κ increases with Dzz. For the other two directions, the values of κ decrease as D increases. This last tendency is stronger for β ~ 1 but almost null for β ~ 10. The secondary peak in the distribution of D values might represent quasi-laminar flows forming part of very large vortices, correct detection and description of which is beyond the scope of this study.


INTRODUCTION
Since the beginnings of the Space Age, it is known that various regions of the Earth's magnetosphere are filled with turbulent plasmas (see, for example, Sonett, 1960;Ness et al., 1961). Such dynamics are somewhat expected if the Earth's magnetic field is considered as an obstacle in the path of the incoming turbulent, supersonic, and superalfvenic solar wind, which leads to turbulent wake formation. Nonetheless, the magnetotail presents a more complex structure than an ordinary wake, as the plasma density and beta (β) parameter in the plasma sheet are much higher than in the tail lobes, and the plasma sheet is much more turbulent (e.g. Antonova, 1985;Montgomery, 1987;Angelopoulos et al., 1993;Borovsky et al., 1997). This behavior constitutes one of the many unsolved problems in magnetospheric physics (Borovsky et al., 2020).
Detailed studies of the properties of plasma sheet turbulent flows in the magnetotail started at the end of the 90's. For instance, Borovsky et al. (1997Borovsky et al. ( , 1998 showed that unpredictable random fluctuations of flow velocity are observed in the plasma sheet and that their amplitudes are much larger than the mean flow velocities. They also obtained the probability distributions and autocorrelation times of flow velocity and concluded that the observed properties correspond to magnetohydrodynamic (MHD) turbulence. The flow velocities were normally shown to be sub-Alfvenic, which is especially important for frozen-in condition applicability and the observed properties of turbulent fluctuations. Borovsky and Funsten (2003) later confirmed that the observed turbulence is essentially formed by eddies (eddy turbulence) and determined that the typical scale of such eddies is 4,000-10,000 km, or even greater, with a mean value of about 1.6 Earth Radii (R E ). Later, though Weygand et al. (2005) obtained somewhat smaller estimates of 6, 000 km for the eddy sizes. Successful direct observations of eddy vortices in the plasma sheet are only possible if measurements are carried out fairly simultaneously by two or more satellites separated by long distances, such as THEMIS satellite mission (Time History of Events and Macroscale Interactions during Substorms; Angelopoulos, 2008). Plasma vortices of such scales were observed by Keika et al. (2009), Keiling et al. (2009), Panov et al. (2010 using THEMIS data. The eddy diffusion coefficients are a useful way to quantify MHD turbulence. Borovsky et al. (1998) proposed a method for calculating this coefficient using data from a single satellite as a function of the root mean square velocity and autocorrelation time.
Many authors have carried out measurements of eddy diffusion coefficients in the plasma sheet using data from different satellites. In general, most results point to D αα ∼ 10 5 km 2 /s (e.g. Ovchinnikov et al., 2000;Troshichev et al., 2002). Subsequent studies of turbulent eddy diffusion in the plasma sheet have shown that eddy-diffusion coefficients increase with Earth distance toward the tail (Stepanova et al., 2009Pinto et al., 2011), but are nearly constant at mid distances and slightly decrease in the distant tail (Troshichev et al., 2002). The eddy diffusion was also found to be sensitive to the geomagnetic activity, increasing during geomagnetic substorms, especially during the expansion phase process (Stepanova et al., 2009. Studies of the stability of the plasma sheet commonly consider laminar flows only (see Petrukovich et al., 2015;Sitnov et al., 2019). Nonetheless, the presence of turbulence requires to modify this approach. In this case, the stability of a turbulent plasma sheet can be reached due to a balance between regular and turbulent transports: Large-scale dawndusk electric field allows plasma advects toward the tail center; meanwhile, turbulent transport could compensate for this advection, thereby forming a stable plasma sheet (Antonova and Ovchinnikov, 1996, 1999, 2001Antonova, 2002;Stepanova and Antonova, 2011). The same approach can be applied to explain the dawn-dusk asymmetry in the thickness of the turbulent low latitude boundary layer: a balance between the regular transport and transport due to turbulent eddy diffusion is assumed in combination with a dependence of this thickness on the orientation of the interplanetary magnetic field (IMF) (Antonova, 2005).
However, strictly speaking, the MHD approach for describing magnetospheric plasmas is valid only when the particle distribution functions are Maxwellian. Meanwhile, this may not be the case, as the energy distributions of particles in the magnetosphere, especially in the plasma sheet, are generally not Maxwellian. In fact, enhanced high energy tails are often observed in the energy distributions (e.g. Christon et al., 1989;Espinoza et al., 2018).
Kappa distribution functions (Livadiotis, 2017) have proven to be useful representations of the differential particle fluxes in the magnetosphere (e.g. Olbert, 1968;Christon et al., 1989;Haaland et al., 2010;Navarro et al., 2015;Runov et al., 2015;Viñas et al., 2015;Yoon and Livadiotis, 2017;Kirpichev and Antonova, 2020). In general, these functions will recover more information than a Maxwellian function, albeit the addition of one more parameter: the kappa index (κ). While high kappa values imply a higher slope of the distribution at the highest energies (hence a more Maxwellian behavior); low kappa values indicate the presence of a substantial population of energetic particles. Their presence can partially be understood by the fact that particle collisions are rare in these plasmas, and there seems not to be an efficient mechanism to thermalize them (evolve the distribution toward a Maxwellian function) over timescales shorter than the main dynamical scales that govern the systems. In this regard, although the derivation of the Kappa distribution from first principles is under strong debate, the Kappa function might be understood as a generalized case of the Maxwellian distribution (Olbert, 1968;Vasyliunas, 1968;Tsallis, 1988Tsallis, , 2004Tsallis, , 2009Collier, 1993Collier, , 1999Nauenberg, 2003Nauenberg, , 2004Livadiotis and McComas, 2013;Livadiotis, 2015;Treumann and Baumjohann, 2020).
The general form of the Kappa distribution is denoted by f : where f is the phase space density, E is the kinetic energy, n is the particle density, m is the particle mass, Ŵ is the Euler gamma function, and κ and E c are the κ-parameter and characteristic or core energy, respectively. The κ index is a measure of the energy spectrum slope of the suprathermal particles forming the tail of the velocity distribution function. Indeed, for κ → ∞, Equation (1) becomes identical to the Maxwellian distribution. The Kappa distributions given by Equation (1) exhibit a thermal core with characteristic energy E c and suprathermal tails, such that the total characteristic particle kinetic energy E total is given by which enables a straightforward comparison between Kappa and Maxwellian distributions, and to outline the effects of suprathermals as shown by Lazar et al. (2015Lazar et al. ( , 2016. Over the last decades, it has become clear that Kappa Distributions can successfully describe ion distribution functions in the plasma sheet. For instance, Christon et al. (1989Christon et al. ( , 1991 used the particle instruments onboard the International Sun-Earth Explorer 1 (ISEE 1) and found that kappa-index ranges between 4 and 8 for both ions and electrons, with a most probable value between 5 and 6. Later, Haaland et al. (2010) found that the κ-index ranges between 3 and 6, using data from the Cluster satellites. Espinoza et al. (2018) also used the Kappa approximation to model ions and electrons flux spectra along the plasma sheet. Their statistical results reveal that the electrons have smaller kappa κ e than ion kappa κ i , which suggests electrons non-thermal properties are stronger than ions. Besides, their results in the relative numbers of energetic ions show a significant dawn-dusk asymmetry, being lower in the dusk-side, which increases during substorms. This is consistent with the previous study of Wing et al. (2005).
Stepanova and Antonova (2015) utilized Kappa distributions to fit ion and electron flux spectra for five events in which the THEMIS satellites were oriented in the plasma sheet. They obtained snapshots of kappa index properties that show a tendency for the electron κ-index to be high toward tailward direction; while the situation with the ion κ-index is less clear. In two of the five events analyzed, the kappa-index decreased toward the tail. To explain this, it was proposed that the stochastic acceleration responsible for the diffusion in energy space could lead to hardening of particle spectra (low kappa values) near the Earth region. Simultaneously, the diffusion in velocity space could lead to softening of particle spectra (Collier, 1999). Considering that for some events, the relaxation took place further away from the Earth, it was proposed that plasma transport from the Earth to the tail takes enough time for the relaxation of spectra (aging).
Nonetheless, it is well-known that the dawn-dusk electric field drives a regular transport toward the Earth. On the other hand, the bursty bulk flows mainly produce much of the particle transport earthward, which is faster than regular convection. The dipolarization fronts also produce plasma transport in the earthward direction (Runov et al., 2012). At the same time, the form of ion spectra fitted by Kappa distribution is conserved during the dipolarization, and on average, no non-adiabatic acceleration of ions in dipolarization flux bundles were present (Runov et al., 2015). Consequently, none of these sources of plasma transport could explain the gradual increase of kappa toward the tail due to tailward plasma transport accompanied by the "aging" of the distribution function.
Stepanova and Antonova (2015) proposed that particle transport from the inner magnetosphere toward the tail could be attributed to turbulent eddy diffusion. In this case, the time available for Kappa distribution function relaxation to Maxwellian will depend on balance between the processes leading to both transport directions. The characteristic time of turbulent transport between 3 and 12 h were estimated. At the same time, they calculated the average bulk velocity and found that most velocities were directed toward the Earth. However, there are points at which the bulk velocity was directed tailward, thus reducing the time available for relaxation. Whereas, the low kappa values were found precisely for these events, suggesting that Kappa distribution evolution is mainly due to relaxation in the velocity space. It could be that MHD turbulence is what modifies the particle distribution functions. To date, there are no specific studies about a possible relation between the kappa index and MHD turbulence-considering turbulent transport, local turbulent acceleration, or conversely turbulent mixing; that might lead to the Maxwellization of the distribution functions. In this study, we explore the relation between the ion kappa index and the eddy diffusion coefficients. This result contributes to a better understanding of the possible influence that MHD turbulence may have on the departure of the ion distribution functions from Maxwellian.
The organization of this paper is as follows: In section 2, we describe the instruments and methodology used to obtain the ion Kappa distribution parameters; In section 3, we present the results of the analyses and explore the relationship between kappa and three components of the eddy diffusion tensor. In section 4, we discuss the results; and in section 5, we summarize and conclude our findings.

INSTRUMENTATION, DATA SELECTION, AND ANALYSIS
We utilized data taken by the THEMIS mission (Time History of Events and Macroscale Interactions during Substorms; Angelopoulos, 2008) during the years 2008-2009, downloaded via the THEMIS ftp website. 1 Quasi-static magnetic field measurements were taken using the Flux Gate Magnetometer (FGM; Auster et al., 2008). Ion spectra were obtained by combining the measurements of the Electrostatic Analyzers (ESA; McFadden et al., 2008), which operates at lower energies, from a few eV to 25 keV, and the Solid State Telescopes (SST; Angelopoulos, 2008), which is sensitive to energies above 30 keV. We used level 2 full particle energy flux spectrogram. The angular distributions are not considered in the measured data since the pitch angle dependency has been averaged. Taking into consideration that ions in the plasma sheet are typically isotropic, the loss of anisotropy information is not critical. Moreover, due to the average process the number of counts in each energy channel increases, which is beneficial for the study of energy distributions. A mass-spectrometer was not included in the THEMIS instruments and it is impossible to differentiate ions of different species; therefore, we refer to them as just ions.
We restricted our analyses to ions energy range of 1.75-210.5 keV, to ensure that the actual fits of the tail and core parts of ion distribution function are obtained following the method used by Espinoza et al. (2018). We discarded lower energies due to contamination from the photoelectrons and spacecraft potential. Similarly, we discarded higher energies due to contamination by solar cosmic rays, energetic electrons, and a low number of counts. This energy range does not contain relativistic ions, hence we use the ordinary κ−distribution function, without the relativistic corrections suggested by Scherer et al. (2018Scherer et al. ( , 2019, Lazar et al. (2020).
Example of ion energy flux spectra measured by combining both particle instruments (ESA and SST) onboard THEMIS  Figure 1. The circles on the plots are the average of the spectra obtained for the 12 min time windows. The open circles represent measurements from the ESA, while the filled circles represent the SST. For the average flux data, the error bars represent the spread between the maximum and minimum values observed. The error bars were found to vary significantly between the ESA and SST data, so they were normalized in the same way as Espinoza et al. (2018). The inverse squared of the error bars is used to define weights for the fits, which were performed over the averaged data using a non-linear least-squares method combined with the Levenberg-Mardquart algorithm. We inspected hundreds of spectra visually and decided to work only with the fits that give a reduced chi-squared χ 2 < 100.

satellites (solid lines) is illustrated in
The observed differential energy spectra of ions were fitted to the model obtained by transforming the Kappa distribution from Equation (1) to differential energy fluxes: This expression agrees with the expression of the differential particle flux (I(E)) used, for example, by Vasyliunas (1968), Christon et al. (1989), Olsson and Janhunen (1998). As in Stepanova and Antonova (2015), Espinoza et al. (2018), here we use the differential energy flux, that is F(E) = E·I(E) (Lyons et al., 1985;Baumjohann and Treumann, 1997). The differential energy flux was used since this is how the information related to particle distribution functions is given by the THEMIS team. Following Borovsky et al. (1998, see their appendix), to estimate the eddy-diffusion coefficient D zz , the eddy transport is assumed to be a Markov process (i.e., each plasma displacement z is independent of previous displacements). This is valid when displacements are separated in time by the autocorrelation time τ auto of the flow velocity. In the Markovian picture, the diffusion coefficient is given as; A typical z displacement which has a z-component of fluid velocity is z = V z τ auto . With this, the turbulent transport is evaluated by determining the eddy diffusion coefficient tensor denoted as D αα : where τ αα is the decay time of the autocorrelation function of the ion bulk velocity component V α ; which acts as a measure of the persistence of a fluctuating bulk velocity, and V rms is the root mean square (rms) of the α component of the velocity fluctuations around a mean value, which can be determined from Both the diffusion coefficients and the parameters of the kappa fits were obtained from 12-min intervals shifted every 6 min. This allowed us to use about 240 bulk velocity measurements in each interval to calculate V rms and τ αα , as described in Stepanova et al. (2011). All measurements were constrained to the following Geocentric Solar Magnetospheric (GSM) coordinate system: To ensure that the measurements are performed in plasma sheet-like plasmas, the time intervals were selected based on the following criteria: ion number density (n i ), ion temperature (T i ), x components of the magnetic field (B x ), bulk velocity (V x ) were restricted to n i ≥ 0.1 cm −3 , T i ≥ 1 keV, |B x | < 20 nT, and |V x | < 300 km/s. As well, we considered only cases in which β > 0.1; where β is the ion plasma parameter defined as the ratio of the ion plasma pressure to magnetic pressure: β = ( P i B 2 /2 µ 0 ), where B is the total magnetic field strength and µ 0 is the magnetic permeability, and P i is the plasma pressure. The values of P i and ion bulk velocity were taken directly from the level 2 data products containing the ESA macro-parameters (see the SPEDAS THEMIS documentation 2 ). The contribution of ions with energies above 25 keV (detected by the SST, Angelopoulos, 2008) to the bulk velocity has been evaluated and found to be negligible (Lee and Angelopoulos, 2014), hence the influence of SST on macro parameters, such as P i , is not strong. Thus, we utilized only the ESA measurements to calculate beta parameter and rms velocities. Figure 2 illustrates the data analysis described above, using data obtained from THEMIS C satellite on February 26, 2008, between 05 to 08 UT during quiet geomagnetic conditions. In this case, the satellite was inside the plasma sheet until ∼7:45 UT. Both the eddy diffusion coefficients and the κ index were obtained while |B x | < −20 nT. As mentioned at the beginning of this section, we used 12-min time intervals to calculate the eddy diffusion and kappa parameters. Nonetheless, we acknowledge that this choice for the duration of the interval might interfere with the calculation of the eddy diffusion coefficients. To estimate this effect we present two hodograms in Figure 3, which show the evolution of particles velocities in the (V x ,V y ) and (V x , V z ) planes. We calculated both the average and the root mean square velocities (V rms ) for a 3-h long interval, as well as for selected 12 2 http://themis.igpp.ucla.edu/software_docs.shtml min intervals. As it can be seen, the average velocities calculated over 12-min intervals are of the same order as the average value over 3-h. This suggests that the 12-min intervals are enough to recover eddy diffusion coefficients. Figure 4 illustrates the spatial coverage of the selected intervals projected in the X-Ygsm and X-Zgsm planes in the geocentric solar magnetospheric system. Each interval is represented by a dot, which is color-coded with the measured plasma β parameter. Furthermore, the intervals selected correspond to the quiet time geomagnetic conditions, which were identified using the criteria applied by Stepanova et al. (2011). They are based on the variation of the high resolution (1-min) auroral electrojet lower index (AL) obtained from the OMNI database. This was achieved by subjecting the 1-min resolution of the AL index to AL ≥ −100 nT and the slope absolute value |s| to ≤ 1/2 nT/min, which indicate how fast the AL-index changed during a minute (more details about this analysis are given in Stepanova et al., 2011). Figure 2, for relatively low β the values of κ in some intervals anti-correlate with the eddy diffusion coefficients. However, whether this behavior persists for all orbits and for all values of β is unclear. Thus we made an extensive statistical study of a set of hundreds of thousands quiet time 12-min-long intervals, for which we have calculated D, κ, and β. Considering that both D and β cover a few orders of magnitude, we define a grid in the logarithmic (β, D) space using a cell size of log 10 β = 0.2, and log 10 D = 0.5 within the range −1 < log 10 β < 2 and 0 < log 10 D < 7. This was used to create the color-coded plots shown in Figure 5. Figure 5 shows the number of measurements (N) and the average values of κ i , in each bin. The empty bins contain less than ten measurements. As it can be seen in the left column [panels (A), (C), and (E)], for all eddy diffusion coefficients components, and for all values of β, the distributions of eddy diffusion coefficients have two maxima. The main peak corresponds to turbulent plasma, with D between 10 4 and 10 6 km 2 /s. These high values are the main contribution to the values of the eddy diffusion coefficients averaged within each β bin, which are shown with a white solid line in the figure. The secondary maximum might correspond to nearly laminar flows (D ∼ 10 2 km 2 /s), which could be part of vortices larger than the maximum vortex size that we are able to detect (determined by the 12-min intervals chosen for this study). In such case, the rms velocity would be low, and the autocorrelation time would exceed 12 min, hence could not be measured correctly. In other words, this method may be insensitive to very large vortices, especially considering that the satellite velocity is often of the same order of magnitude as the plasma velocity averaged over turbulent eddies. Therefore, we restrict the following analysis of a possible relation between κ and D αα to turbulent flows only.

As seen in
In order to understand whether there is a relation between κ and D αα for a fixed value of β (as suggested by the right panels of Figure 5), we use the linear function κ = A log 10 D + B to fit κ    and log 10 D. Figure 6 shows these dependencies for some selected values of β. The fitted data, corresponding to turbulent flows, are plotted with filled circles (D > 10 3 km 2 /s). The vertical error bars represent the standard deviation of the κ values within each bin; meanwhile, the horizontal error bars reflect the D bin's width.
Despite the strong dispersion of the κ values in each bin, there are systematic trends in the behavior of the measured slopes (A) and intercept (B), as β increases (Figure 7). These trends are different for each D αα component. For the x-component the κ indices almost do not depend on D: A xx ∼ 0 except for β ∼ 1, for which A xx is minimum and negative; thereby implying that the κ-index decreases with D xx . The trend for Y is similar to the one observed for X but clearer. The behavior of A zz is different: it monotonically increases with β and changes its sign, implying FIGURE 6 | Dependence of ion κ indices on mean log 10 D for different plasma β. The first column is for D xx , while second and third columns are for D yy and D zz , respectively. The function κ = A log 10 D + B was fitted to the points plotted with filled circles, and the results are plotted with a solid line.
that while κ decreases with D for low β, κ increases with D for β ≥ 2. The latter suggests that eddy diffusion in plasmas with high β might make the distribution functions more Maxwellian.

DISCUSSION
Our study shows that MHD turbulence in the plasma sheet might have an intermittent character. As seen in Figure 5, for most β values the eddy diffusion coefficients vary over a wide range (10 1 -10 6 km 2 /s) and present two peaks. The majority of the coefficients concentrate around 10 4 -10 5 km 2 /s, which correspond to medium scale vortices with scales of ∼ 10, 000 km that contribute to turbulent transport.
The eddy diffusion coefficients lower than 10 3 km 2 /s may correspond to very large vortices for which the coefficients are underestimated due to the 6-min limit in the determination of the autocorrelation time (half of the 12-min-long windows used for the analyses). The duration of the time interval used can also affect the values of the average V rms over the interval. Borovsky et al. (1997) used a 2-h interval and obtained average bulk velocity values of the order of 10 km/s. In our study, these averages range between a few km/s to a few hundreds km/s, while the typical values are a few tens km/s. Inside big vortices the average V rms will be low, because within the vortex the velocity of the vortex is subtracted and only small fluctuations remain.
On the other hand, the presence of stable quasi-laminar flows would agree with the results of Angelopoulos et al. (1999), who analyzed the properties of velocity fluctuations and concluded that the geomagnetic tail is a system that has properties of intermittent turbulence and exhibits sporadic variability. This fact is also reflected in multi-scale features of magnetic fluctuations in the near-Earth magnetosphere (Lui, 2001;Consolini et al., 2005). The intermittent behavior can be observed over all the plasma sheet for example see Figure 2, also Stepanova et al., 2011). However, the average values of the eddy diffusion coefficients are 10 4 -10 5 km 2 /s (solid white lines in Figure 5), indicating a strong presence of turbulence, and suggesting that turbulent transport may play an important role in magnetospheric dynamics.
The evolution observed of κ toward the geomagnetic tail has been associated with the transport of particles induced by turbulent eddy diffusion. For slow plasma transport from the Earth to the tail, there is enough time for the Kappa distribution to relax toward Maxwellian distribution due to diffusion in the velocity space (aging). This mechanism of relaxation (proposed by Collier, 1999) explains also the increase of the κ-index with the E c that has been observed (Kirpichev and Antonova, 2020;Eyelade et al., 2021).
Our results show that low values of κ are observed for the highest eddy diffusion coefficients D xx and D yy . In this case, turbulent transport in the X and Y directions is faster, and the available time for Kappa distribution functions to relax to a Maxwellian is shorter. Thus the Kappa distributions cannot thermalize, which is reflected in the low κ indices.
However, some energetic tails might also be the consequence of reconnection outflows, as observed in many experiments and also modeled by El-Alaoui et al. (2010. These outflows appear as part of turbulent cascades during MHD simulations, when the computer codes combine both low resistivity and small grid spacing, thereby obtaining comparatively large values of the magnetic Reynolds number. El-Alaoui et al. (2010, argued that the formation of localized reconnection regions is the main process driving turbulence in the plasma sheet. These processes would also lead to a decrease in κ-index. On the other hand, for the Z direction the eddy diffusion seems to increase with the κ-index, for high values of β. This could be related to turbulent mixing, which is commonly observed in the plasma sheet. The mixing length of the average plasma sheet is ∼10,000 km (Borovsky et al., 1998), leading to equalization of temperature across the plasma sheet, which is commonly observed (Huang and Frank, 1994). The plasma sheet extends to thousands of Earth's radii toward the tail in the X direction, while in the Y direction to tens of radii and just about ten Earth's radii in the direction perpendicular to the neutral sheet Z. This might explain why our results show that turbulent plasma mixing might prevail in the Z direction; meanwhile the fast turbulent transport and local reconnection is relevant in the X and Y directions. The efficiency of all the aforementioned processes could depend on the plasma β parameter, as reflected in our study. Nonetheless, more exhaustive statistical studies are necessary to untangle these effects.

CONCLUSIONS
Observations of ion Kappa distribution made by the multisatellite THEMIS mission in the magnetotail plasma sheet were statistically explored in conjunction with eddy diffusion coefficients D and the plasma β parameter. Our study reveals the presence of turbulent flows (D ∼ 10 4 − 10 6 km 2 /s) alternated with quasi-laminar flows (D ∼ 10 2 km 2 /s), which might belong to large vortexes that are beyond the detection limits of our method. For turbulent plasmas, several processes related to MHD turbulence lead to either an increase or decrease of the κ index, depending on the value of β and the direction of the turbulent transport with respect to the plasma sheet.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All the research analysis in this study was carried out by AE, CE, and MS, while IK assisted in analyzing Kappa distribution functions. IO was engaged in the analysis of eddy diffusion transport, while EA was involved in the theoretical interpretation of the results. All authors contributed to the manuscript revisions and approved the final version of the manuscript.

FUNDING
This work was supported by Agencia Nacional de Investigación y Desarrollo de Chile (ANID) grants 21181777 and 1191351. MS acknowledges support from University of Santiago Chile through DICYT grant number 042031S. CE and MS acknowledge support from AFOSR (FA9550-19-1-0384).