Nonextensive Footprints in Dissipative and Conservative Dynamical Systems

Abstract

Despite its centennial successes in describing physical systems at thermal equilibrium, Boltzmann–Gibbs (BG) statistical mechanics have exhibited, in the last several decades, several flaws in addressing out-of-equilibrium dynamics of many nonlinear complex systems. In such circumstances, it has been shown that an appropriate generalization of the BG theory, known as nonextensive statistical mechanics and based on nonadditive entropies, is able to satisfactorily handle wide classes of anomalous emerging features and violations of standard equilibrium prescriptions, such as ergodicity, mixing, breakdown of the symmetry of homogeneous occupancy of phase space, and related features. In the present study, we review various important results of nonextensive statistical mechanics for dissipative and conservative dynamical systems. In particular, we discuss applications to both discrete-time systems with a few degrees of freedom and continuous-time ones with many degrees of freedom, as well as to asymptotically scale-free networks and systems with diverse dimensionalities and ranges of interactions, of either classical or quantum nature.

1. Introduction

Statistical mechanics constitutes one of the pillars of contemporary theoretical physics. It was introduced in the 19th century by L. Boltzmann and J.W. Gibbs, and the name was coined by Gibbs himself. It is based on mechanics (classical, quantum, relativistic), electromagnetism, and theory of probabilities. Probabilities enter through the so-called entropic functional S, whose generic form for discrete stochastic variables is given by

$$S\left(\left\{p_i\right\}\right)=kF\left(\left\{p_i\right\}\right)\left(\sum _{i=1}^W{p}_i=1\right),$$

(1)

where $F\left(\left\{p_i\right\}\right)$ is an appropriate generic functional, k being typically equal either to unity or to the Boltzmann constant $k_B$. Historically, Boltzmann and Gibbs used continuous variables ($p\left(x\right)$ instead of $p_i$). The corresponding discrete form is given by

$$S_{BG}\left(\left\{p_i\right\}\right)=-k\sum _{i=1}^W{p}_{i}ln{p}_i.$$

(2)

Later on, J. von Neumann extended this functional by introducing quantum mechanical operators, and C.E. Shannon made an important connection to the theory of communications. The BG entropy (2) is additive in the following sense [1]: if we have two probabilistically independent systems A and B, i.e., $p_{ij}^{A+B}=p_i^A{p}_j^B,\forall (i,j)$, we straightforwardly verify that

$$S_{BG}(A+B)=S_{BG}\left(A\right)+S_{BG}\left(B\right).$$

(3)

The powerful BG statistical mechanics is grounded on entropy (2), and correctly describes a wide number of thermal properties of uncountable systems. However, it fails when generic long-range space–time correlations are present in the system, which is typically the case of the so called complex natural, artificial, and social systems (see [2] and references therein). A paradigmatic such situation is well known to occur in gravitation.

In order to overcome this crucial difficulty for vast classes of systems, the BG theory was generalized in 1988 [3] by building on the following generalized entropic form

$$S_q\left(\left\{p_i\right\}\right)=k\frac{1-{\sum }_{i=1}^W{p}_i^q}{q-1}(q\in \mathbb{R};S_1=S_{BG}),$$

(4)

hereafter referred to as the q-entropy. This entropic functional satisfies

$$\frac{S_q(A+B)}{k}=\frac{S_q\left(A\right)}{k}+\frac{S_q\left(B\right)}{k}+(1-q)\frac{S_q\left(A\right)}{k}\frac{S_q\left(B\right)}{k}.$$

(5)

It is, consequently, nonadditive for $q\ne 1$.

The entropy (4) can be alternatively written in the following forms:

$$S_q=k\sum _{i=1}^W{p}_i{ln}_q\frac{1}{p_i}=-k\sum _{i=1}^W{p}_i^q{ln}_q{p}_i=-k\sum _{i=1}^W{p}_i{ln}_{2-q}{p}_i,$$

(6)

where the q-logarithmic function is defined as follows:

$${ln}_{q}z\equiv \frac{z^{1-q}-1}{1-q}({ln}_{1}z=lnz).$$

(7)

Its inverse function, the q-exponential, is defined as follows,

$$e_q^z\equiv {[1+(1-q)z]}_{+}^{\frac{1}{1-q}}(e_1^z=e^z),$$

(8)

where ${[\cdots ]}_{+}=[\cdots ]$ if $[\cdots ]>0$, and zero otherwise. For the particular instance of equal probabilities, i.e., $p_i=1/W,\forall i$, we obtain

$$S_q=k{ln}_{q}W,$$

(9)

which recovers, for $q=1$, the celebrated formula $S_{BG}=klnW$, carved on Boltzmann’s tombstone in Vienna.

Under appropriate canonical constraints, the entropy (4) is extremized by the following distribution,

$$p_i=e_q^{-{\beta }_q(E_i-{\mu }_q)},$$

(10)

which, for $q=1$, recovers the ubiquitous BG weight

$$p_i=\frac{e^{-\beta E_i}}{Z}\equiv \frac{e^{-\beta E_i}}{{\sum }_{j=1}^W{e}^{-\beta E_j}},$$

(11)

Z being the usual BG partition function.

The real quantities ${\beta }_q\equiv 1/k{T}_q$ and ${\mu }_q$ play, respectively, the roles of q-generalized inverse temperature and chemical potential.

The BG theory is undoubtedly successful in a myriad of physical systems at thermal equilibrium. However, as mentioned above, it fails, or it is even ill-defined, for describing vast classes of natural, artificial and social complex systems at nontrivial stationary or quasistationary states. The goal of the present review is to focus on the most typical systems for which the BG theory fails, whereas the present q-generalized theory (also currently referred to as nonextensive statistical mechanics or q-statistics in the literature [4]) satisfactorily describes various relevant thermostatistical properties. Illustrations of this fact include distributions emerging in high-energy physics (e.g., at LHC/CERN and Brookhaven energies) [5,6], outer space [7], cold atoms in dissipative optical lattices [8,9], and granular matter [10], among many others.

In general, the set of probabilities $\left\{p_i\right\}$ can be directly defined within the realm of a purely probabilistic model (e.g., [11]), or through a dynamical model within a phase space, in such a way that the dynamics itself creates, along time, the necessary set of probabilities (e.g., [2] and references therein). The present paper mainly, though not exclusively, focuses on this latter possibility.

In Section 2 we discuss dynamical systems (either dissipative or conservative) with few degrees of freedom. In Section 3, we discuss dynamical systems (either dissipative or conservative) with many degrees of freedom and variable dimensionality and range of interaction. In Section 4, asymptotically scale-free networks are addressed. In Section 5, several clues about the domains of validity of either BG statistical mechanics or q-statistics are presented. Some final remarks and conclusions close the paper with Section 6.

2. Few Degrees of Freedom

In the study of the statistical characterization of nonlinear dynamical systems, models with few degrees of freedom are very popular due to the relative ease of use compared to more complex systems. Here, we will focus first on the most paradigmatic model for the dissipative case with few degrees of freedom, namely the logistic map [12,13,14,15], and we will review its statistical characterizations, which are available to date in the literature.

Then, we will do the same for the standard map [16], which can be considered as the most paradigmatic conservative (area-preserving) model studied in dynamical systems theory.

2.1. Dissipative Models

The logistic map is defined as

$$x_{t+1}=1-a{x}_t^2,$$

(12)

where $0<a\le 2$ is the control parameter and $x_t\in [-1,1]$, with $t=0,1,2,\cdots$. This map attains its first edge of chaos, denoted by $a_c$, at $a_c=1.40115518909\cdots$, which can be reached from the periodic region by period doubling through the accumulation of $2^{\infty }$ periods, or from the chaotic region by the merging of an infinite number of bands.

2.1.1. Sensitivity to Initial Conditions

We will briefly review the three different methods presented in the literature to study the sensitivity to initial conditions of the logistic map.

Sensitivity Function Method: The first method studies the sensitivity to initial conditions that use the sensitivity function $\xi \left(t\right)$, defined as

$$\xi \left(t\right)\equiv \underset{\Delta x\left(0\right)\to 0}{lim}\frac{\Delta x\left(t\right)}{\Delta x\left(0\right)},$$

(13)

where $\Delta x\left(0\right)$ and $\Delta x\left(t\right)$ are the discrepancies of the initial conditions at times 0 and t [17,18]. The functional form suggested for handling both the strongly chaotic region and the chaos threshold is

$$\xi \left(t\right)={\left[1+(1-q){\lambda }_{q}t\right]}^{1/(1-q)}(q\in \mathbb{R}),$$

(14)

where ${\lambda }_q$ is the generalized Lyapunov coefficient. It is easy to verify that the standard exponential form $exp\left({\lambda }_{1}t\right)$ is recovered for $q=1$ (here, ${\lambda }_1$ is the standard Lyapunov exponent). Generically, when $q\ne 1$ we have a power-law behaviour and if ${\lambda }_q<0$ and $q>1$ (${\lambda }_q>0$ and $q<1$), the system is said to be weakly insensitive (sensitive) to the initial conditions. This is clearly different from the standard case where we have strong insensitivity (sensitivity) for ${\lambda }_1<0$ (${\lambda }_1>0$). We calculate $\xi \left(t\right)$ through

$$ln\xi \left(t\right)=\sum _{t}ln\left[\frac{df\left(x_t\right)}{dx}\right],$$

(15)

with $f\left(x\right)=1-a{x}_t^2$ and $x_0=0$; $\xi \left(t\right)$ exhibits, at the chaos threshold, a power-law divergence in the form $\xi \propto t^{1/(1-q_{sens})}$, from where the value of $q_{sens}$ can be calculated by measuring, on a log–log plot, the upper bound slope $1/(1-q_{sens})$. We plot our results in Figure 1, from where the value $q_{sens}=0.245\cdots$ [18] is obtained. This constitutes the first method in the literature to evaluate $q_{sens}$ and it has already been used for several map families [18,19,20,21,22].

Multifractal Spectrum Method: The second method in the literature for estimating the value of $q_{sens}$ is based on the geometrical aspects of the attractor of the map at the chaos threshold. It is known [23] that the multifractal singularity spectrum $f\left(\alpha \right)$ reflects the fractal dimension of the subset with singularity strength $\alpha$. In this function, the two points where $f\left(\alpha \right)$ vanishes, namely ${\alpha }_{min}$ and ${\alpha }_{max}$, characterize the scaling behavior of the most concentrated and most rarefied regions on the attractor. By using these points, the scaling relation

$$\frac{1}{1-q_{sens}}=\frac{1}{{\alpha }_{min}}-\frac{1}{{\alpha }_{max}}(q_{sens}<1)$$

(16)

has been proposed, from which one can, once again, deduce the $q_{sens}$ value from a completely different method [19]. This relation can be arranged as $1/(1-q_{sens})=ln2/ln{\alpha }_F$, which clearly exhibits that this scaling is dependent on the Feigenbaum constant ${\alpha }_F$. Indeed, this enables us to estimate $q_{sens}$ value with an extraordinary precision because ${\alpha }_F$ is available in the literature with over 10,000-digit precision [24].

Entropy Increase Method: Finally, the third method proposed for the estimation of $q_{sens}$ is based on the generalization of the Kolmogorov–Sinai (KS) entropy [25,26,27]. Although the KS entropy is defined in principle through a single trajectory in phase space, an ensemble-based procedure is much simpler computationally. The KS entropy rate is defined through the limit $K_q\equiv {lim}_{t\to \infty }{lim}_{W\to \infty }{lim}_{N\to \infty }{S}_q\left(t\right)/t$, where t is the time, W is the number of regions in the partition of the phase space, N is the number of initial conditions (all chosen at $t=0$ within one region among the W available ones) that are evolving in time, and the entropy is given by

$$S_q\left(t\right)=\frac{1-{\sum }_{i=1}^W{\left[p_i\left(t\right)\right]}^q}{q-1}.$$

(17)

The method is based on the conjecture that there is a unique value $q_{sens}$, such that $K_q$ is finite for $q=q_{sens}$, vanishes for $q>q_{sens}$, and diverges for $q<q_{sens}$ [25]. This value coincides with the one previously obtained by the other two methods.

2.1.2. Relaxation Dynamics

Now let us focus on the dynamic evolution, in phase space, of an ensemble of initial conditions uniformly distributed over the phase space and explore its relaxation dynamics. Typically, a partition of the phase space on $N_{box}$ cells of equal size is performed and the evolution in time of a set of M identical copies of the system is implemented. Then, there is a proper entropy $S_q$, with $q=q_{rel}$, which evolves at a constant rate such that

$$K_{q_{rel}}=\underset{N_{box}\to \infty }{lim}[S_{q_{rel}}\left(t\right)-S_{q_{rel}}\left(0\right)]/t$$

(18)

goes to a constant value as $t\to \infty$. The number of occupied cells evolves in time as

$$W\left(t\right)={[W{\left(0\right)}^{1-q_{rel}}+(1-q_{rel})K_{q_{rel}}t]}^{1/(1-q_{rel})},$$

(19)

with the exponent $\mu =-1/(q_{rel}-1)<0$ governing the asymptotic power-law decay [28,29]. Although early numerical assumptions predicted a value of $\mu$ close to $-0.71$, more precise approaches estimated this value as $\mu =-0.800138\cdots$, from which $q_{rel}=2.249784\cdots$ is obtained [30,31,32]. Numerical results are shown in Figure 2. In order to reach the analytical result, we need to take the limits $N_{box}\to \infty$ and $M\to \infty$. However, our computer limitations allow us now to achieve $M=265\times {10}^6$, and this order of magnitude yields the approximative value $\mu =0.800\cdots$ for small time steps, provided that the ratio $r=M/N_{box}$ is large enough. In addition, due to the fact that the number of boxes is very small, the curve saturates. On the other hand, as we increase the number of boxes, because M stays constant due to computer limitation, this ratio decreases. This makes the initial part of the curve slightly depart from the theoretical result though the final part goes farther before saturation and stays closer to the analytic curve (black dotted line).

2.1.3. Central Limit Behavior

Now we can finally discuss a third q parameter, known as $q_{stat}$ (“stat” stands for stationary state), which comes from the analysis of the central limit behavior. In order to achieve this task, let us focus on the probability distribution of the sum of iterates of the logistic map (or any other dynamical system) given by

$$y:=\sum _{i=1}^N\left(x_i-\langle x\rangle \right),$$

(20)

where N is the number of iterations, $x_i$ are the iterates of the logistic map and $\langle \cdots \rangle$ stands for the time average. It can be analytically shown that, for strongly chaotic systems, the distribution becomes Gaussian for $N\to \infty$ [33,34] due to the standard central limit theorem (CLT). However, several complex systems exist where this tendency fails and, consequently, the standard CLT is violated. Defining

$$\langle x\rangle =\frac{1}{M}\frac{1}{N}\sum _{j=1}^M\sum _{i=1}^N{x}_i^{\left(j\right)},$$

(21)

(where M is the number of randomly chosen initial conditions) we can study the behavior of such systems. It should be noted here that the average is not only taken over a large time series but also over a large number of initial conditions. The first map studied in this way was the logistic map in the vicinity of chaos threshold [35,36]. It was shown that, if the chaos threshold is appropriately approached (see [37] for details) by using the Huberman–Rudnick law [38], it is possible to verify that the limit distribution seems to converge to a q-Gaussian defined as

$$P\left(y\right)=\left\{\begin{array}{cc}P\left(0\right){\left[1-\beta (1-q)y^2\right]}^{\frac{1}{1-q}}\hfill & for\beta (1-q)y^2<1\hfill \\ 0\hfill & otherwise\hfill \end{array}\right.,$$

(22)

where $q<3$ and $\beta >0$ are parameters. This distribution recovers the Gaussian case as $q\to 1$. Results for the probability distribution of y for the logistic map are given in Figure 3. It is clear that a q-Gaussian behavior gradually develops, with periodic oscillations, as we get closer to the chaos threshold. This is not only correct for the standard period-doubling route to chaos but also for any other periodic cycles available in the chaotic region [37].

2.2. Conservative Models

After all these findings for dissipative maps, it is natural to study conservative maps on similar grounds. A very suitable model is the standard map [16,39], given by

$$\begin{array}{ccc}\hfill p_{i+1}& =& p_i-Ksin{x}_i\hfill \\ \hfill x_{i+1}& =& x_i+p_{i+1},\hfill \end{array}$$

(23)

where $p_i$ and $x_i$ ($i=1,2,3\cdots$) are taken as modulo $2\pi$, and $K\ge 0$. This map is a low-dimensional area-preserving one, i.e., it is a conservative system like a Hamiltonian. Therefore, it is a natural system to be analysed in terms of a central limit behavior, exactly in the same sense that it has already been done for the dissipative case.

The phase space portrait of the standard map is shown in Figure 4 for four representative values of the map parameter. It is evident that the chaotic sea grows inside the stability islands as the value of K increases.

Because the chaotic sea dominates the full phase space for the $K=10$ case, one would expect that, for the probability distribution of the random variable y defined in Equation (20), the usual CLT holds. This is indeed verified in Figure 5. A Gaussian trend is clearly visible, as expected.

However, at the other end, where the full phase space is dominated by the stability islands (namely, for $K=0$ case), a completely non-Gaussian trend is evident from Figure 6. It was shown in [40] that the probability distribution can be approximated by a q-Gaussian with $q\simeq 1.94$. Later on, the asymptotic value of q was analytically found to be $q=2$ [41]. Obviously, this value can only be approached numerically as $N\to \infty$ and $M\to \infty$ limits are taken simultaneously.

For the values of the K parameter between these two extreme cases, the phase space is typically a mixture of the chaotic sea with periodic islands. For those cases, the behavior of the probability distribution has been found to be a linear combination of the Gaussian and the q-Gaussian cases that rule the chaotic and stability island regions of the map, respectively [40,42].

After these findings for the standard map, very recently other conservative maps have already been discussed in this perspective in the literature [43,44,45].

3. Many Degrees of Freedom

Long-range interacting systems with many degrees of freedom have been intensively studied in recent years and new methodologies have been developed in an attempt to understand their intriguing features. One of the most promising directions is the combination of statistical mechanics tools and methods adopted in dynamical systems [46]. In particular, phase transitions have been extensively explored in both conservative and dissipative long-range systems.

3.1. Coupled Pendula Models

In this context, the Hamiltonian Mean Field (HMF) model [47,48,49,50] and the Kuramoto model [51,52] represent two interesting toy models of coupled pendula, the former conservative and the latter dissipative, which are paradigmatic for many real systems with long-range forces and have several applications. Both models share the same order parameter and display a spontaneous phase transition from an homogeneous/incoherent phase to a magnetized/synchronized one.

The HMF model describes the dynamics of N classical spins or inertial rotors, characterized by the angles ${\theta }_i\in [-\pi ,\pi [$ and the conjugate momenta $p_i\in ]-\infty ,\infty [$, which can also be represented as particles moving on the unit circle. In its ferromagnetic version the Hamiltonian of the model is given by

$$H=K+V=\sum _{i=1}^N\frac{p_i^2}{2m}+\frac{1}{2N}\sum _{i,j=1}^N\left[1-cos({\theta }_i-{\theta }_j)\right],$$

(24)

where the mass m is usually set to 1. The potential term of Equation (24) reveals the mean field nature of the model, because each rotor can interact with all the others. Such a nature becomes more evident if we define as order parameter the magnetization $\mathbf{M}=M{e}^{i\varphi }=\frac{1}{N}{\sum }_{j=1}^N{e}^{i{\theta }_j}$, where $0<M<1$ and $\varphi$ are the modulus and the global phase. Within this assumption, the Hamilton equations of motion can be written as

$$\ddot{{\theta }_i}=\frac{1}{N}\sum _{j=1}^{N}sin({\theta }_j-{\theta }_i)=Msin(\varphi -{\theta }_i),i=1,\cdots ,N,$$

(25)

which correspond to the equations of single pendula in a mean field potential. In other words, each particle can be considered as moving in a mean-field potential determined by the instantaneous positions of all the other particles.

The Kuramoto model is considered one of the simplest models exhibiting spontaneous collective synchronization. It describes a large population of coupled limit-cycle oscillators, each one characterized by a phase ${\theta }_i$ and a natural frequency ${\omega }_i$, whose dynamics is given by

$${\dot{\theta }}_i={\omega }_i+\frac{K}{N}\sum _{j=1}^{N}sin({\theta }_j-{\theta }_i)i=1,\cdots ,N,$$

(26)

where $K\ge 0$ is the coupling strength. The natural frequencies are time-independent and are randomly chosen from a symmetric, unimodal distribution $g\left(\omega \right)$. Usually, uniform or Gaussian $g\left(\omega \right)$ distributions are considered. As in the case of HMF model, one can imagine the oscillators as particles moving on the unit circle. The order parameter of the Kuramoto model is perfectly equivalent to the magnetizaton in the HMF model and it is given by $\mathbf{r}=r{e}^{i\varphi }=\frac{1}{N}{\sum }_{j=1}^N{e}^{i{\theta }_j}$, where $\varphi$ is, again, the average global phase corresponding to the centroid of the phases of the oscillators and the modulus $0<r<1$ represents the degree of synchronization of the population. In terms of the variables r and $\varphi$, Equation (26) can be rewritten as

$${\dot{\theta }}_i={\omega }_i+Krsin(\varphi -{\theta }_i),i=1,\cdots ,N,$$

(27)

where the mean field character of the system becomes obvious, as was also the case for the HMF model.

As anticipated, HMF and Kuramoto models can be considered as limiting cases, respectively conservative and overdamped, of a more general model of coupled pendula [53,54]. Let us consider the following mean field equations describing a system of driven and damped pendula (with unit mass):

$$\ddot{{\theta }_i}+B\dot{{\theta }_i}+CMsin({\theta }_i-\varphi )=\Gamma ,i=1,\cdots ,N,$$

(28)

where the C is the coupling term, M the mean field order parameter, B the damping coefficient and $\Gamma$ the torque term.

In the conservative case, putting $B=\Gamma =0$ and $C=1$, from Equation (28) one recovers the mean field equations of motions of HMF model:

$$\ddot{{\theta }_i}=-Msin({\theta }_i-\varphi )i=1,\cdots ,N.$$

(29)

On the other hand, considering the overdamped limit of Equations (28), obtained for $B>>1$, the second derivatives can be neglected and, putting $C/B=K$, $M=r$ and $\Gamma /B={\omega }_i$, they reduce to the Kuramoto rate equations

$$\dot{{\theta }_i}={\omega }_i-Krsin({\theta }_i-\varphi )i=1,\cdots ,N,$$

(30)

where the constant natural frequencies play the role of the constant torque.

In the next two sections, we briefly summarize some interesting results obtained for both models, also showing that their common origin is actually reflected in some analogous dynamical features characterized by an anomalous behavior.

3.2. The Kuramoto Model

The main feature of the Kuramoto model, which makes it paradigmatic for a large class of phenomena [55,56,57,58,59], is its capacity to produce the spontaneous synchronization of a population of incoherent oscillators provided that a critical value $K_C$ of the coupling is exceeded. In fact, for small values of K, the oscillators act as if they were uncoupled and each oscillator tends to run independently and incoherently with its own frequency. Instead, for $K>K_C$, the coupling term in each of Equations (26) tends to synchronize each oscillator with all the others, and the system exhibits a spontaneous transition from the previous incoherent state to a synchronized one, where all the oscillators rotate at the same frequency $\Omega$ (a value which corresponds to the average frequency of the system, preserved by the dynamics). As shown by Kuramoto itself [51], the critical value of the coupling depends only on the central value $g(\omega =0)$ of the distribution $g\left(\omega \right)$ in accordance with the expression $K_C=\frac{2}{\pi g\left(0\right)}$.

For a given value of K, as the population becomes more coherent, the order parameter r grows and the effective coupling $Kr$ increases. In this regime of partial synchronization, as predicted by the solutions of Equation (27), two kinds of oscillators coexist depending on the size of $\left|{\omega }_i\right|$ relative to $Kr$: (i) oscillators with $\left|{\omega }_i\right|>Kr$, called $drifting-oscillators$, that run incoherently around the unit circle, and (ii) oscillators with $\left|{\omega }_i\right|\le Kr$, called $locked-oscillators$, which are trapped in a rotating cluster. The dynamic interplay between these two kinds of oscillators is probably at the root of the microscopic chaotic behavior which, as we will show, characterizes the regime of partial synchronization. On the other hand, when the effective coupling $Kr$ becomes strong enough, all the oscillators rotate in the same cluster at the frequency $\Omega$ and any fingerprints of chaos disappears: in fact, in this case, the system behaves like a single giant oscillator and becomes thus integrable.

It is interesting to study the shape of the phase transition in the Kuramoto model together with the behavior of the largest Lyapunov exponent (LLE), both as a function of the chosen distribution $g\left(\omega \right)$. Ref. [54] explored what happens by using bounded Gaussian distributions with $\omega \in [-3,3]$ and a standard deviation $\sigma$ going from $1.0$ to $5.0$, being $g\left(\omega \right)$ in this latter case almost uniform (see the original paper for details about the numerical integration and the adopted software). By calculating the order parameter and the LLE as function of K for several values of $\sigma$, one obtains the curves shown in Figure 7 for N = 20,000. In this figure, for increasing values of $\sigma$, a transformation from a continuous second-order phase transition (lower panel, on the right) toward an abrupt first-order-like one (on the left) is clearly visible, together with an analogous change in the peak of the largest Lyapunov exponent (upper panel), whose maximum value is also related with $g\left(\omega \right)$. In correspondence, the critical value $K_C$, initially depending on $\sigma$, moves from right to left toward a final value which does not depend on $\sigma$ any more and is very close to the theoretical value $K_C=12/\pi$ predicted for a true uniform $g\left(\omega \right)$ with $\omega \in [-3,3]$ (which in principle would strictly request $\sigma =\infty$). At the same time, the region of partial synchronization, which is mainly situated after the phase transition for small values of $\sigma$, progressively shifts before the phase transition for increasing values of $\sigma$.

As we will see in the next section, the Kuramoto second-order phase transition obtained for a Gaussian $g\left(\omega \right)$ is very similar to that observed for the HMF model, where an analogous behavior of the LLE has also been found. In the following, let us focus on a very peculiar dynamical feature, which is present in both models, even if with different characteristics.

In Figure 8, again for a system of N = 20,000 oscillators, we draw the temporal evolution of the order parameter $r\left(t\right)$ for several single runs for three values of $K/K_C$ near the phase transition, for both Gaussian (left column) and uniform (right column) $g\left(\omega \right)$. It clearly appears that in the latter case (panels (e) and (f)) stationary states with high and low asymptotic values of r coexist, whereas in the former case (panels (a), (b), and (c)) only partially synchronized stationary states are visible. Such a result reinforces the distinction between the first-order-like and second-order-like dynamical phase transitions, occurring in the Kuramoto model depending on the different $g\left(\omega \right)$ distributions, which seem to play a very crucial role. As already noticed in Ref. [54], in some cases (see for example panels (a) and (f)) metastable states appear for both $g\left(\omega \right)$ distributions. As we will see in the next section, this phenomenon can be also recognized in the HMF model, where analogous metastable quasistationary states (QSS) have been observed when the system starts from out-of-equilibrium initial conditions.

Another interesting feature, common to both Kuramoto and HMF models, is the violation of the standard CLT in conditions of weak chaos (i.e., when LLE $\sim 0$) (see Ref. [60]; see the original paper for details about the numerical integration and the adopted software). In the Kuramoto model, one can use the angles of the N oscillators as stochastic variables and build the corresponding rescaled sums y’s by picking out, for each oscillator, n values of the angle ${\theta }_i$ at fixed intervals of time $\delta$ along the deterministic time evolution imposed by Equation (26), i.e.:

$$y_i=\frac{1}{\sqrt{n}}\sum _{k=1}^n{\theta }_i\left(k\delta \right),\mathrm{with}i=1,2,\cdots ,N.$$

(31)

In such a way, the product $\delta \times n$ will give the total simulation time over which the sum is calculated. The resulting PDFs are shown in the bottom panels of Figure 9 for both uniform (left column) and Gaussian (right column) $g\left(\omega \right)$ and for values of the coupling below the critical threshold ($K=0.1$ and $K=0.6$, respectively).

Fat-tailed distributions clearly appear, which can be well fitted by q-Gaussians with $q=1.7$, as expected because both the order parameter (top panels) and the LLE (middle panels) stay very close to zero for the whole time interval considered. This confirms that, for these values of K, the system results to be incoherent and weakly chaotic regardless of the velocity distribution adopted. In the next section we will see that, again, a similar feature can be observed also in the context of the HMF model.

3.3. The HMF Model

A general Hamiltonian for the HMF model was introduced originally in Refs. [61,62] and then initially investigated from a dynamical point of view in Refs. [47,63,64]. The Hamiltonian can be written as

$$H=K+V=\sum _{i=1}^N\frac{p_i^2}{2}+\frac{\epsilon }{2N}\sum _{i,j=1}^N[1-cos({\theta }_i-{\theta }_j)].$$

(32)

As anticipated in Section 3.1, this model can be seen as N classical $XY$-spins with infinite range couplings, or also as a system of N particles moving on the unit circle. In the latter case the coordinate ${\theta }_i$ of particle i represents its position on the circle and $p_i$ its conjugate momentum. For $\epsilon >0$, particles attract each other or, in the other case, spins tend to align (ferromagnetic case, see Equation (24)), whereas for $\epsilon <0$, particles repel each other and spins tend to antialign (antiferrorromagnetic case) [47]. At short distances, particles cross each other or they collide elastically because they have the same mass, equal to 1.

As we have already seen, it is useful to introduce the mean field vector $\mathbf{M}=M{\mathrm{e}}^{i\varphi }=\frac{1}{N}{\sum }_{i=1}^N{\mathbf{m}}_i$, where ${\mathbf{m}}_i=(cos{\theta }_i,sin{\theta }_i)$. Here, M and $\varphi$ represent the modulus and the phase of the order parameter, which characterizes the degree of clustering in the particle interpretation, while it is the magnetization of the $XY$ spins. Thus, the potential energy can be rewritten as a sum of single-particle potentials $v_i$

$$V=\frac{1}{2}\sum _{i=1}^N{v}_i\mathrm{with}{v}_i=1-Mcos({\theta }_i-\varphi ).$$

(33)

The motion of each particle is, therefore, coupled to all the others, because the variables M and $\varphi$ are determined at each time t by the instantaneous positions of all particles.

The mean-field nature of the system can be also appreciated looking at Equation (25), presented in Section 3.1. These equations can be integrated numerically at fixed energy (see for example Refs. [47,62,63,64] for technical details). Considering the ferromagnetic case and defining the temperature T of the system as $2<K>/N$ (where $<K>$ denotes the time averaged kinetic energy), it is easy to verify that, starting from equilibrium initial conditions for both angles and velocities, one can recover the so-called caloric curve, i.e. the following relation between the temperature and the energy per particle $U=H/N$ predicted by the Boltzmann–Gibbs (BG) statistical mechanics in canonical ensemble [47,62],

$$U=\frac{1}{2\beta }+\frac{1}{2}(1-M^2),$$

(34)

where $\beta =1/T$ is the inverse temperature. The theoretical equilibrium solution for the caloric curve is reported as a solid line in panel (a) of Figure 10, and it is compared with the results of microcanonical numerical simulations for two different sizes of the system, namely $N={10}^2$ and $N={10}^3$. The agreement is clearly very good, and the same holds for the behavior of the magnetization M as function of U, shown in panel (b). In both plots, a second-order phase transition from a ferromagnetic (condensed) phase to a paramagnetic (homogeneous) one is also visible, in correspondence with the critical value of either the energy density, $U_c=0.75$, or the critical temperature, $T_c=0.5$, as predicted by the theory. On the other hand, if the system is started with strong out-of-equilibrium initial conditions in an energy range going from $U=0.5$ up to $U_c$, the model shows a nontrivial dynamical relaxation to the BG equilibrium. The class of out-of-equilibrium initial conditions we consider, called water bag initial conditions, consists of ${\theta }_i=0\forall i$ and the momenta uniformly distributed (according to the total energy density U). In Figure 11 we report, for $U=0.69$ and $N=500$, the time evolution of $2<K>/N$ (where $<K>$ denotes the time-averaged kinetic energy), a quantity that coincides with the temperature T. As expected, the system does not relax immediately to the BG equilibrium, but rapidly reaches a quasistationary state (QSS) consisting of a plateau associated with a temperature value $T_{\mathrm{QSS}}=2U-\epsilon$, lower than the canonical prediction, corresponding to the continuation of the linear part of the caloric curve to lower energies. The T vs. U plot (caloric curve), shown in inset (c) for the QSS regime, confirms a large disagreement with the equilibrium prediction (inset d). Furthermore, it turns out that the system remains trapped in such a state for a time that diverges with the size N of the system [65]. This means that if the thermodynamic limit is performed before the infinite-time limit, the QSS become stable and the system never relaxes to the BG equilibrium, exhibiting different equilibrium properties characterized by non-Gaussian velocity distributions (see inset a) [65]. Such velocity distributions have been fitted in Ref. [65] and have been shown to be in agreement with the prediction of the Tsallis’ generalized thermodynamics.

The characteristics of the QSS have been studied in several papers: vanishing Lyapunov spectrum [62] negative microcanonical specific heat [63], dynamical correlations in phase-space [49,50], Lévy walks and anomalous diffusion [64], nonergodicity [66], and glassy behaviour [50,67]. Let us here focus on the violation of the standard central limit theorem (CLT) due to the weak chaos observed in the QSS [66]. Analogous to what we saw in the previous section for the Kuramoto model, for the HMF model we can construct the probability density functions of quantities expressed as a finite sum of stochastic variables and select these variables along the deterministic time evolutions of the N rotors. More precisely, we can study the distribution quantity y defined as

$$y_i=\frac{1}{\sqrt{n}}\sum _{k=1}^n{p}_i\left(k\right)fori=1,2,\cdots ,N,$$

(35)

where $p_i\left(k\right)$, with $k=1,2,\cdots ,n$, are the rotational velocities of the ith rotor taken at fixed intervals of time $\delta$ along the same trajectory obtained by integrating the HMF equations of motions. The product $\delta \times n$ gives the total simulation time over which the sum of Equation (35) is calculated. In panel (a) of Figure 12 we plot the ensemble average PDF of the velocities calculated (over 100 different realizations) at $t=200$, i.e., at the beginning of the QSS regime, and after a very long transient, at t = 250,000 (full circles), for a system with N = 50,000 rotors. In panel (b), we plot the time average PDF for the normalized variable y with $n=5000$ and $\delta =100$, after a transient of 200 time units and over a simulation time of 500,000 units along the QSS. The shape of the time average PDF (b) results in a robust q-Gaussian, with $q=1.4$, not only in the tails, but also in the center (see inset). Because the time average PDF is completely different from the ensemble average PDF, we can also confirm the inequivalence between the two kinds of averages and the existence of a q-Gaussian attractor in the QSS regime of the HMF model [66].

3.4. Classical Inertial Rotors in d Dimensions

A generalization of the HMF (or XY) model, the so-called $\alpha$-XY model, was first introduced in Ref. [68]. Its Hamiltonian is

$$H_{\alpha }=K+V_{\alpha }=\sum _{i=1}^N\frac{p_i^2}{2}+\frac{\epsilon }{2\tilde{N}}\sum _{i\ne j}\frac{1-cos({\theta }_i-{\theta }_j)}{r_{ij}^{\alpha }},$$

(36)

where the parameter $\alpha$ determines the range of the interactions and $r_{ij}=|{\mathbf{r}}_i-{\mathbf{r}}_j|$ is the distance, in lattice units, between rotors i and j, which are located in a d-dimensional hypercubic lattice ($d=1,2,3$). The prefactor $1/\tilde{N}$, with

$$\tilde{N}=\sum _{i\ne j}\frac{1}{r_{ij}^{\alpha }}$$

(37)

in front of the potential energy term in (36) allows for an extensive total potential energy $V_{\alpha }$ for all values of $\alpha /d$ and reduces to $\tilde{N}=N-1\sim N$ for $\alpha =0$ (thus recovering the HMF model) and to $\tilde{N}=2d$ for $\alpha \to \infty$ (nearest-neighbours interactions).

Throughout the whole long-range regime, $0⩽\alpha /d⩽1$, the $\alpha$-XY model not only shares with the HMF model the appearance of QSSs, but also the same equation of state, given by Equation (34).

Finally, the equations of motion for the $\alpha$-XY model are given by

$$\left\{\begin{array}{cc}{\dot{\theta }}_i\hfill & =p_i\hfill \\ {\dot{p}}_i\hfill & =\epsilon {\mathbf{m}}_i\times {\mathbf{M}}_i\hfill \end{array},\right.$$

(38)

with ${\mathbf{M}}_i=\frac{1}{\tilde{N}}{\sum }_{j\ne i}\frac{{\mathbf{m}}_j}{r_{ij}^{\alpha }}$.

Equations (38) have been solved by means of a fourth-order symplectic algorithm [69], together with a fast Fourier transform algorithm, with an integration step $h=0.2$, which provided conservation of the energy per particle within a relative precision of ${10}^{-5}$ throughout all the calculations.

The scaling with the system size of the largest Lyapunov exponent, in the high-energy regime, as a function of the interaction range for the $\alpha$-XY model has been studied in Ref. [70]. A scaling in the form $\lambda \sim N^{-\kappa (\alpha /d)}$ was found, with a $\kappa$ dependence on $\alpha$ and d only through the quotient $\alpha /d$. The authors numerically obtained $\kappa \left(0\right)=1/3$, as analytically estimated in Ref. [71] for the HMF model and a decaying trend of $\kappa$ with $\alpha /d$ down to the value $\kappa (\alpha /d⩾1)=0$. Nevertheless, recent numerical results, to be published elsewhere, indicate that positive values of $\kappa$ are also obtained in the region $\alpha /d>1$.

In Ref. [72], probability distribution functions of time averaged momenta and individual energies for the $\alpha$-XY model were obtained, as shown in Figure 13 for the Boltzmann–Gibbs regime (that is, after the transition from the QSS to the final state characterized by the Boltmann–Gibbs temperature). Far from the usual Maxwellians and exponentials, $q_p-$Gaussian distributions were obtained for the averaged momenta and $q_E-$exponential distributions for the individual energies, with $q_p$ and $q_E$ depending on the quotient $\alpha /d$ as shown in Figure 14. We will come back to these results in Section 5.2.

The duration $t_{\mathrm{QSS}}$ of the quasistationary state of the $\alpha$-XY model as a function of N, $\alpha$ and d was studied in Ref. [73]. The relation $t_{\mathrm{QSS}}\sim N^{A(\alpha /d)}{e}^{-B{(\alpha /d)}^2}$ was obtained, in which the dependence of $t_{\mathrm{QSS}}$ on $\alpha$ and d is again through the quotient $\alpha /d$.

In addition, the duration of the QSS of the $\alpha$-XY model as a function of the total energy per particle U was studied in Ref. [74]. The QSS, which exists along the whole long-range interaction regime, provided $0.61<U<U_c=\frac{3}{4}$, disappears as $U\to U_c$ and this duration was found to go through a critical phenomenon, namely $t_{\mathrm{QSS}}\propto {(U_c-U)}^{-\xi }$, with a universal value $\xi \simeq 5/3$, independent of $\alpha /d$.

A further generalization consists of considering rotors of two degrees of freedom, instead of planar ones. The Hamiltonian of the so-called $\alpha$-Heisenberg model, first introduced in [75], is

$$H_{\alpha }=K+V_{\alpha }=\sum _{i=1}^N\frac{{\mathbf{L}}_i^2}{2}+\frac{\epsilon }{2\tilde{N}}\sum _{i\ne j}\frac{1-{\mathbf{S}}_i·{\mathbf{S}}_j}{r_{ij}^{\alpha }},$$

(39)

where now the three component spins ${\mathbf{S}}_i$, normalized to unity, with associated angular momenta ${\mathbf{L}}_i={\mathbf{S}}_i\times {\dot{\mathbf{S}}}_i$ (we consider unit moments of inertia), rotate freely in space and the order parameter is again the total magnetization, now given by

$$\mathbf{M}=\frac{1}{N}\sum _{i=1}^N{\mathbf{S}}_i.$$

(40)

The equation of state for the $\alpha$-Heisenberg model in the long-range regime is given by

$$U=T+\frac{\epsilon }{2}(1-M^2)$$

(41)

and the critical values of the temperature and the energy per particle for the para-ferro transition to occur are $(T_c,U_c)=(\frac{1}{3},\frac{5}{6})$.

Finally, the equations of motion for the $\alpha$-Heisenberg model are

$$\left\{\begin{array}{cc}{\dot{\mathbf{S}}}_i\hfill & ={\mathbf{L}}_i\times {\mathbf{S}}_i\hfill \\ {\dot{\mathbf{L}}}_i\hfill & =\epsilon {\mathbf{S}}_i\times {\mathbf{M}}_i\hfill \end{array},\right.$$

(42)

with ${\mathbf{M}}_i=\frac{1}{\tilde{N}}{\sum }_{j\ne i}\frac{{\mathbf{S}}_j}{r_{ij}^{\alpha }}$. Equations (42) have been integrated by means of the standard fourth order Runge–Kutta algorithm, together with a fast Fourier transform algorithm, with an integration step, $h=0.2$, such that the energy was conserved within a relative precision ${10}^{-5}$.

As in the case of the $\alpha$-XY model, $q_L-$Gaussian and $q_E-$exponential distributions were obtained for the time averaged angular momenta and energies, respectively, for the $\alpha$-Heisenberg model in [76], as shown in Figure 15. Also in this paper, a scaling of the largest Lyapunov exponent with the system size in the form $\lambda \sim N^{-\kappa (\alpha /d)}$ was obtained, as shown in Figure 16, with positive values of $\kappa$ beyond the long-range regime (that is, for $\alpha /d>1$). We will come back to these results in Section 5.2.

The duration $t_{\mathrm{QSS}}$ of the quasistationary state of the $\alpha$-Heisenberg model as a function of N, $\alpha$ and d has also been studied in Ref. [77]. As in the $\alpha$-XY model, a relation in the form $t_{\mathrm{QSS}}\sim N^{A(\alpha /d)}{e}^{-B{(\alpha /d)}^2}$ has also been obtained.

4. Asymptotically Scale-Free Networks

Rather unexpectedly a priori, some ubiquitous classes of growing networks—usually referred to as scale-free ones—are closely related [78,79,80,81,82,83,84,85,86] with various of the previous complex many-body systems. The relationship is neatly caused by the assumption of preferential attachment along the network growth. To illustrate this, we refer here to the d-dimensional model focused on in [85]. The attachment probability for the newly arrived site i is assumed to be given by

$$\begin{array}{c}\hfill {\Pi }_{ij}\propto \frac{{\epsilon }_i}{d_{ij}^{{\alpha }_A}}({\alpha }_A\ge 0),\end{array}$$

(43)

where $d_{ij}$ is the Euclidean distance between i and j, where j runs over all preexisting sites. The local energy ${\epsilon }_i$ is defined as

$$\begin{array}{c}\hfill {\epsilon }_i\equiv \sum _{j=1}^{k_i}\frac{w_{ij}}{2}(w_{ij}\ge 0),\end{array}$$

(44)

where the link weight $w_{ij}$ satisfies the probability distribution

$$\begin{array}{c}\hfill P\left(w\right)=\frac{\eta }{w_0\Gamma \left(\frac{1}{\eta }\right)}{e}^{-{(w/w_0)}^{\eta }}(w_0>0;\eta >0),\end{array}$$

(45)

which satisfies ${\int }_0^{\infty }dwP\left(w\right)=1$. As particular cases of Equation (45) we have $\eta =1$, which corresponds to an exponential distribution, $\eta =2$, which corresponds to a half-Gaussian distribution, and $\eta \to \infty$, which corresponds to an uniform distribution within $w\in [0,w_0]$. See Figure 17.

The numerical simulations strongly suggest, for all values of the model parameters, $p\left(\epsilon \right)\propto e_q^{-{\beta }_q\epsilon }$ where

$$\begin{array}{c}\hfill q=\left\{\begin{array}{cc}\frac{4}{3}\hfill & if0\le {\alpha }_A/d\le 1\hfill \\ \frac{1}{3}{e}^{1-{\alpha }_A/d}+1\hfill & if{\alpha }_A/d>1.\hfill \end{array}\right.\end{array}$$

(46)

See Figure 18. The long-range region corresponds to $0\le {\alpha }_A/d\le 1$, the intermediate region corresponds to $1<{\alpha }_A/d<\infty$, and the short-range region is attained only at the ${\alpha }_A/d\to \infty$ limit. If the attachment probability was assumed to decay exponentially with distance as ${\Pi }_{ij}\propto e^{-d_{ij}/\Delta }(\Delta >0)$, instead of as the power-law (43), the short-range region (i.e., $q=1$) would emerge for any $\Delta >0$.

5. Clues Concerning the Domains of Validity of BG and q-Statistics

In what follows, we present several suggestive clues that are available in the literature enabling us to conjecture the domains of validity and failure of the BG statistics and of q-statistics.

5.1. Clue I—Asymptotically Scale-Free Networks

In Section 4, we have focused on a d-dimensional network growth model which, by definition of its preferential attachment rule, is essentially scale-free [85]; its interaction decays as $1/r^{{\alpha }_A}({\alpha }_A\ge 0)$. We have obtained an energy distribution given by

$$p\left(\epsilon \right)\propto e_q^{-{\beta }_q\epsilon },$$

(47)

which is really not particularly surprising given the hypotheses under which the model has been constructed. The corresponding value of q is given by Equation (46). The value $q=4/3$ precisely corresponds to the divergence of its second moment $\langle {\epsilon }^2\rangle \equiv {\int }_0^{\infty }d\epsilon {\epsilon }^{2}p\left(\epsilon \right)$. Indeed, this moment is finite for all values of $q<4/3$ and diverges logarithmically precisely at $q=4/3$. For ${\alpha }_A/d>1$, the value of q starts decreasing but does not drop to $q=1$ until the limit $\alpha /d\to \infty$ is attained. This is because the interaction $1/r^{{\alpha }_A}$ cannot be properly considered as local unless $\alpha /d\to \infty$. This interaction is integrable for $\alpha /d>1$. This is necessary but not sufficient for the BG statistics to be justified. Integrability is equivalent to the finiteness of the zeroth moment of r, but the validity of the BG theory requires that all momenta be finite. In other words, ${\int }_1^{\infty }dr{r}^{d-1}{r}^n{r}^{-{\alpha }_A}$ is finite for all values of $n=0,1,2,\cdots$. This does not happen unless ${\alpha }_A\to \infty$. In contrast, if the interaction was given, for example, by $e^{-r/r_0}(r_0>0)$, or if it was nonvanishing only for, say, first, second, or third neighboring sites, all momenta would be finite, and we would naturally expect BG statistics to be verified. To summarize, the region $0\le {\alpha }_A/d\le$ can be properly called a long-range regime, but the region $1<{\alpha }_A/d<\infty$ cannot be properly be called a short-range regime. Indeed, it is a sort of intermediate-range regime. The proper short-range regime only emerges for ${\alpha }_A/d\to \infty$.

5.2. Clue II—Momenta and Energy Distributions of Classical Many-Body Hamiltonians

Let us first focus on the $\alpha$-XY inertial ferromagnet addressed in Section 3.4 (see Figure 13 and Figure 14).

For the momenta distributions, we verify that $q_p$ appears to be equal to 5/3 for $0\le \alpha /d\le 1$ and to possibly decrease exponentially down to $q_p=1$ for $\alpha /d$ increasing above unity. The second moment $\langle p^2\rangle$ of a $q_p$-Gaussian is finite for $q_p<5/3$ and diverges for $5/3\le q_p<3$.

For the energy distributions, we verify that $q_E$ appears to be equal to 4/3 for $0\le \alpha /d\le 1$ and to possibly decrease exponentially down to $q_E=1$ for $\alpha /d$ increasing above unity. The second moment $\langle E^2\rangle$ of a $q_E$-exponential is finite for $q_E<4/3$ and diverges for $4/3\le q_E<2$.

Let us focus now on the $\alpha$-Heisenberg inertial ferromagnet (see Figure 15). We see that $q_L$ appears to be compatible with $q_L=4/3$ within the interval $0\le \alpha /d\le 1$. Analogously, $q_E$ appears to be numerically compatible with $q_E=4/3$ within the same interval.

Let us finally focus on the $\beta$ version of the $\alpha$-Fermi–Pasta–Ulam–Tsingou d-dimensional oscillators (see Figure 19). In this case, the numerical $\alpha$-dependence of the index $q_p$ of the $q_p$-Gaussian distribution of velocities follows a different path. Indeed, it appears to vanish, like for the two previous Hamiltonian models ($\alpha$-XY and $\alpha$-Heisenberg), neatly above the value $\alpha =1$, but it is not constant within the interval $0\le \alpha \le 1$. However, these preliminary results are relatively old and not totally consistent. Therefore, we should be cautious when revisiting these simulations with higher precision.

5.3. Clue III—Maximal Lyapunov Exponent of the Classical $\alpha$-Heisenberg Inertial Ferromagnet

This property has been focused on in [76] for the $d=1,2,3$ models. From the available data, it has been possible to plot Figure 16.

Clues I, II, and III enable, for classical many-body Hamiltonian systems, the conjecture indicated in Figure 20 for each of the $q_p$ and $q_E$ indices for the momenta $q_p$-Gaussians and the energy $q_E$-exponential distributions, respectively.

Another, more generic, conjecture is possible as well, as indicated in Figure 22. This alternative conjecture recovers that of Figure 20 as ${(\alpha /d)}_c\to \infty$ is its particular case.

5.4. Clue IV—Viscous-Fluid Spherical Capacitor

A spherical capacitor whose interior is constituted by a dissipative fluid containing a large number of equally charged small masses has recently been focused on in [90] within molecular dynamics. The interaction between the particles is a repulsive Coulombian one. The radial density of the stationary state exhibits a q-exponentially decaying profile with $q\simeq 1.35$ rather than the commonly assumed exponentially decaying profile within the Debye–Huckel theory for electrolites (or the Yukawa theory in nuclear physics) (see Figure 23).

The repulsive Coulombian interaction between the particles corresponds to $\alpha =1$. Therefore, given that $d=3$, we have $\alpha /d=1/3$, which definitively corresponds to long-range interactions. This result is numerically consistent with the value $q=4/3$ suggested in Figure 20.

5.5. Clue V—Overdamped Many-Body Systems

Several types of many-body fluid confined d-dimensional systems with repulsive two-body interactions between particles whose motion is overdamped are susceptible to analytical study (see [91,92], to mention but a few).

If the interaction decays exponentially with distance (without exclusion of the origin) we typically obtain, for the $d=1$ stationary state distribution for the spatial distribution associated with a parabolic confinement, a q-Gaussian with $q=0$. The same distribution is verified for the scaled velocity distribution.

If the repulsive interaction decays as the power-law $\propto 1/r^{\alpha }$ (with exclusion of the origin, and $\alpha >d$), the stationary-state d-dimensional spatial distribution associated with a parabolic confinement is once again proven to be a q-Gaussian with $q=1-\alpha /d$. The same class of distribution is approximately observed for the scaled velocity distribution.

Both classes of systems are indicated in Figure 24.

The clue that is provided by the latter of these two examples is that $q\ne 1$, and not $q=1$, within the intermediate region $\alpha /d>1$.

5.6. Clue VI—Kinetics of Point Defects in Short-Range-Interacting Hamiltonians

The kinetics of point defects in short-range interacting d-dimensional Hamiltonians with the $O\left(n\right)$ symmetry for $n=d$ has been focused in [93,94]. It was found that the distribution of the point defect velocities is a q-Gaussian with

$$q_p=\frac{d+4}{d+2}.$$

(48)

It happens that this value of $q_p$ is the upper bound for a d-dimensional $q_p$-Gaussian to have a finite second moment. For $d=1$ we obtain $q_p=5/3$, which coincides with the value conjectured in Figure 20 for the plateau associated with $\alpha /d\in [0,1]$.

5.7. The Intriguing Case Of The Lennard–Jones’ Two-Body Potentials For Modeling Real Gases

Mie introduced in 1903 [95] a model for real gases ($d=3$) with the two-body potential

$$V\left(r\right)=Cϵ\left[{\left(\frac{\sigma }{r}\right)}^n-{\left(\frac{\sigma }{r}\right)}^m\right](0<m<n;ϵ>0)$$

(49)

with

$$C\equiv \frac{n}{n-m}{\left(\frac{n}{m}\right)}^{\frac{m}{n-m}}.$$

(50)

Two decades later, Lennard–Jones studied [96,97,98] the particular case $(n,m)=(12,6)$ on physical and thermodynamical grounds. The Lennard–Jones potential corresponds therefore to

$$V_{LJ}\left(r\right)=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],$$

(51)

where the attractive term is interpreted as a van der Waals interaction, and the repulsive term is a phenomenological one representing a hard, strong interaction. For nearly one century, this model for a real fluid has been handled on BG grounds with qualitatively successful comparisons with experiments. However, it corresponds to $\alpha /d=2$, which, according to the conjecture indicated in Figure 20, is included within the intermediate-range region, for which we expect $q\ne 1$. How are we supposed to handle this puzzling situation? A possible outcome might be argued as follows. A relation of the type

$$q_E(\alpha /d)\simeq 1+\kappa (\alpha /d)$$

(52)

has some degree of plausibility. Indeed, for, say, the $\alpha =0$ XY ferromagnet (MFA model), we might have $q_E\simeq 4/3$ and $\kappa \simeq 1/3$. In that case, the $\alpha /d=2$ case would lead (see caption of Figure 16) to a value similar to $q\simeq 1.001$. This value surely is numerically indistinguishable from $q=1$ at the precision within which comparisons of empirical and theoretical (Lennard–Jones) results are performed for real gases.

6. Final Remarks and Conclusions

The basic footprints of Boltzmann–Gibbs statistical mechanics are the energy distribution, the velocity distribution, and the time and size dependencies of the entropy. They are, respectively, given by exponentials (the BG weight, as well as the time sensitivity to the initial conditions, and phenomena relaxing with time), Gaussians (the Maxwellian distribution of velocities, as well as similar fluctuating quantities), and $q_{entropy}=1$ (i.e., $S_{BG}$). Within the q-generalized theory that has been focused on here, these footprints become q-exponentials, q-Gaussians, and q-entropies, respectively. These are precisely the properties that have been reviewed above, for systems with both low and large number of degrees of freedom, either conservative or dissipative, with dimensionalities $d=0,1,2,3$, diverse ranges of interactions (long, intermediate, short), involving both discrete and continuous time, with either classical or quantum nature, among other specificities. It has been profusely exhibited, in diverse complex dynamical systems, the various situations in which the q-concepts satisfactorily replace the BG concepts that are intimately related to ergodicity, mixing, symmetry of homogeneous occupancy of phase space, and similar ones.

Naturally, several interesting sets of open questions arise, which remain challenges for the future. One of them deserves a special mention at this stage, namely the energy and velocity distributions of classical many-body Hamiltonian d-dimensional systems with nonfrustrating two-body interactions of arbitrary range. The coupling constant is assumed to decay with the two-body distance r as $1/r^{\alpha }$: $0\le \alpha /d\le 1$ corresponds to long-range interactions, $1<\alpha /d<\infty$ corresponds to intermediate-range, and, finally, $\alpha /d\to \infty$ corresponds to short-range interactions, precisely the case which is most currently studied within the BG theory. First-principle dynamical numerical studies appear to indicate that, for $0\le \alpha /d<\infty$, two physically different regions might exist in the $(N,t)\to (\infty ,\infty )$ limit, one of them corresponding to $N/t^{\gamma }$ below some $(\alpha /d)$-depending threshold (with $\gamma >0$), and the other one corresponding to $N/t^{\gamma }$ above that same threshold. The correct description in the former case appears to be BG statistics, whereas in the latter case it appears to be q-statistics. In the region of validity of q-statistics, the energy distribution would be of the $q_E$-exponential form, with $q_E=4/3$ for $0\le \alpha /d\le 1$, and would exponentially decay down to $q_E=1$ for $\alpha /d$ increasing above unity up to infinity. The momenta distribution could be of the $q_m$-Gaussian form with $(q_m-1)/(q_E-1)=2$ for all values of $\alpha /d$. At the $\alpha /d\to \infty$ limit, the correct description follows BG statistics for any $(N,t)\to (\infty ,\infty )$ ordering. In other words, in the $\alpha /d\to \infty$ limit, the abovementioned threshold diverges.

These and related remarks suggest the following features.

-

The spatially averaged two-body potential ${\int }_1^{\infty }dr{r}^{d-1}{r}^{-\alpha }$ is finite for $\alpha /d>1$, and diverges for $0\le \alpha /d\le 1$. Such finiteness is necessary but not sufficient for all the BG thermostatistical quantities to be finite. Consistently, the total internal energy is thermodynamically extensive for $\alpha /d>1$, and superextensive for $0\le \alpha /d\le 1$.

-

The finiteness of the spatially averaged two-body potential is necessary for BG statistical mechanics to be applicable but it is not sufficient. Its full applicability requires also the finiteness of all the associated momenta, i.e., ${\int }_1^{\infty }dr{r}^{d-1}{r}^{-\alpha }{r}^n$ must also be finite for $n=1,2,3,\cdots$. Such a strong requirement is satisfied only in the $\alpha /d\to \infty$ limit of the present power-law models, or for Hamiltonians involving interactions only among relatively close neighbors (first, second, and third neighbors, for instance).

-

The maximal Lyapunov exponent appears to decay with the number N of elements as $N^{-\kappa (\alpha /d)}$ with $\kappa (\alpha /d)\ge 0$. It is possible that roughly $\kappa (\alpha /d)\simeq q_E(\alpha /d)-1$, for all values of $\alpha /d$. If so, we can guarantee strong chaos (hence, mixing in phase-space, hence ergodicity) in the $N\to \infty$ limit only for $\alpha /d\to \infty$. In all other cases, i.e., $0\le \alpha /d<\infty$, we would have, in the $N\to \infty$, weak chaos, and therefore ergodicity and mixing will not be guaranteed. This is consistent with the failure of the BG theory which is observed (nonexponential energy distribution, and non-Gaussian momenta distribution).

-

The fact that a BG partition function, as well as other thermostatistically relevant quantities (e.g., equations of states, energy and velocity distributions) are computable (within analytical mean-field methods, for example) is necessary but not sufficient for the BG theory to satisfactorily describe the system.

Further insights and analytical results in the statistical mechanics of complex systems are certainly very welcome in order to corroborate the above scenario involving some conjectural issues.