Statistical analysis of stochastic magnetic fluctuations in space plasma based on the MMS mission

Based on the Magnetospheric Multiscale (MMS) mission we look at magnetic field fluctuations in the Earth's magnetosheath. We apply the statistical analysis using a Fokker-Planck equation to investigate processes responsible for stochastic fluctuations in space plasmas. As already known, turbulence in the inertial range of hydromagnetic scales exhibits Markovian features. We have extended the statistical approach to much smaller scales in space, where kinetic theory should be applied. Here we study in detail and compare the characteristics of magnetic fluctuations behind the bow shock, inside the magnetosheath, and near the magnetopause. It appears that the first Kramers- Moyal coefficient is linear and the second term is quadratic function of magnetic increments, which describe drift and diffusion, correspondingly, in the entire magnetosheath. This should correspond to a generalization of Ornstein-Uhlenbeck process. We demonstrate that the second order approximation of the Fokker-Planck equation leads to non-Gaussian kappa distributions of the probability density functions. In all cases in the magnetosheath, the approximate power-law distributions are recovered. For some moderate scales we have the kappa distributions described by various peaked shapes with heavy tails. In particular, for large values of the kappa parameter this shape is reduced to the normal Gaussian distribution. It is worth noting that for smaller kinetic scales the rescaled distributions exhibit a universal global scale-invariance, consistently with the stationary solution of the Fokker-Planck equation. These results, especially on kinetic scales, could be important for a better understanding of the physical mechanism governing turbulent systems in space and astrophysical plasmas.


INTRODUCTION
Turbulence is a complex phenomenon that notwithstanding progress in (magneto-)hydrodynamic simulations is still a challenge for natural sciences (Frisch 1995), and physical mechanisms responsible for turbulence cascade are not clear (Biskamp 2003). Fortunately, collisionless solar wind plasmas can be considered natural laboratories for investigating this complex dynamical system (Bruno & Carbone 2016). Fluctuations of magnetic fields play an important role in space plasmas, leading also to a phenomenon known as magnetic reconnection (e.g., Burlaga 1995;Treumann 2009).
Turbulent magnetic reconnection is a process in which energy can proficiently be shifted from a magnetic field to the motion of charged particles. Therefore, this process responsible the redistribution of kinetic and magnetic energy in space plasmas is pivotal to the Sun, Earth, as well as to other planets and generally across the whole Universe. Reconnection also impedes the effectiveness of fusion reactors and regulates geospace weather which can affect contemporary technology such as the Global Positioning System (GPS) navigation, modern mobile phone networks, including electrical power grids. Electric fields responsible for reconnection in the Earth's magnetosphere has been observed on kinetic scales by the Magnetospheric Multiscale (MMS) mission (Macek et al. 2019a,b). Certainly, reconnection in the magnetosphere is related to turbulence at small scales (Karimabadi et al. 2014).
Basically, in a Markov process, given an initial probability distribution function (PDF), the transition to the next stage can fully be determined. It is also interesting here that we can prove and demonstrate the existence of such a Markov process experimentally. Namely, without relying on any assumptions or models for the underlying stochastic process we are able to extract the differential equation of this Markov process directly from the collected experimental data. Hence this Markov approach appears to be a bridge between the statistical and dynamical analysis of complex physical systems. There is a substantial evidence based on statistical analysis that stochastic fluctuations exhibits Markov properties (Pedrizzetti & Novikov 1994;Renner et al. 2001). We have already proved that turbulence has Markovian features in the inertial range of hydromagnetic scales (Strumik & Macek 2008a,b). Admittedly, for turbulence inside the inertial region of magnetized plasma, the characteristic spectra should be close to standard Kolmogorov (1941) power-law type with exponent: −5/3 ≈ −1.67 (Kolmogorov 1941) and Kraichnan (1965) power-law spectrum with exponent: −3/2, but surprisingly, the absence of these classical spectra, especially on smaller scales, seems to be the rule.
Moreover, we have also confirmed clear breakpoints in the magnetic energy spectra in the Earth's magnetosheath, which occur near the ion gyrofrequencies just behind the bow shock, inside the magnetosheath, and before leaving the magnetosheath. Namely, we have observed that the spectrum steepens at these points to power exponents in the kinetic range from -5/2 to -11/2 for the magnetic field data of the highest resolution available within the MMS mission (Macek et al. 2018). Therefore, we would like to investigate the Markov property of stochastic fluctuations outside this inertial region of magnetized plasma on small scales, when the slopes are consistent with kinetic theory.
It should also be noted that based on the measurements of magnetic field fluctuations in the Earth's magnetosheath gathered onboard the MMS mission, we have recently extended these statistical results to much smaller scales, where kinetic theory should be applied (Macek et al. 2023).
Here we compare the characteristics of stochastic fluctuations behind the bow shock, inside the magnetosheath, and near the magnetopause. In this paper, we therefore present the results of our comparative analysis, where we check whether the solutions of the Fokker-Planck (FP) equation are consistent with experimental PDFs in various regions of the magnetosheath.
In Section 2, a concise description of the MMS mission and the analyzed data is provided, while the Section 3 outlines used stochastic mathematical and statistical methods. The vital results of our analysis are presented in Section 4, which demonstrates that the solutions of the FP equation are in good agreement with empirical PDFs. Finally, Section 5 emphasizes the significance of stochastic Markov processes in relation to turbulence in space plasmas, which exhibit a universal global scale-invariance across the kinetic domain.

DATA
The MMS mission, which begun in March 12, 2015 and is still in operation, delves into the connection and disconnection of the Sun's and Earth's magnetic fields. Four spacecraft (namely MMS 1 -MMS 4), which carry identical instrument suites, are orbiting near the equator to observe magnetic turbulence in progress. They are formed into a pyramid-like formation. Each satellite has an octagonal form that is around 3.5 meters in breadth and 1.2 meters in height. The satellites rotate at Three Revolutions Per Minute during scientific operations. Position data is supplied by ultra-precise GPS apparatus, while attitude is sustained by four stellar trackers, two accelerometers, and two solar sensors. All of the spacecraft have identical instruments to measure plasmas, magnetic and electric fields, and other particles with remarkably high (milliseconds) time resolution and accuracy. This allows us to figure out if reconnection takes place in an individual area, everywhere within a broader area simultaneously, or traversing through space. The MMS studies the reconnection of the solar and terrestrial magnetic fields in both the day and night sides of Earth's magnetosphere, which is the only place where it can be directly observed by spacecraft. In our previous studies we have observed reconnection in the Earth's magnetosphere involving small kinetic scales (Macek et al. 2019a,b).
We have further examined fluctuations of all components of the magnetic fields B = ( , , ) in the Geocentric Solar Ecliptic (GSE) coordinates, with the magnitude strength = |B| (square root of the sum of the squares of the field components), which have been taken from the MMS Satellite No. 1 (MMS 1), located just beyond the Earth's bow shock (BS). In this way, we have shown that magnetic fluctuations exhibits Markov character also on very small kinetic scales (Macek et al. 2023). Moreover, we have noticed that in all components the Markovian features are quite similar. Here, we would like to further investigate statistical properties of magnetic fluctuations in various regions of the magnetosheath. The spacecraft trajectories in the magnetosheath, in three different regions, namely: • (a) just behind the bow shock (BS), • (b) deep inside the magnetosheath (SH), and • (c) near the magnetopause (MP), which have been shown in Fig. 1 of (Macek et al. 2018).
Therefore, we would like to look at the measurements of the magnetic field strength = |B|, but now at various regions of the magnetosheath. To investigate magnetosheath stochastic fluctuations, now we have chosen the same three different time intervals samples, which correspond to Table 1 (List of Selected MMS 1 Interval Samples) of Ref. (Macek et al. 2018). In cases (a) and (c) of approximately 5 minutes and 1.8 minutes respective intervals, we use burst type observations from the FGM (FluxGate Magnetometer) sensor with the highest resolution (Δ ) of 7.8 ms (128 samples per second) with 37,856 and 13,959 data points, correspondingly. However, in the other case (b), between the bow shock and the magnetopause, where only substantially lower resolution 62.5-125 ms in survey mode (8-16 samples per second) data are available, we have a much longer interval lasting 3.5 h with 198,717 data points with Δ = 62.5 s. All of the data are publicly available on the website: http://cdaweb.gsfc.nasa.gov, which is hosted by NASA.
Admittedly, the gaps in time series, which commonly appear in the data gathered from space missions, can have a considerable impact on the conclusions that can be derived from statistical analysis based on experimental data. One of the powerful but simple tools used to cope with this problem is a linear interpolation method between points, which we have used, if necessary, to fill these gaps in the analyzed data sets. Therefore, in Fig. 1 on the upper side of each case (a) -(c) from left to right, we have depicted time series of the magnetic field strength = |B|. Whereas on the bottom side of each case, we have shown the respective power spectral density (PSD) obtained using the method proposed by Welch (1967). 28 .12.2015, 01:48:04.251 -01:52:59.997 28.12.2015, 06:19:00.043 -09:45:59.964 27.12.2016, 11:30:24.003 -11:32:13.052  Table 1 of Ref. (Macek et al. 2018).
The calculated average ion and electron gyrofrequencies are as follows: in case (a) = 0.25 Hz, and = 528 Hz; case (b) = 0.24 Hz and = 510 Hz; case (c) = 0.29 Hz and = 609 Hz (Macek et al. 2018). In addition, employing the hypothesis according to Taylor (1938), relating time and space scales in this way: = sw · , where is a spatial scale and sw is the mean velocity of the solar wind flow in the magnetosheath, we estimate characteristic inertial frequencies for ions and electrons: in case (a) = 0.55 Hz and = 24.5 Hz; case (b) = 0.41 Hz and = 18.1 Hz; case (c) = 0.45 Hz and = 20.1 Hz. We have marked these values on each graph of power spectral density. In case (a) the obtained spectral exponent is about −2.60 ± 0.06 somewhat steeper, before the = 24.5 Hz threshold and undoubtedly more steepen than the Kolmogorov (1941) (-5/3) or Kraichnan (1965) (-3/2) slopes. On the other hand, outside the inertial range of scales large spectral exponents has been reported from the Cluster multi-spacecraft mission (Sahraoui et al. 2009), the WIND data (Bruno et al. 2014), including the proposed explanation of nature of solar wind magnetic fluctuations on kinetic scales based on the missions (e.g., Lion et al. 2016;Roberts et al. 2016). Owing to unprecedented high 7.8 ms time resolution of magnetometer data in the MMS mission available in burst mode, we also see that in case (a) the slope is of −2.60 ± 0.06 (close to -5/2) above = 24.5 Hz. This is further followed by an even steeper spectrum with the slope of −5.59 ± 0.32 (close to -11/2 or -16/3). Because of a substantially lower survey data resolution of 62.5 ms in case (b) the spectrum with −2.24 ± 0.09 (≈ −7/3) is steeper than the Kolmogorov (1941) (−5/3) spectrum only after the visible breakpoint in the slope, which lies at = 0.12 Hz, i.e. near the ion gyrofrequency = 0.24 Hz, while more gentle slope of −0.77 ± 0.06 is observed before this breakpoint. Finally, in case (c), similarly as in case (a) using burst data, the spectral exponent of −2.75 ± 0.05 is again steeper before, and even more with the exponent −3.82 ± 0.06 (close to -7/2) after the observed breakpoint that lies at around the electron Taylor's (1938) shifted frequency = 20 Hz, as discussed by Macek et al. (2018). This shows that the observed stochastic nature of fluctuations in the sub-ion scale could be due to the interaction between coherent structures (Perrone et al. 2016(Perrone et al. , 2017, and a very high slope of -16/3 is possibly related to the dissipation of the kinetic Alfvén waves (e.g., Schekochihin et al. 2009).

METHODS OF DATA ANALYSIS
As usual, we use the fluctuations of the magnetic fields = |B|, which describe this turbulent system at each time > 0. Therefore, with a given time scale > 0 ∀ , one can typically define the increments of this quantity as follows: and, assuming an arbitrary > 0, it can be labeled as or for simplicity in the following sections. We assume that the fluctuations of increment in a larger time scale are transferred to smaller and smaller scales. Therefore, stochastic fluctuations may be regarded as a stochastic process in scale with the -point joint (transition) conditional probability density function denoted by ( 1 , 1 | 2 , 2 , . . . , , ). In this case, the conditional probability density function is defined by default as: with the marginal (unconditional) probability density function, ( , ), and the joint probability function, ( , ; , ), of finding the fluctuations at a scale and at a scale , for 0 < < . In the same way, we may construct the conditional probability densities for any longer sequences of increments .
The stochastic process is Markovian if the conditional probability function depends only on the initial values 1 and 2 at the time scales 1 and 2 , but not on 3 at the next larger scale 3 , and so on, i.e., for any = 1, . . . , we have: for 0 < 1 < 2 < . . . < . Basically, the Markov process can be determined by the initial conditional probability function ( 1 , 1 | 2 , 2 ). Strictly speaking, the future states of the process are conditionally independent of past states. Because of this relation, the conditional probabilities are also called transition probabilities, while the property of Eq. (3) is known as a memorylessness.
One of the generalizations of Eq. (3) is called the Chapman-Kolmogorov (CK) condition, which is given by the equation (Risken 1996): for 1 < ′ < 2 . This equation can be interpreted in the following way: the transition probability from 2 at a time scale 2 to 1 at a time scale 1 is the same as a product of the transition probability from 2 at a time scale 2 to ′ at a time scale ′ , and the transition probability from ′ at a time scale ′ to 1 at a time scale 1 , for all possible ′ 's. Let us emphasize here, that such a generalization is a necessary condition for a stochastic process to be the Markov process. Next, from the CK condition of Eq. (4), by using a standard series expansion, one can derive a corresponding Kramers-Moyal backward expansion with an infinite number of terms. Backward expansions are equations of evolution of probability ( , | ′ , ′ ), where we differentiate with respect to . This equation has the following differential form (Risken 1996, Section 4.2): where it is important to note that the differential symbol acts on both ( ) ( , ) and ( , | ′ , ′ ) coefficients. Since the solutions of the forward and backward KM equations are equivalent, then without loss of generality, we can label it as KM expansion. Formally, ( ) ( , ) are called KM coefficients, which in this way are defined as the limit at → ′ of -th power of conditional moments (see Risken 1996): Ideally, using the conditional moments ( ) ( , , ′ ), the KM coefficients can be evaluated, though they cannot be obtained directly from the analyzed data. While these conditional moments can be calculated from the empirical observations, the ( ) ( , ) coefficients can only be obtained by extrapolation in the limit → ′ according to Eqs. (6) and (7), but these formulae can not be applied explicitly. One of the popular extrapolation methods for this problem is a use of piecewise linear regression model with breakpoints. This is a type of regression models, which allows multiple linear models to fit to the analyzed data. The crucial objective of this method is an accurate estimation of a number of breakpoints. First, in order to estimate the best breakpoint position, we have evaluated every value within a specified interval and looked at the value of logarithmic transformation of the likelihood function (also known as log-likelihood function) of each adjusted model. Naturally, the highest value of this function provides the optimal breakpoint. Further, to select (and estimate) the best possible number of breakpoints of the segmented relationship, we have used the standard Akaike (1973) Information Criterion (AIC) and Bayesian Information Criterion (BIC) (Schwarz 1978). Nonetheless, the truly similar results are obtained when the lowest time resolution is taken. Thus, in our case, we have a simple approximation of the KM coefficients, which is given by: where a Δ is a given lowest time resolution of the time series. It is also interesting to note that ( ) ( , ) coefficients show the same dependence on as ( ) ( , , ′ ). This simplification substantially decrease the time required to obtain the results numerically. Now, in order to find the solution of Eq. (5), it is necessary to determine the number of terms of the right hand side (RHS) of this equation that need to be considered. According to Pawula's theorem, the KM expansion of a positive transition probability ( , | ′ , ′ ) may end after the first or second term (e.g., Risken 1996, Section 4.3). If it does not end after the second term, then the expansion must contain an infinite number of terms. On the other hand, if the second term is the last one, namely ( ) ( , ) = 0 for ⩾ 3, then the KM expansion of Eq. (5) leads to the following particular formula: with the well-known FP operator L FP ( , ) in the squared parenthesis (e.g., Risken 1996, Eqs. 5.1 and 5.2) governing the evolution of the probability density function ( , | ′ , ′ ) and is called the FP equation (also known as a forward Kolmogorov equation). It has been primarily used for the Brownian motion of particles, but now Eq. (9) defines a generalized Ornstein-Uhlenbeck process. Strictly speaking, this is a linear second-order partial differential equation of a parabolic type. By solving the FP equation, it is possible to find distribution functions from which any averages (expected values) of macroscopic variables can be determined by integration. If the relevant time-dependent solution is provided, this equation can be used to not only describe stationary features, but also the dynamics of systems. The first term (1) ( , ) and a second term (2) ( , ) > 0 determining the FP Equation (9) are responsible for the drift and diffusion processes, respectively. The former process accounts for the deterministic evolution of the stochastic process (as a function of and ). The latter process modulates the amplitude of the -correlated Gaussian noise Γ( ) (which is known as the Langevin force -the fluctuating force ( ) per unit mass ), that fulfills the normalization conditions: where is a Dirac delta function and ⟨Γ( )⟩ = 0 (see Risken 1996). Thus, in the equivalent approach another complementary equation arises: which is formally called the Langevin equation. Here we have used the Itô (1944) definition, that is missing a spurious drift (e.g., Risken 1996, Section 3.3.3), hence the drift coefficient (1) occurs directly, unlike in the Stratonovich (1968) definition. Admittedly, the Itô (1944) definition is more difficult to interpret and analyze, because of the new rules for integration and differentiation that must be used. Although, owing to a powerful apparatus, which is the Itô Lemma, it allows us to deal with stochastic processes analytically. Anyway, here again, all higher KM coefficients ( ) for ⩾ 3 are equal to zero. Note that the negative signs on the left hand side (LHS) of Eqs. (9) and (10) show that the corresponding transitions proceed backward to smaller and smaller scales. Next, because the differentiating in the FP operator in Eq. (9) should act on both the KM coefficients and the conditional probability density ( , | ′ , ′ ) by performing relatively simple transformations, it can be rewritten in the following expanded form (Risken 1996, Eq. 45.3a): Formally, Eq. (11) resulting from the FP Equation (9) is the second-order parabolic partial differential equation. It is also worth mentioning that this equation is the generalization of the case of thermal conductivity diffusion equation, which can be solved with the initial and boundary conditions ( , = 0| ′ , ′ = 0) = ( , ′ ) and ( = 0, | ′ = 0, ′ ) = 0, respectively, using the method of separation of variables. The solution of nonstationary FP Equation (11) can well be approximated numerically, i.e. by the difference method. The master curve for the probability density function ( , ) of Eq. (11) can readily be evaluated by the stationary solution ( , ) of Eq. (9), which is given by that results from comparing the LHS of Eq. (9) with zero.

RESULTS
In order to inspect processes responsible for stochastic fluctuations in space plasma, we have applied the methods described in Section 3 to small-scale in cases (a) and (c) and medium-scale in case (b) fluctuations of the magnetic field = |B| in the Earth's magnetosheath. In general, the approach presented in this paper could be applied under a few important conditions that should be tested as preliminary procedures (see Rinn et al. 2016). The first condition is that time series data must be stationary. If they were non-stationary, then the conditional moments given by Eq. (7) are not essentially meaningful. The second condition is that the process should be Markovian, i.e., the present state should only depend on the preceding state. The third condition is that the Pawula's theorem must hold, as discussed in Section 3.
Having this in mind, we have started with the brief analysis and description of the relevant time series and the corresponding graphs of power spectral densities (PSD). Next, we have checked stationarity of all analyzed time series (see, e.g., Macek 1998). To show that a Markov processes approach is suitable in our situation, we have moved forward to the verification of the necessary CK condition, through estimation of the KM coefficients, and then have checked the validity of the Pawula's theorem. This lets us to apply the reduced formula of the FP Equation (9), which describes evolution of the probability density function ( , ). Following our initial discussion, we must now verify whether the data time series under study is stationary. Generally, if a time series exhibits no trend, has a constant variance over time, and a consistent autocorrelation function over time, then it is classified as stationary. Such time series are also much easier to model. There are a variety of ways to evaluate this feature of any time series. One of such method is the Augmented Dickey & Fuller (1979) test. This test uses the following null and alternative hypotheses: H 0 : The time series is non-stationary, vs. H 1 : The time series is stationary. When the p-value is less than 0.05, then the null hypothesis can be rejected and it can be concluded that the time series is stationary. In fact, after performing such a statistical test, we have determined that in cases (a) and (b), the respective p-values are < 0.01, indicating that the null hypothesis can be rejected. Thus, these magnetic field strength = |B| time series are stationary. However, in case (c) where a much smaller data sample is available the p-value is equal to 0.154, hence we have failed to reject the null hypothesis. The result suggests that the time series is non-stationary and has some time-dependent structure with varying variance over time.
Once again, there are various methods of eliminating trends and seasonality, which define non-stationary time series. Trends can cause the mean to fluctuate over time, while seasonality can lead to changes in the variance over time. The most straightforward approach to address this issue is the differencing technique, a common and frequently used data transformation that is applied for making time series data stationary. Differencing is achieved by subtracting the previous observation from the current one. Following notation in Eq. (1), this can simply be written as ( ) = ( ) − ( − 1). To reverse this process, the prior time step's observation must be added to the difference value. The practice of computing the difference between successive observations is referred to as a lag-1 difference. The number of times that differencing is carried out is referred to as the order of differentiation. Fortunately, in our case (c), applying the lag-1 (order 1) difference operation has been sufficient to get rid of non-stationarity. The augmented Dickey-Fuller test has yielded a p-value of less than 0.01, thus the null hypothesis could be rejected, indicating that the analyzed = |B| time series is stationary.
We have used one of the exploratory data analysis approaches called unsupervised binning method (compare with normalized histogram method) to make bins (histogram's boxes) and to obtain the empirical conditional probability density functions ( 1 , 1 | 2 , 2 ), for 0 < 1 < 2 directly from the analyzed data. First, we have estimated the empirical joint PDF ( 1 , 1 ; 2 , 2 ) by counting the number of different pairs ( 1 , 2 ) on a 2-dimensional grid of equal-width data bins (small intervals). This bins integer should be neither too large, such that each bin no longer contains a significant quantity of points, nor too small, such that any dependency of the drift and diffusion coefficients on the state variable cannot be detected. Next, we have performed the normalization such that the integral over all bins is equal to 1 (note that the sum will not be equal to 1 unless bins of unity width are chosen). Similarly, the empirical one-dimensional PDF ( 2 , 2 ) can be estimated with the use of a one-dimensional grid of bins (and carrying-out the normalization), and the empirical conditional PDFs are obtained using Eq. (2) directly (in a numerical sense).
In such a way, we have found the empirical conditional probability density functions from the analyzed data, which are shown by red continuous contours in Fig. 2. They are compared here with the theoretical conditional PDFs that are solutions of the CK condition of Eq. (4) displayed by blue dashed contours, which are 2-dimensional representation of 3-dimensional data. Such a comparison is seen in Fig. 2 for the magnetic field increments , at the various scales: in cases (a)    subsequent isolines correspond to the following decreasing levels of the conditional PDFs, from the middle of the plots, for following magnetic field increments : case (a) 2, 1.1, 0.5, 0.3, 0.05, 0.01; case (b) 5, 1, 0.7, 0.45, 0.3, 0.22, 0.15, 0.1, 0.05; and case (c) 7, 3.3, 1.3, 0.3, 0.08, 0.06. This is rather evident that the contour lines corresponding to these two empirical and theoretical probability distributions are nearly matching for all three cases. Thus, it appears that the CK condition of Eq. (4) is sufficiently well satisfied. Next, in the corresponding Fig. 3, we have verified again the CK condition of Eq. (4). Intuitively speaking (and somehow informally), what we see in Fig. 2 is just a view 'from the top' of the 3-dimensional shape, while in Fig. 3 the orthogonal cuts are depicted. Again, we have compared these cuts through the conditional probability density functions for particular chosen values of parameter 2 , which can be seen at the top of each plot. As is evident, the cuts through the empirical probability density functions coincide rather well with the cuts through the theoretical probability density functions, providing good fits in all of the analyzed cases. Admittedly, only in case (b) for 2 = 0 [nT] the cuts    points deviate from the lines in tails, but it seem to be caused by the data outliers, which can eventually be further eliminated. It is necessary to mention that after such a comparison for different values of ( 1 , ′ , 2 ), we have found that the CK condition of Eq. To verify Pawula's theorem, which states that if the fourth-order coefficient is equal to zero, then ( ) ( , ) = 0, ⩾ 3, it is necessary to estimate the (1) ( , ), (2) ( , ) and (4) ( , ) coefficients using our experimental data. The standard procedure for calculating these values is to use an extrapolation method such as a piecewise linear regression to estimate the respective limits in Eq. (6). However, as already mentioned in Section 3, the similar results are obtained by simplifying the problem of finding these coefficients, by using Eq. (8), which enables us to estimate these values using the adequately scaled ( ) ( , , ′ ) coefficients. In our situation, the time resolution Δ is equal to 7.8 ms in case (a) and (c), while in case (b) it is 62.5 ms. Thus, given the conditional probabilities ( 1 , 1 | 2 , 2 ), for 0 < 1 < 2 , we have calculated these central moments directly from Eq. (7), using the obtained empirical data by counting the numbers ( ′ , ) of occurrences of two fluctuations ′ and . Given that the errors of ( ′ , ) might be simply determined by √︁ ( ′ , ), then, in a similar way, it is possible to calculate the errors for the conditional moments ( ) ( , , ′ ). Consequently, scaling these values according to Eq. (8), we have obtained the empirical KM coefficients. By examination of the ( ) ( , , ′ ) and ( ) ( , ) coefficients, we can observe that they both exhibit the same dependence on . The results of this analysis are shown in Fig. 4, where on the upper part we have depicted the first order coefficient depending on , while at the bottom we have shown the second and fourth order coefficients depending on , for all three cases (a), (b), and (c). Moreover, for each case, we have provided the calculated confidence intervals (error bars). It is demonstrated, that the fit for (1) ( , ) coefficient is a linear function of and for (2) ( , ) is a quadratic function of , for Δ = 0.0078 s in cases (a) and (c), and Δ = 0.0625 s in case (b). In fact, we have checked that the same fits are reasonable up to even 150Δ for all three analyzed cases. This means that in this instance, there should be no difficulties with fitting the polynomials for different Δ .
As seen at the bottom part of Fig. 4 of cases (a) and (c), it is evident that the Pawula's theorem is clearly satisfied. On the other hand, in case (b) it might be not so obvious. For instance, for ≈ −6.2 nT, we can see that the value of (4) ( , ) is somewhat greater than zero. In this case, we can use the somewhat weaker version of this theorem, which states that it is sufficient to check if (4) ( , ) ≪ [ (2) ( , )] 2 , for all (see Risken 1996;Rinn et al. 2016). Thus, in this situation, we have [ (2) ( , )] 2 ≈ 1225, which is significantly larger than (4) ( , ) ≈ 1, for ≈ −6.2 nT. Therefore, it is reasonable to conclude that the Pawula's theorem is sufficiently well fulfilled in all of the analyzed cases. Hence we can assume that the Markov process is described by the FP Equation (9).
In order to find the analytical solution of the FP Equation (9), we have proposed certain approximations of the lowest order KM coefficients. As previously discussed (see Fig. 4), it is straightforward that (1) ( , ) exhibits a linear dependence, whereas (2) ( , ) displays a quadratic dependence on . Consequently, it is reasonable to assume the following parametrization: where the relevant parameters 1 > 0, 2 > 0, and 2 > 0 depend on temporal scale > 0. Moreover, it appears that all of these parameters exhibit a power-law dependence on temporal scale : where the values for all of the logarithmized parameters , , ∈ R, as well as the , , ∈ R are given in Table 1.
It is important to emphasize that the functional dependencies of the coefficients 1 ( ), 2 ( ), and 2 ( ) on given by Eq. (14) are merely parametrizations of the empirical results. In fact, here power-laws have been selected, because they have adequately described the observed values with sufficient accuracy. Nevertheless, some alternative theoretical analyses may lead to slightly different functional dependence (see Renner et al. 2001). Admittedly, it turned out that the values of the fitted parameters can slightly be different from those that fit exactly the solution of the FP Equation (9). Renner et al. (2001) has also highlighted the asymmetry of the fit (2) ( , ) on , which is also present in our analysis (especially in case (c), and to a lesser degree in case (a)).
The obtained fits to the MMS observations in the magnetosheath are depicted in Fig. 5, for each case (a), (b), and (c), showing the dependence of KM coefficients parameters on scale > 0. Since our data contain a multitude of relatively low values and a few exceedingly large values, which would render a linear graph rather unreadable, instead of using a standard linear graph, we have decided to employ logarithmic scales for both the vertical and horizontal axes (so called: log − log plot). Such a procedure is rather straightforward. For example, for the first row of Eq. (14), taking the logarithm of both sides one obtains: log( 1 ( )) = log( ) + log( ), which is a special case of a linear function, with the exponent corresponding to the slope of the line. The value of log( ) corresponds to the intercept of a log( 1 ( ))-axis, while the log( )-axis is intercepted at log /(− ). We have opted for this approach to enhance the clarity of the presentation. Therefore, since we have used both the logarithmic scales the respective power-laws appear as straight lines in Fig. 5. Similarly, the graphical representations for all the parameters 1 , 1 , and 2 of Eqs. (13) and (14), which we have provided, are helpful for identifying correlations and determining respective constants , , and < 0, > 0, < 0 in Table 1 (cf. Macek et al. 2023). After performing a careful analysis of the MMS magnetic field magnitude data, our findings indicate that the power-law dependence is applicable for the values of: ≲ 100Δ = 0.78 s in case (a); ≲ 150Δ = 9.375 s in case (b); ≲ 50Δ = 0.39 s in case (c), and for some larger scales, say ≳ , the shapes of the probability density functions appear to be closer to Gaussian. However, despite the satisfactory results obtained at these small kinetic scales, a more intricate functional dependence (possibly polynomial fits) is characteristic for much higher scales, in particular, in the inertial domain (Strumik & Macek 2008a,b).
As a result of our investigations, we are able to obtain analytical stationary solutions ( ) given by Eq. (12) following from the FP Equation (9), which can be expressed by a continuous kappa distribution (also known as Pearson's type VII distribution), which exhibits a deviation from the normal Gaussian distribution. The probability density function (PDF) of this distribution is of a given form: where, for 2 ( ) ≠ 0, 0 ( ) ≠ 0, we have a shape parameter = 1 + 1 ( )/[2 2 ( )] and 2 0 = 2 ( )/[ 2 ( ) + 1 ( )/2], while 0 = (0) satisfies the normalization ∫ ∞ −∞ ( ) = 1. By substituting ( ) to this integral we find that: (τ) . These values of would correspond to the nonextensivity parameter of the generalized (Tsallis) entropy = 1 − 1/ (e.g., Burlaga & Viñas 2005). In our case this is given by = 1 ( ) 1 ( )+2 2 ( ) and is equal to 0.341 in case (a), 0.273 in case (b), and 0.716 in case (c). The extracted values of the and parameters provide robust measures of the departure of the system from equilibrium. We see that these values are similar to ∼ 0.5 for ∼ 2 reported for the Parker Solar Probe (PSP) data by Benella et al. (2022). Now, by using the system of Eqs. (13) with Eq. (9), we have arrived at the following formula (Macek et al. 2023): This implies that the FP Equations (11) and (18) are expressed in terms of a second-order parabolic partial differential equation. Thus, through the implementation of the numerical Euler integration scheme, which has been verified for stationary solution ( , ) = 0, we are able to successfully solve the non-stationary FP equation numerically. Our results are in line with those obtained by Rinn et al. (2016) using the statistical modeling package in programming language R. Figure 6 shows the findings resulting from our analysis based on the MMS data. Here we compare the solutions of the FP Equation (9) with the empirical probability density functions of ( , ): (a) near the bow shock (BS), (b) inside the magnetosheath (SH), and (c) near the magnetopause (MP) for various scales (not greater than ). The displayed plotted curves, in each case, are as follows: the stationary solution (denoted by open circles), the non-stationary solutions (marked with dashed lines), and the empirical PDFs (indicated with various colored continuous lines).
Further, in cases (a) and (c) the corresponding time scales are: = 0.0078, 0.04, 0.078, 0.12, 0.2, 0.39, and 0.78 s; whereas in case (b) these scales are equal: = 0.0625, 0.3125, 0.625, 0.9375, 1.5625, 3.125, and 9.375 s. The corresponding curves are shifted in the vertical direction from bottom to top for even better clarity of presentation. It is also worth noting that we have used the semi-logarithmic scale , what is useful when dealing with data that covers a broad range of values. On this scale, the vertical scale is logarithmic (base 10) axis, which means that the separation between the ticks on the graph is proportional to the logarithm of PDF, while the horizontal -axis is a standard linear scale, and the ticks are evenly spaced.
What is important to note from this picture are the peaked leptokurtic shapes of PDFs and corresponding stationary solutions. Namely, in case (a) the peak (with large kurtosis) is present for scales up to ∼ 0.5 s; in case (b) up to about ∼ 3 s; and in case (c) up to ∼ 0.25 s. For these levels selected for each case the PDF becomes closer to Gaussian (i.e., approximately parabolic shape on the graph with the semi-logarithmic scales), as expected for large values of . In case (c) we can see more jumps in fluctuations, i.e., the curves are not so smooth. Fluctuations are quite evident in both the empirical curves and the theoretical solutions, so it seems that some numerical noise is present in the tails of the PDFs. Admittedly, reducing noise is a tricky issue, although the easiest way is to artificially smooth using the simple moving average. Therefore, we have tried this procedure for = 1, 2, 3 steps and it has appeared that the = 3 choice is sufficient. Figure 7 depicts finally the probability density functions (PDFs) of fluctuations of the strength of the magnetic field rescaled by the standard deviations ( ) in the following way: PDF( ) −→ ( ) · PDF( ).
In this way, we can define a master curve for the shape of the PDFs. Again, we have used the logarithmic scale on the vertical axis. We also see that the rescaled curves are consistent with the stationary solutions of Eq. (15), as marked with open circles in Fig. 6. It should be noted that all the curves in Fig. 7 are very close to each other for small scales. However, for larger = 50 or 100Δ these shapes deviate from the master curve and naturally tend to the well known Gaussian shape. We see that the shape of the PDFs obtained from the MMS data exhibits a global scale-invariance in the magnetosheath up to scales of ∼ 2 s. A similar collapse has also been reported with the PSP data at subproton scales A collapse of probability density functions of (compare Fig. 6), which are scaled by the corresponding standard deviations (see Eqs. (19) and (20)), for the small time scales stopping at approximately: ∼ 0.4 s. in case (a), ∼ 2.0 s in case (b), and ∼ 0.25 s in case (c). (Benella et al. 2022). Figure 7 shows that fluctuations in the magnetosheath are described by a stochastic process. Admittedly, the mechanism of generation of these magnetic fluctuation at small kinetic scale is not known, but the this results suggest some universal characteristics of the processes. An alternative point of view has recently been proposed by Carbone et al. (2022).

CONCLUSIONS
Following our studies in the space plasmas at large inertial scales (Strumik & Macek 2008a,b), we have examined time series of the strength of magnetic fields in different regions of the Earth's magnetosheath, where the spectrum steepens at subproton scales (Macek et al. 2018). With the highest resolution available on the MMS, the data samples just after the bow shock and near the magnetopause are stationary and for somewhat lower resolution deep inside the magnetosheath the deviations from stationarity are small and could well be eliminated. Basically, in all these cases the stochastic fluctuations exhibits Markovian features. We have verified that the necessary Chapman-Kolmogorov condition is well satisfied, and the probability density functions are consistent with the solutions of this condition.
In addition, the Pawula's theorem is also well satisfied resulting in the Fokker-Planck equation reduced to drift and diffusion terms. Hence, this corresponds to the generalization of Ornstein-Uhlenbeck process. Further, the lowest Kramers-Moyal coefficients have linear and quadratic dependence as functions of the magnetic field increments. In this way, the power-law distributions are well recovered throughout the entire magnetosheath. For some moderate scales we have the kappa distributions described by various peaked shapes with heavy tails. In particular, for large values of the kappa parameter these distributions are reduced to the normal Gaussian distribution.
Similarly as for the PSP data, the probability density functions of the magnetic fields measured onboard the MMS rescaled by the respective standard deviations exhibit a universal global scale-invariance on kinetic scales, which is consistent with the stationary solution of the Fokker-Planck equation. We hope that all these results, especially those reported at small scales, are important for a better understanding of the physical mechanism governing turbulent systems in space and laboratory.