© The Authors 2025

1. Introduction

The term “solar wind” refers to a low-density, high-speed stream of charged particles that emanates from the Sun and permeates the heliosphere. With the advent of space exploration, numerous missions have been dedicated to measuring the parameters and fields of the solar wind. These missions offer a unique opportunity to gain valuable insight into the nature and behavior of this phenomenon.

Many processes occurring in the solar wind are inherently nonlinear. To study these processes, it is essential to consider the temporal variations in the characteristics of the solar wind. While most traditional methods of classical physics are primarily suited for stationary or quasi-stationary phenomena, the analysis of dynamic regimes, fluctuations, and self-similar scaling requires the application of nonlinear dynamics. In this context, the development of methods based on fractal geometry to describe the temporal behavior of the solar wind is of particular interest.

Turbulence is a fundamental feature of the solar wind and is commonly observed in both neutral fluid and plasma flows. It acts as a mechanism for transferring energy from large scales, where it is initially injected, to smaller scales. At these microscopic levels, dissipative and dispersive processes convert the transferred energy into other forms, such as heat or particle acceleration. The resulting fluctuations exhibit power-law scaling, a behavior that is derived from the inherent scale invariance and self-similarity of the system. This power-law spectrum has been consistently observed in the magnetic-fieldfluctuations of the solar wind since the early days of space exploration, making it a well-established characteristic (Coleman 1968; Bruno & Carbone 2013). Fat-tailed distributions have been observed in various solar wind parameters; for example, numerous studies have analyzed fluctuations in the heliospheric magnetic field strength, B (Burlaga 1991). These analyses have led to the identification of three key implications: (i) fat-tailed (non-Gaussian) distributions, (ii) slow relaxation processes indicative of long-term memory effects and, (iii) multifractal structure in the time series.

(i) In the case of magnetic field strength fluctuations, fat-tail distributions emerge as a result of intermittent energy transfer in turbulent systems, which increases the probability of extreme fluctuations. Fat-tailed behavior is also a prominent feature in the velocity distribution profiles of plasma particles (Maksimovic et al. 1997a; Shan & Saleem 2017; Yoon et al. 2024). In particular, the suprathermal tails observed in electron distributions in the solar wind are closely associated with the energetic particles associated with energetic particle events. The nonuniform nature of energy transfer in turbulent systems causes energy to concentrate in localized spatial regions, leading to the emergence of highly energetic fluctuations. As a result, the tails of the probability distribution functions (PDFs) become populated by these energetic particles, giving rise to the so-called long or fat-tailed distributions. Electron distribution functions in the solar wind consistently exhibit three distinct components: a thermal core and a suprathermal halo, both present at all pitch angles, and a sharply field-aligned “strahl” component, typically directed anti-sunward (Štverák et al. 2009). Although Coulomb collisions can account for the relative isotropy of the core population, the origin of the halo, particularly its sunward-directed portion, remains poorly understood (Maksimovic et al. 2005). Furthermore, non-Gaussian distributions have been reported in studies of electron temperature anisotropy in the solar wind (Štverák et al. 2008). Higher statistical moments become particularly relevant when the distribution exhibits heavy tails, especially in regimes where fluctuations are comparable to or exceed the average value of the measured quantity. In such regimes, the statistical moments – mean, variance, and higher-order moments – become essential tools for characterizing the system. For instance, when a variable exhibits no fluctuations, its distribution collapses into a Dirac delta function, making the mean value a sufficient description. However, with small fluctuations, distributions tend to approximate a Gaussian shape, and the variance becomes a relevant measure. When fluctuations are even larger – as in the case of fat-tailed distributions – higher-order statistics are required for proper characterization. Several works (e.g., Burlaga & F.-Viñas 2005a) used fourth-order moments (kurtosis) to demonstrate non-Gaussian behavior since a Gaussian distribution would exhibit a kurtosis of exactly 3. This highlights the importance of employing a more detailed statistical framework, such as the one used in the present analysis. This kind of distribution was also observed in the heliospheric magnetic field strength fluctuations (Burlaga & F.-Viñas 2004a,b, 2005a; Burlaga & Ness 2009; Burlaga et al. 2024a,b).

(ii) There is growing evidence that the transition to a quasi-equilibrium, non-Gaussian state in the solar wind involves inherently slow relaxation processes. Physically, a relaxation process describes how a system evolves toward a steady or equilibrium state after being disturbed. In classical systems, this evolution is often exponential and rapid. However, in complex or turbulent systems like the solar wind, relaxation can be much slower and follow non-exponential (e.g., power-law) decay patterns. These slower dynamics reflect the presence of long-range correlations and memory effects in the plasma, suggesting that the system retains information about its past states. The timescale and nature of these relaxation processes are central to understanding how turbulence and energy transfer operate across different scales in space plasmas (Servidio et al. 2014; Verscharen et al. 2019). As suggested by both hydrodynamic theory and recent magnetohydrodynamic numerical simulations, these relaxation processes can occur during the turbulent cascade and manifest as localized patches with equilibrium-like configurations. The coupling of processes across multiple scales plays a crucial role in shaping the global dynamics and thermodynamics of the solar wind. In particular, the presence of slow relaxation processes is often associated with the emergence of fat-tailed distributions (Zamora & Tsallis 2022).

(iii) The solar wind is a highly turbulent medium, exhibiting strong field fluctuations on a broad range of scales. These include an inertial range where a turbulent cascade is believed to be active. Notably, the solar wind cascade displays intermittency, although the degree of intermittency may vary depending on solar wind conditions. Intermittency can be interpreted as a manifestation of the multifractal nature of the turbulent cascade. A multifractal structure in the magnetic field strength (B) has been observed at various heliocentric distances and across different phases of the solar cycle (Burlaga 1991, 2004; Burlaga et al. 2003). The foundational theory of multi-fractals has been explored extensively in the literature; see, for example, Mandelbrot (1972) and Anselmet et al. (1984). The origin of multi-fractality in the solar wind can be attributed to the extension of intermittent turbulence to larger spatial scales at greater distances from the Sun, or it may arise from the nonlinear evolution and interaction of large-scale structures such as corotating streams, ejecta, and shocks. Although solar wind plasma is often treated as almost incompressible, observed correlations between velocity, temperature, and density (Elliott et al. 2016; Borovsky et al. 2021) have raised the question of whether similar nonlinear or multifractal structures might also be present in proton density. In fact, spectral analysis has revealed that proton density fluctuations exhibit a Kolmogorov-like power-law behavior (Shaikh & Zank 2010; Chen et al. 2011). More recently, small-scale fluctuations in the solar wind proton density have been shown to exhibit multifractal properties (Sorriso-Valvo et al. 2017), highlighting the need for different intermittency measures to fully characterize the small-scale cascade.

The paper is structured as follows. In Sect. 2 we present the theoretical background of non-extensive statistical mechanics and its relation to multifractal structures, fat-tailed distributions, and slow relaxation processes. Section 3 details the methodology used for extracting the q-triplet parameters from solar wind proton density data. The results of the 17-year data analysis and the validation of the q-triplet are presented in Sect. 4. Finally, Sect. 5 offers a discussion of the implications of our findings and outlines future directions for research.

2. Multi-fractals, fat-tail distributions, and slow relaxation processes under the view of non-extensive statistics

In this section we summarize the key theoretical concepts that form the basis of our data analysis methodology. We first present an overview of non-extensive statistical theory. Based on these theoretical foundations, we then describe the data analysis methodology and the algorithm employed to produce the novel results, which are discussed in detail in the next section.

The statistical theory of Boltzmann and Gibbs (BG) is grounded in the molecular chaos hypothesis, which assumes that the system exhibits ergodic motion in its microscopic phase space. In other words, the system can explore all microscopic states allowed with equal probability. In such cases, the probability distributions are Gaussian, and the observed time series exhibit fluctuations consistent with normal diffusion processes. Equilibrium dynamics corresponds to physical states characterized by uncorrelated or weakly correlated noise.

In contrast, nonequilibrium nonlinear dynamics can exhibit strong, long-range correlations. In such regimes, Gaussian statistics are inadequate to describe the observed behavior, as the underlying phenomena follow non-Gaussian statistics and violate the assumptions of the classical central limit theorem and the law of large numbers (Umarov et al. 2008, 2010). The standard BG statistical theory relies on two foundational assumptions: ergodicity and thermodynamic equilibrium. However, in systems where the dynamics are chaotic, exhibit sensitivity to initial conditions, possess memory effects, or involve long-range interactions, these assumptions no longer hold. As a result, the applicability of BG statistics is limited in such contexts. Specific theoretical difficulties on these kinds of systems are related to the fact that the parts interact with many others at long distances, so it is impossible to cut the system into almost independent pieces. Therefore, there is no distinction between bulk and surface, and consequently these systems are non-additive and non-ergodic (phase-space is not occupied uniformly). As a result, a new kind of statistic is necessary.

Since the early 1990s, non-extensive statistical mechanics has been applied in a wide range of scientific fields, demonstrating remarkable versatility and yielding multiple applications (Wilk & Włodarczyk 2000; Gell-Mann & Tsallis 2004; Tsallis 2009a,b; Vignat & Plastino 2009). It has proven particularly useful in the context of astrophysics (Plastino & Plastino 1993; Chavanis & Sommeria 1998; Scarfone et al. 2008; Sahu & Tribeche 2012; Rosa et al. 2013; Pavlos et al. 2018; Zamora et al. 2018, 2020). In particular, it has been found that the non-Gaussian distributions of magnetic field strength increments and other solar wind parameters are accurately described by the q-Gaussian distributions predicted by non-extensive statistical mechanics (Burlaga & F.-Viñas 2004a,b, 2005a; Burlaga et al. 2020).

Non-extensive statistical mechanics is based on a generalized measure of entropy S_q introduced in Tsallis (1988). S_q is defined as

$$\begin{array}{c}\hfill S_q=\frac{k_B}{q-1}[1-{\displaystyle \int p}{(x)}^{q}dx],\end{array}$$(1)

where p is the probability, and q is the non-extensivity parameter. For q = 1, the non-extensive entropy reduces to the standard BG entropy. The q-logarithm function is defined as

$$\begin{array}{c}\hfill {ln}_q(x)=\frac{x^{1-q}-1}{1-q},x>0.\end{array}$$(2)

It is easy to verify that ln_q = 1(x) = ln(x). The q-logarithm satisfies the following property:

$$\begin{array}{c}\hfill {ln}_q(x_A{x}_B)={ln}_q(x_A)+{ln}_q(x_B)+(1-q){ln}_q(x_A){ln}_q(x_B).\end{array}$$(3)

This function generalizes the natural logarithm. It follows directly that S_q can be expressed as S_q = k_B∫ln_q(ρ) dx. This expression resembles the BG entropy.

The inverse function of Eq. (2) is defined as the q-exponential function, given by

$$\begin{array}{c}\hfill e_q(x)={[1+(1-q)x]}_{+}^{1/(1-q)}.\end{array}$$(4)

This function generalizes the standard exponential: if q = 1, then e_q = 1(x) = e^x. The notation [ ]₊ means that the function is defined so that it vanishes for negative arguments inside the brackets, that is, [x]₊ = max(x, 0).

Within the non-extensive theory framework, three key features, namely, non-Gaussian distributions, slow relaxation processes, and multifractal structures, are interconnected through the so-called q-triplet. This concept was first introduced in Tsallis (2004, 2005), providing a unifying framework for describing complex, nonequilibrium systems such as the solar wind.

The q-triplet has proven to be a valuable tool for analyzing time series in atmospheric and space plasma environments. It has been applied to the study of solar activity using the AE and D_st indices (Gopinath et al. 2018), sunspot dynamics (Pavlos et al. 2012a), nonlinear analysis of the solar flare index (Karakatsanis et al. 2013), magnetospheric self-organization processes (Pavlos et al. 2012b), and nonequilibrium phase transitions in solar wind plasma dynamics during calm and shock periods (Pavlos et al. 2015).

As evidenced by the bibliography cited so far, substantial progress has been made in this field, especially in recent years. However, the PDFs of the solar wind parameters, turbulence, and transport of energetic particles remain open questions to this day (Viall & Borovsky 2020).

Empirically derived non-Gaussian distributions are becoming increasingly prevalent in space physics, as the power-law nature of various suprathermal tails is combined with more classical quasi-Maxwellian cores. In fact, q-Gaussian distributions have been used in plasma sciences long before under the name kappa distributions, which were independently proposed (Maksimovic et al. 1997b; Livadiotis 2016; Yoon 2019; Lazar & Fichtner 2021; Louarn et al. 2021). However, it can be shown that the two are equivalent through a suitable transformation (Livadiotis & McComas 2009). Nevertheless, the Tsallis statistical framework provides a set of mathematical and conceptual tools that go far beyond a mere modification of the distribution, making its implementation highly enriching for plasma theory in atmospheric and space environments. These non-Gaussian distributions arise naturally within the framework of non-extensive statistical mechanics, which offers a robust theoretical foundation for describing and analyzing complex systems out of equilibrium. Given the strong correspondence between empirically observed non-Gaussian distributions and the predictions of non-extensive statistics, the full suite of non-extensive statistical tools becomes available to the space physics community for investigating the non-Gaussian characteristics of particle and energy distributions observed in space (Livadiotis & McComas 2009). Moreover, the applicability of these methods extends beyond the solar wind. For example, Tsallis statistics have been shown to be effective in studying ionospheric plasma (Chernyshov et al. 2015; Ogunsua 2018), auroral glow (Chernyshov et al. 2024), and magnetospheric dynamics (Pavlos et al. 2011; Gopinath et al. 2018).

2.1. Quasi-stationary attractors and the q_stat parameter

Contrary to BG statistical mechanics, where the function of energy describing a thermal equilibrium state is characterized by a Gaussian function, a correlated quasi-equilibrium physical process can be described by the following nonlinear differential equation (Tsallis 2009a):

$$\begin{array}{c}\hfill \frac{d(p_{i}Z)}{d{E}_i}=-\beta {(p_{i}Z)}^{q_{\mathrm{stat}}},\end{array}$$(5)

with solution

$$\begin{array}{c}\hfill p_i=\frac{e_{q_{\mathrm{stat}}}^{-\beta E_i}}{Z},\end{array}$$(6)

where

$$\begin{array}{c}\hfill \beta =\frac{1}{k_{B}T},Z={\displaystyle \sum _j}{e}_{q_{\mathrm{stat}}}^{-\beta E_j}.\end{array}$$(7)

The PDF is then given by

$$\begin{array}{c}\hfill p(x)\propto {[1-(1-q_{\mathrm{stat}})\beta x^2]}^{\frac{1}{1-q_{\mathrm{stat}}}},\end{array}$$(8)

for continuous variables. The above distribution function is the so-called q-Gaussian function, and corresponds to the attracting stationary solution associated with the nonlinear dynamics of the system. The stationary solutions p(x) describe the probabilistic nature of the dynamics in the attractor set in the phase space. The stationary parameter, q_stat, varies accordingly as the attractor changes.

2.2. Relaxation processes and the q_rel parameter

Boltzmann–Gibbs statistics are associated with the exponential relaxation of macroscopic quantities to thermal equilibrium, ie., one expects an exponential decay with a relaxation time τ. If ΔS denotes the deviation of entropy from its equilibrium value (S₀), then the probability of a proposed fluctuation is given by

$$\begin{array}{c}\hfill p\sim exp(\mathrm{\Delta }S/k_B).\end{array}$$(9)

At the macroscopic level, the relaxation toward equilibrium of a dynamical observable O(t), which describes the system’s evolution in phase space, can be modeled by the general form:

$$\begin{array}{c}\hfill \frac{d\mathrm{\Omega }}{\mathit{dt}}\simeq -\frac{1}{\tau }\mathrm{\Omega },\end{array}$$(10)

where

$$\begin{array}{c}\hfill \mathrm{\Omega }(t)\equiv \frac{[O(t)-O(\infty )]}{[O(0)-O(\infty )]}\end{array}$$(11)

is a normalized measure of the deviation of O(t) from its stationary state value. Under the non-extensive generalization, the standard exponential relaxation process is replaced by a meta-equilibrium formulation governed by

$$\begin{array}{c}\hfill \frac{d\mathrm{\Omega }}{\mathit{dt}}=-\frac{1}{\tau }{\mathrm{\Omega }}^{q_{\mathrm{rel}}},\end{array}$$(12)

where q_rel characterizes the degree of non-extensivity in the relaxation process. The solution to this equation is

$$\begin{array}{c}\hfill \mathrm{\Omega }(t)=e_{q_{\mathrm{rel}}}^{-t/\tau },\end{array}$$(13)

where e_q^x is the q-exponential function.

2.3. Sensibility to initial conditions and the q_sens parameter

In BG statistical mechanics, systems typically exhibit exponential sensitivity to initial conditions. This behavior, known as strong chaos, is characterized by exponential divergence of nearby trajectories and quantified by one or more positive Lyapunov exponents (Ott 2002).

In contrast, non-extensive statistical mechanics is associated with q-exponential sensitivity to initial conditions, a hallmark of weak chaos. This regime is described by a q-exponential growth governed by the non-extensivity parameter q_sens.

The entropy production process is intimately connected with the structure of the system’s attractor in phase space. This structure can be characterized by its multi-fractality and by thesensitivity to initial conditions, which can be modeled by the following differential equation:

$$\begin{array}{c}\hfill \frac{d\xi }{\mathit{dt}}={\lambda }_1\xi +({\lambda }_q-{\lambda }_1){\xi }^{q_{\mathrm{sens}}},\end{array}$$(14)

where ξ(t) quantifies the divergence between nearby trajectories and λ₁ is the largest Lyapunov exponent. For λ₁ > 0 (λ₁ < 0), the system is strongly chaotic (regular), while for λ₁ = 0 it is at the edge of chaos. ξ(t) is defined through

$$\begin{array}{c}\hfill \xi \equiv \underset{\mathrm{\Delta }x(0)\to 0}{lim}\frac{\mathrm{\Delta }x(t)}{\mathrm{\Delta }x(0)},\end{array}$$(15)

with Δx(t) representing the distance between neighboring trajectories in phase space (Tsallis 2002).

The solution to Eq. (14) is given by

$$\begin{array}{c}\hfill \xi (t)={[1-\frac{{\lambda }_q}{{\lambda }_1}+\frac{{\lambda }_q}{{\lambda }_1}{e}^{(1-q_{\mathrm{sens}}){\lambda }_{1}t}]}^{\frac{1}{1-q_{\mathrm{sens}}}}.\end{array}$$(16)

This expression captures the nonlinear sensitivity of the system to initial conditions, and the parameter q_sens serves as a quantitative measure of the degree of deviation from standard exponential sensitivity.

According to Lyra & Tsallis (1998), the scaling properties of the most rarefied and most concentrated regions of multifractal dynamical attractors can be used to estimate the divergence ξ of nearby orbits, according to the first-order approximation:

$$\begin{array}{c}\hfill \xi =e_{q_{\mathrm{sens}}}^{{\lambda }_{q}t}={[1+(1-q_{\mathrm{sens}}){\lambda }_{q}t]}^{\frac{1}{1-q_{\mathrm{sens}}}}.\end{array}$$(17)

If the Lyapunov exponent λ₁ ≠ 0 then q_sens = 1 (strongly sensitive if λ₁ > 0, strongly insensitive if λ₁ < 0). If the Lyapunov exponent λ₁ = 0 (weakly sensitive) then q_sens < 1.

2.4. The q-triplet

Consider the three distinct features of nonlinear systems discussed earlier. The set (q_stat, q_rel, q_sens) constitutes what is known as the q-triplet (also occasionally referred to as the q-triangle Gell-Mann & Tsallis 2004). The values of the q-triplet characterize the attractor set of the dynamics in phase space. In the case of equilibrium (i.e., BG statistics), the q-triplet takes the values (q_stat = 1, q_rel = 1, q_sens = 1).

These indices are interrelated, as they all arise from the particular way in which the system explores its phase space. In the case of the solar wind, the following relationships are expected to hold under the framework of non-extensive statistical mechanics (Gazeau & Tsallis 2019):

$$\begin{array}{cc}\hfill \frac{1}{q_{\mathrm{rel}}-1}& =\frac{1}{q_{\mathrm{sens}}-1}+1,\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill \frac{1}{q_{\mathrm{stat}}-1}& =\frac{1}{q_{\mathrm{sens}}-1}+2.\hfill \end{array}$$(19)

Hence, only one of the q-triplet indices is independent. These relationships arise from the generalized thermodynamic structure applicable to systems that exhibit long-range correlations, memory effects, or fractal space-time constraints – features commonly found in turbulent plasmas. The solar wind, particularly in its turbulent and intermittent regimes, provides a natural testbed for such nonequilibrium conditions. However, the exact applicability of these relations may depend on solar wind conditions (e.g., fast vs. slow wind and shocked vs. quiet intervals), and understanding their variability across these regimes remains an active area of research. The proposed values of the q-triplet for the solar wind, based on the analysis in Burlaga & F.-Viñas (2005b) and Burlaga & Ness (2013), are: $q_{\mathrm{stat}}=\frac{7}{4}$, q_rel = 4, and $q_{\mathrm{sens}}=-\frac{1}{2}$. If we define the auxiliary quantities as

$$\begin{array}{cc}\hfill a_{\mathrm{sens}}& :=\frac{1}{1-q_{\mathrm{sens}}}=\frac{2}{3},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill a_{\mathrm{stat}}& :=\frac{1}{q_{\mathrm{stat}}-1}=\frac{4}{3},\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill a_{\mathrm{rel}}& :=\frac{1}{q_{\mathrm{rel}}-1}=\frac{1}{3},\hfill \end{array}$$(22)

we also verify that

$$\begin{array}{c}\hfill a_{\mathrm{rel}}+a_{\mathrm{stat}}-a_{\mathrm{sens}}=1.\end{array}$$(23)

The q-triplet thus leads to a interesting mathematical structure. If we define ϵ ≡ 1 − q, the q-triplet becomes equivalent to the set: ${ϵ}_{\mathrm{stat}}=-\frac{3}{4}$, ϵ_rel = −3, and ${ϵ}_{\mathrm{sens}}=\frac{3}{2}$. These values satisfy the following relationships:

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{stat}}& =\frac{{ϵ}_{\mathrm{sens}}+{ϵ}_{\mathrm{rel}}}{2}(\mathrm{arithmetic}\mathrm{mean}),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{sens}}& ={\left({ϵ}_{\mathrm{stat}}{ϵ}_{\mathrm{rel}}\right)}^{1/2}(\mathrm{geometric}\mathrm{mean}),\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{rel}}^{-1}& =\frac{{ϵ}_{\mathrm{stat}}^{-1}+{ϵ}_{\mathrm{sens}}^{-1}}{2}(\mathrm{harmonic}\mathrm{mean}).\hfill \end{array}$$(26)

The interpretation of these intriguing relationships in terms of some underlying symmetry or analogous physical principle remains an open question (Gazeau & Tsallis 2019).

The aim of this work is to investigate and verify these relationships. To this end, we performed a systematic analysis of large-scale fluctuations in the solar wind proton density using data collected by several spacecraft located at the L1 point. Our study focuses on identifying multifractal structures, probability distributions, and relaxation processes. Subsequently, we analyzed the correlations among these three phenomena.

The novelty of this study lies in the fact that this is the first systematic investigation of the q-triplet in solar wind proton density, based on continuous data spanning 17 consecutive years. Previous works have already provided evidence supporting the q-triplet framework in astrophysical and atmospheric systems, but such studies have typically been restricted to specific years or conditions (see, e.g., Burlaga & F.-Viñas 2005b; Ferri et al. 2010). Were, we interpreted our results in the context of non-extensive statistical mechanics, which appears to be consistent with the observed nonlinear structure of the data.

3. Data analysis

Let us now consider some specific observations of the fluctuations of proton density in solar wind. The data we utilized in this study were taken from the OMNI directory¹ (King & Papitashvili 2005), which contains the hourly mean values of the interplanetary magnetic field and solar wind plasma parameters measured by various spacecraft near the Earth’s orbit. We used the non-shifted low resolution dataset, which is primarily a 1963-to-current compilation of hourly averaged, near-Earth solar wind magnetic field and plasma parameter data from several spacecraft in geocentric or L1 (Lagrange point) orbits. In particular, since 2004, the priority data are taken from two spacecrafts: Wind (Kasper 2002) and ACE (McComas et al. 1998). As an example, Fig. 1 shows observations of the hourly averages of proton density N_p in solar wind from day 1 to 365, year 2022.

Fig. 1. Time series, Np(t), of hourly averages of the proton density as a function of time, year 2022. The data were taken from the OMNI directory.

As can be seen, the fluctuations in N_p are large during this interval, that is, the amplitudes of the fluctuations are larger than the mean. For each year between 2008 and 2024, we wanted to deduce the parameters q_stat, q_rel, and q_sens.

3.1. Determination of q_stat

The value of q_stat is derived from a PDF. The successive fluctuations in N_p can be described by the PDFs of

$$\begin{array}{c}\hfill d{N}_p(i)\equiv N_p(i+1)-N_p(i),\end{array}$$(27)

properly normalized using the moving average $⟨N_p(i)⟩=\frac{N_p(i+1)+N_p(i)}{2}$. Our statistical analysis is based on the algorithm described in Ferri et al. (2010). The range of dN_p is subdivided into small “cells” (a data binning process) of width δN_p, in order to evaluate the frequency of dN_p values falling within each bin. The choice of bin width is a crucial step in the algorithmic process and is equivalent to solving the binning problem: a proper initialization of the bin size can significantly accelerate the statistical analysis and promote convergence of the algorithm toward the correct solution. In our case, we used the Sturges’ method.

The PDF observed for the year 2022 is shown in Fig. 2a as an example. The q-Gaussian distribution provides an excellent fit to all observed PDFs across the years 2008–2024.

Fig. 2. (a) PDFs of relative hourly changes in the proton density for the year 2022. The circles are the data points. The solid red curve is a nonlinear fit of the data with a q-Gaussian, and the dashed blue curve is a Gaussian distribution. (b) Linear correlation between lnq(p) and (dNp/⟨Np⟩)2. The red line is the best fit with qstat = 1.64 ± 0.01.

For an initial assessment, we performed a fast nonlinear fit of the PDF using a q-Gaussian (Eq. (8)) to obtain a preliminary estimate, q′. Since this method typically yields an error of around 20%, we reduced the uncertainty by linearizing the PDF. To do so, we considered the plot of ln_q(p) versus (dN_p/⟨N_p⟩)², as shown in Fig. 2b.

To refine the estimate, we varied q in steps of δq = 0.01 around the initial value q′, performing a linear regression at each step and calculating the corresponding correlation coefficient (CC). The value of q that yields the highest CC was selected as the best estimate of q_stat.

3.2. Determination of q_rel

To estimate q_rel, one can analyze the decay of specific observables Ω>(t), such as the autocorrelation function C(τ) or the mutual information I(τ). The value of q_rel can be determined from a autocorrelation coefficient C(τ), defined asfollows:

$$\begin{array}{c}\hfill C(\tau )\equiv \frac{⟨[N_p(t_i+\tau )-⟨N_p⟩]·[N_p(t_i)-⟨N_p⟩]⟩}{⟨{[N_p(t_i)-⟨N_p⟩]}^2⟩},\end{array}$$(28)

where N_p(t_i) represents the proton density at time t_i, and the average ⟨N_p⟩ denotes a time average over the full dataset. According to non-extensive statistics, it should decay as an asymptotic power law, i.e., log C(τ) = a + slog τ, where the slope s = 1/(1 − q_rel), and q_rel characterize a relaxation process. In Fig. 3 we show an example (year 2022) where the relaxation exhibits a power-law decay.

Fig. 3. Autocorrelation coefficient (C(τ)) vs. scale (τ) computed from hourly averages of proton density for the year 2022. The solid red line is the best fit to the first five data points with a q exponential (qrel) of 4.4 ± 0.5.

3.3. Determination of q_sens

The sensibility to initial conditions parameter, q_sens, can be derived from the multifractal spectrum f(α) of the attractor associated with the nonlinear dynamical system. The sensitivity to initial conditions in nonlinear systems is described by a q-exponential distribution with q = q_sens, rather than an exponential distribution, as is the case for strong chaos.

To investigate the presence of a multifractal structure in the time series, we plotted the moments of N_p at various timescales τ = 2ⁿ, n = 0, 1, 2, 3, …. For a given value of τ, we calculated the mobile averaged value ⟨N_p⟩ over the time interval τ. From this series, we constructed the moments, N_p^k, when k is any positive or negative number. In standard multifractal analysis, the notation q is used for these moments; however, we used k here to avoid confusion with the non-extensivity parameter, q.

To calculate the k-th order statistical moments of the proton density, N_p, we used a sliding-window (or “moving average”) technique over the normalized time series. Specifically, we first normalized the full dataset by dividing the raw proton density, N_p(t), by its global average, ⟨N_p⟩. Then we computed a moving average ⟨N_p⟩_τ(t) using windows of size τ days. The k-th moment at scale τ was then defined as

$$\begin{array}{c}\hfill {⟨N_p^k⟩}_{\tau }=\frac{1}{M}{\displaystyle \sum _{i=1}^M}{\left(\frac{1}{\tau }{\displaystyle \sum _{j=0}^{\tau -1}}\frac{N_p(t_i+j)}{⟨N_p⟩}\right)}^k,\end{array}$$(29)

where M is the number of moving windows over which the average is computed, and t_i marks the starting point of each window. This method captures how the distribution of the proton density evolves across multiple timescales and moment orders.

The result is a curve of the k-th moment of N_p as a function of scale. Finally, we repeated this procedure for multiple values of k, which yielded a family of curves, one for each value of k, as shown in Fig. 4. These curves are straight lines on a log-log plot, and the slope increases with the magnitude of k, indicating the presence of a multifractal structure over the analyzed range of scales.

Fig. 4. k-th moments of various mobile averages of Np as a function of scale for the year 2022. A range of scales is observed in which the points for a given moment (k) lie close to a straight line. Each colored line corresponds to a different k-order moment. The absolute value of the slope increases with increasing |k|, indicating the existence of multifractal structure.

This yields a set of slopes (k_i, s_i), which can be described by a nonlinear function s(k). In other words, if the proton density profile exhibits a multifractal structure, then

$$\begin{array}{c}\hfill ⟨N_p^k⟩\sim {\tau }^{s(k)}.\end{array}$$(30)

The function s(k) characterizes the specific multifractal structure. The set of observed points (k_i, s_i) can be approximated well by a polynomial function s(k), so just a few coefficients are sufficient to describe the multifractal, as shown in Fig. 5a. According to Mandelbrot (1972, 1989), s(k) is a quadratic polynomial when the time series follows a log-normal distribution. The variance of the log-normal distribution obeys a scaling symmetry (Gupta & Waymire 1991). In our case, the data deviate from the quadratic fit (see the solid blue line in Fig. 5a), which is expected since our distribution is not log-normal, but rather q-Gaussian. In practice, one fits the lowest-degree polynomial that provides a good fit to the data. In our example from the year 2022, we use a fourth-degree polynomial.

Fig. 5. (a) Points (ki, si) with error bars, derived from the slopes in Fig. 4. The solid curve is a fourth-degree polynomial fit. (b) Multifractal spectrum f(α) derived from the same data. The solid curve is a quadratic polynomial fit used to determine the zeros αmax and αmin.

It is useful to introduce two additional descriptions. The first is the “generalized dimension” (D_k(k); Hentschel & Procaccia 1983), which is related to s(k) by the equation

$$\begin{array}{c}\hfill D_k(k)=1+\frac{s(k)}{k-1}·\end{array}$$(31)

The D_k describes the Rényi generalized dimensions, defined as

$$\begin{array}{c}\hfill D_k=\frac{1}{k-1}\underset{\lambda \to 0}{lim}\frac{log{\displaystyle {\sum }_{i=1}^N{p}_i^k}}{log\lambda },\end{array}$$(32)

where p_i is the local probability at location i in phase space, and λ is the local scale. The Rényi k indices (typically denoted q, but we use k here to avoid confusion) can take values across the entire real line, ( − ∞, +∞).

The second description is given in terms of the multifractal spectrum (f(α); Halsey et al. 1986), defined by the relations

$$\begin{array}{cc}\hfill \alpha & =\frac{d}{\mathit{dk}}\left[(k-1)D_k(k)\right],\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill f(\alpha )& =k\alpha (k)-(k-1)D_k(k),\hfill \end{array}$$(34)

where α is the Hölder exponent. Using the coefficients of the fitted polynomial s(k), we calculated a set of points in the multifractal spectrum f(α) using Eqs. (33) and (34). The resulting spectrum (α, f(α)) is shown in Fig. 5b.

The extremes of the spectrum, α_min and α_max, for which f(α) = 0, are related to q_sens (Lyra & Tsallis 1998; Tsallis 2004) according to

$$\begin{array}{c}\hfill \frac{1}{1-q_{\mathrm{sens}}}=\frac{1}{{\alpha }_{\min }}-\frac{1}{{\alpha }_{\max }}·\end{array}$$(35)

To determine α_min and α_max, it is necessary to fit the observations with a suitable function and identify the intersection with the x axis, extrapolating f(α) if necessary. The uncertainties in α_min and α_max propagate to the uncertainty in q_sens, but these are largely influenced by the fitting function chosen. Although the theoretical form of f(α) is not known, it is expected to be a concave function with a single maximum (Beck & Schögl 1993). A quadratic function, shown by the curves in Fig. 5b, provides a good fit to our observations, although the fit is not unique. A cubic fit also performs well over the observed range, but its extrapolation yields an unphysical inflection point (Burlaga & F.-Viñas 2005b). For the year 2022, using a quadratic fit, we obtain a value of q_sens = −0.38 ± 0.02.

4. Results

After performing the analysis described in the last section for the 17 years under consideration (2008–2024), we present in Table 1 the results of the q-triplet for each year, and the average of all of them. The results support the idea that the theoretical conjectures and previous experimental findings are on the right track (Burlaga & F.-Viñas 2005b; Gazeau & Tsallis 2019). Note that the q-triplet values for each individual year do not coincide, within the error bars, with the theoretical predictions; however, the average values do. This reflects the dispersion in the values obtained, particularly in the determination of q_sens, a fact that is evident from the relatively large standard deviation.

Furthermore, in Fig. 6 we show a plot of the quantity a_rel + a_stat − a_sens as a function of the year, according to Eq. (23) should be equal to 1. The average over the range studied here is (1.0 ± 0.2).

Fig. 6. Value of the quantity arel + astat − asens vs. year, in blue. According to theory, the q-triplet holds the relationship arel + astat − asens = 1 for solar wind (dashed red line). Superimposed (dashed black line) is the annual mean F10.7 solar radio flux index, shown on the secondary y-axis as a proxy for solar activity. While a moderate anticorrelation is observed (Pearson correlation coefficient r = −0.47, p-value p = 0.055), the deviation is not statistically significant, suggesting a possible – but not definitive – influence of the solar cycle on the q-triplet structure.

It is important to recognize that the estimation of intermittency is subject to several uncertainties related to measurement quality, the length of the time series, and the spectral characteristics of the fluctuations, as previously noted by other authors (Sorriso-Valvo et al. 2017; Viall & Borovsky 2020). A hypothesis regarding the variability in the values of the q-triplet is its possible dependence on the solar cycle. We plot the F10.7 index together with the quantity α_rel + α_stat − α_sens in Fig. 6 for illustrative purposes. A linear correlation analysis shows that, while a moderate anticorrelation is observed (Pearson correlation coefficient r = −0.47, p-value p = 0.055), the result is not statistically significant, suggesting a possible, but not definitive, influence of the solar cycle on the q-triplet structure. For example, in Pavlos et al. (2015), the q-triplet was studied during both shock and calm periods in the solar wind, revealing different values for each regime. However, the cited work was based on high-resolution data, focusing on small-scale fluctuations, rather than investigating the long-term (large temporal scale) dependence of the q-triplet on solar activity. These differences can be physically understood as the result of differing dynamical conditions: shock periods are associated with abrupt changes in plasma parameters and enhanced turbulence, which intensify non-Gaussian features and long-range correlations. In contrast, calm periods tend to exhibit more regular, near-equilibrium behavior, with statistics closer to Gaussian and weaker correlations. Thus, the divergence in q-values reflects the underlying complexity and deviation from equilibrium in each regime. The study of therelationship between the q-triplet and solar activity remains an active topic of research, and we intend to present our findings on this topic in a forthcoming publication.

5. Conclusions

We have performed a comprehensive, year-by-year analysis of solar wind proton density fluctuations at the L1 point (near 1 AU), covering 17 consecutive years, from 2008 to 2024. Using the framework of non-extensive statistical mechanics, we examined the presence and behavior of three key features of nonlinear dynamical systems: fat-tailed probability distributions, long relaxation processes, and multifractal structures. These correspond, respectively, to the indices q_stat, q_rel, and q_sens of the Tsallis q-triplet.

Our results confirm both theoretical conjectures and earlier empirical studies (Burlaga & F.-Viñas 2005b; Gazeau & Tsallis 2019). Although the individual annual values of the q-triplet fluctuate and do not always match the theoretical expectations within their uncertainties, the average values over the full 17-year period do align with the predicted relationships among the indices. This agreement suggests that the Tsallis triplet structure is indeed a robust description of the solar wind’s complex behavior, and that the variability seen on a yearly basis reflects both measurement limitations and natural dynamical fluctuations.

The q-triplet has been validated against data obtained by astrophysical observations, such as those cited here, atmospherical observations (see, e.g., Ferri et al. 2010, 2017), and seismogenesis observations (Iliopoulos et al. 2012; Pavlos et al. 2014). A good summary of these findings is made in Pavlos et al. (2018). In none of these cases do we have control over the variables, and therefore the measurements are noisy. Future research should aim to obtain less noisy plasma and magnetic field measurements, particularly during transient events, to better capture the dynamics of nonequilibrium structures. Additionally, more systematic statistical analyses could include classifying solar wind intervals by physical regimes (e.g., fast vs. slow wind, shock-driven events, or coronal mass ejections), applying consistent criteria across extended datasets, and conducting time-localized or scale-dependent studies to distinguish persistent structures from transient fluctuations. It would also be beneficial to evaluate the q-triplet not only on a year-by-year basis but across broader and more physically meaningful periods, such as phases of high and low solar activity within the solar cycle. This would allow for a more comprehensive comparison of different measures, leading to a deeper understanding of the nature of the solar wind’s nonlinear character. This, in turn, would further strengthen the empirical support for the predictions of the q-triplet and other non-extensive theoretical frameworks. Another way to improve the results would be to search for experimental evidence in other phenomena in which variable control ispossible.

In particular, we find that the standard deviation is relatively large for q_sens, reflecting the intrinsic difficulty of estimating this index from multifractal spectra, which depend sensitively on the fitting method and extrapolation of f(α). Despite this, the relationship q_sens < 1 < q_stat < q_rel, previously noted in solar wind magnetic field studies, is preserved in our analysis of protondensity.

Our findings are especially significant in light of the differences in the observational context: while previous studies of the q-triplet specifically were based on interplanetary magnetic field measurements at heliocentric distances ranging from 7 to 87 AU, our study focuses on a plasma variable – proton density – measured continuously near Earth. The fact that the nonlinear character of the solar wind (as captured by the q-triplet) persists even at 1 AU highlights the relevance of these dynamics for understanding near-Earth space weather phenomena. We acknowledge that the broader field of turbulence and intermittency in space plasmas has also been extensively investigated using both magnetic and plasma field data from spacecraft such as MMS, Cluster, and Ulysses (see, e.g., Marino et al. 2008; Kiyani et al. 2009, and Yordanova et al. 2021). Although our focus in this work is on the q-triplet formalism, these contributions offer important context and support for understanding multi-scale dynamics in the solar wind. Our results suggest that long-range correlations, multifractal structure, and slow relaxation processes in the solar wind must be accounted for when modeling its interaction with the Earth’s magnetosphere and the resulting space weather effects. Such nonlinear features may influence critical technologies such as satellite navigation (e.g., GPS), communication systems, and power grids.

Future work could focus on applying the same analysis to other plasma parameters (such as velocity or temperature), investigating shorter temporal windows associated with specific events (e.g., coronal mass ejections), or exploring connections between the q-triplet and geomagnetic indices. Additionally, refined statistical techniques – such as ensemble grouping by solar wind regime or solar cycle phase – may help reduce dispersion in the estimated q-values and strengthen the predictive power of this framework.

Acknowledgments

We acknowledge G. Ferri for support with the analysis code, and thank GSFC/SPDF and OMNIWeb for providing access to the data. This work was financially supported by CONICET (Argentina). D.J.Z. is grateful to C. Tsallis for his valuable comments and corrections to the manuscript.

References

All Tables

All Figures

	Fig. 1. Time series, Np(t), of hourly averages of the proton density as a function of time, year 2022. The data were taken from the OMNI directory.
In the text

	Fig. 2. (a) PDFs of relative hourly changes in the proton density for the year 2022. The circles are the data points. The solid red curve is a nonlinear fit of the data with a q-Gaussian, and the dashed blue curve is a Gaussian distribution. (b) Linear correlation between lnq(p) and (dNp/⟨Np⟩)2. The red line is the best fit with qstat = 1.64 ± 0.01.
In the text

	Fig. 3. Autocorrelation coefficient (C(τ)) vs. scale (τ) computed from hourly averages of proton density for the year 2022. The solid red line is the best fit to the first five data points with a q exponential (qrel) of 4.4 ± 0.5.
In the text

	Fig. 4. k-th moments of various mobile averages of Np as a function of scale for the year 2022. A range of scales is observed in which the points for a given moment (k) lie close to a straight line. Each colored line corresponds to a different k-order moment. The absolute value of the slope increases with increasing |k|, indicating the existence of multifractal structure.
In the text

	Fig. 5. (a) Points (ki, si) with error bars, derived from the slopes in Fig. 4. The solid curve is a fourth-degree polynomial fit. (b) Multifractal spectrum f(α) derived from the same data. The solid curve is a quadratic polynomial fit used to determine the zeros αmax and αmin.
In the text

Fig. 6. Value of the quantity arel + astat − asens vs. year, in blue. According to theory, the q-triplet holds the relationship arel + astat − asens = 1 for solar wind (dashed red line). Superimposed (dashed black line) is the annual mean F10.7 solar radio flux index, shown on the secondary y-axis as a proxy for solar activity. While a moderate anticorrelation is observed (Pearson correlation coefficient r = −0.47, p-value p = 0.055), the deviation is not statistically significant, suggesting a possible – but not definitive – influence of the solar cycle on the q-triplet structure.

In the text