{'fraction_non_alphanumeric': 0.04532930781492468, 'fraction_numerical': 0.01238269167836671, 'mean_word_length': 4.613296989808012, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonextensive statistics: Theoretical, experimental and computational evidences and connections

The domain of validity of standard thermodynamics and Boltzmann-Gibbs statistical mechanics is discussed and then formally enlarged in order to hopefully cover a variety of anomalous systems. The generalization concerns {\it nonextensive} systems, where nonextensivity is understood in the thermodynamical sense. This generalization was first proposed in 1988 inspired by the probabilistic description of multifractal geometries, and has been intensively studied during this decade. In the present effort, after introducing some historical background, we briefly describe the formalism, and then exhibit the present status in what concerns theoretical, experimental and computational evidences and connections, as well as some perspectives for the future. In addition to these, here and there we point out various (possibly) relevant questions, whose answer would certainly clarify our current understanding of the foundations of statistical mechanics and its thermodynamical implications


I. INTRODUCTION
A diffuse belief exists, among many physicists as well as other scientists, that Boltzmann-Gibbs (BG) statistical mechanics and standard thermodynamics are eternal and universal. It is certainly fair to say that "eternal", in precisely the same sense that Newtonian mechanics is "eternal", they indeed are. But, again in complete analogy with Newtonian mechanics, we can by no means consider them as universal. Indeed, we all know that, when the involved velocities approach that of light, Newtonian mechanics becomes only an approximation (an increasingly bad one!) and reality is better described by special relativity. Analogously, when the involved masses are as small as say the electron mass, once again Newtonian mechanics becomes but a (bad) approximation, and quantum mechanics becomes necessary to understand nature. Also, if the involved masses are very large, Newtonian mechanics has to be extended into general relativity. In these senses we certainly cannot consider Newtonian mechanics as being universal. I believe that the same type of considerations apply to standard statistical mechanics and thermodynamics. Indeed, after more than one century highly successful applications of the magnificent Boltzmann's connection of Clausius macroscopic entropy to the theory of probabilities applied to the microscopic world , BG thermal statistics can (and should!) easily be considered as one of the pillars of modern science. However, it is unavoidable to think that, like all other products of human mind, this formalism must have physical restrictions, i.e., domains of applicability, out of which it can at best be but an approximation. It seems that BG statistics satisfactorily describe nature if the effective microscopic interactions are short-ranged (i.e., close spatial connections) and the effective microscopic memory is short-ranged (i.e., close time connections) and the boundary conditions are non(multi)fractal. Roughly speaking, the standard formalisms are applicable whenever (and probably only whenever) the relevant space-time (hence the relevant phase space) is non(multi)fractal. If this is not the case, some kind of extension appears to become necessary. Indeed, an everyday increasing list of physical anomalies are, here and there, being pointed out which defy (not to say that plainly violate!) the standard BG prescriptions. A nonextensive thermostatistics, which recovers the extensive, BG one as particular case, was proposed in 1988 [1,2] which might correctly cover at least some of the known anomalies. Although a fair amount of what legitimately looks like being successful applications is nowadays accumulating, further verifications and deep understanding is needed and welcome. Computational work is highly desired since, on various grounds, the analytic discussion frankly appears to be untractable. Needless, of course, to say that more experimental and theoretical work is absolutely relevant to exhibit the applicability and robustness of the ideas I intend to present herein. In the present contribution, I propose some (hopefully relevant!) questions that are right now open to such theoretical, experimental and computational contributions.
Let us be more specific. As mentioned above, it is nowadays quite well known that a variety of physical systems exists for which the powerful (and beautiful!) BG statistical mechanics and standard thermodynamics present serious difficulties or anomalies, which can occasionally achieve the status of just plain failures. Within a long list that will be systematically focused on later on, we may mention at this point systems involving longrange interactions (e.g., d = 3 gravitation) [3], long-range microscopic memory (e.g., nonmarkovian stochastic processes, on which much remains to be known, in fact) [4,5], and, generally speaking, conservative (e.g., Hamiltonian) or dissipative systems which in one way or another involve a relevant space-time (hence, a relevant phase space) which has a (multi)fractal-like structure. For instance, pure-electron plasma two-dimensional turbulence [6], Lévy anomalous diffusion [7], granular systems [8], phonon-electron anomalous thermalization in ion-bombarded solids ( [9] and references therein), solar neutrinos [10], peculiar velocities of galaxies [11], inverse bremsstrahlung in plasma [12] and black holes [13], to cite a few, clearly appear to be (in some cases), or could possibly be (in others), concrete examples. The present status of these and others will be discussed in Sections III, IV and V.
II. FORMALISM 2 A. ENTROPY As an attempt to overcome at least some of these difficulties a proposal has been advanced, one decade ago [1], (see also [14,15]), which is based on a generalized entropic form, namely where k is a positive constant and W is the total number of microscopic possibilities of the system (for the q < 0 case, care must be taken to exclude all those possibilities whose probability is not strictly positive, otherwise S q would diverge; such care is not necessary for q > 0; due to this property, the entropy is said to be expansible for q > 0). This expression recovers the usual BG entropy (−k W i=1 p i ln p i ) in the limit q → 1. The entropic index q (intimately related to and determined by the microscopic dynamics, as we shall mention later on) characterizes the degree of nonextensivity reflected in the following pseudo-additivity entropy rule where A and B are two independent systems in the sense that the probabilities of A + B factorize into those of A and of B (i.e., p ij (A + B) = p i (A)p j (B)). We immediately see that, since in all cases S q ≥ 0 (nonnegativity property), q < 1, q = 1 and q > 1 respectively correspond to superadditivity (superextensivity), additivity (extensivity) and subadditivity (subextensivity). Eq.
(2) exhibits a property which has apparently never been focused before, and which we shall from now on refer to as the composability property. It concerns the nontrivial fact that the entropy S(A + B) of a system composed of two independent subsystems A and B can be calculated from the entropies S(A) and S(B) of the subsystems, without any need of microscopic knowledge about A and B, other than the knowledge of some generic universality class, herein the nonextensive universality class, represented by the entropic index q, i.e., without any knowledge about the microscopic possibilities of A and B nor their associated probabilities. This property is so obvious for the BG entropic form that the (false!) idea that all entropic forms automatically satisfy it could easily install itself in the mind of most physicists. To show counterexamples, it is enough to check that the recently introduced Anteneodo-Plastino's [16] and Curado's [17] entropic forms satisfy a variety of interesting properties, and nevertheless are not composable. The above pseudo-extensivity property can be equivalently written as follows: We come back onto this form later on in connection with Renyi's entropy. Another important (since it eloquently exhibits the surprising effects of nonextensivity) property is the following. Suppose that the set of W possibilities is arbitrarily separated into two subsets having respectively W L and W M possibilities (W L + W M = W ). We define p L ≡ W L i=1 p i and p M ≡ W i=W L +1 p i , hence p L + p M = 1. It can then be straightforwardly established that where the sets {p i /p L } and {p i /p M } are the conditional probabilities. This would precisely be the famous Shannon's property were it not for the fact that, in front of the entropies associated with the conditional probabilities, appear p q L and p q M instead of p L and p M . This fact will play, as we shall see later on, a central role in the whole generalization of thermostatistics. Indeed, since the probabilities {p i } are generically numbers between zero and unity, p q i > p i for q < 1 and p q i < p i for q > 1, hence q < 1 and q > 1 will respectively privilegiate the rare and the frequent events. This simple property lies at the heart of the whole proposal. Santos has recently shown [18], strictly following along the lines of Shannon himself, that, if we assume (i) continuity (in the {p i }) of the entropy, (ii) increasing monotonicity of the entropy as a function of W in the case of equiprobability, (iii) property (2), and (iv) property (4), then only one entropic form exists, namely that given in definition (1). Of course, the generalization of Eq. (4) to the case where, instead of two, we have R nonintersecting subsets (W 1 + W 2 + ... + W R = W ) is straightforward [19]. To be more specific, if we define π j ≡ W j terms p i (j = 1, 2, ..., R) (hence R j=1 π j = 1), Eq. (4) is generalized into where we notice, in the last term, the emergence of what we shall soon introduce generically as the unnormalized q-expectation value (of the conditional entropies S q ({p i /π j }), in the present case). Another interesting property is the following. The Boltzmann-Gibbs entropy S 1 satisfies the following relation: Moreover, Jackson introduced in 1909 [20] the following generalized differential operator (applied to an arbitrary function f (x)): which satisfies D 1 ≡ lim q→1 D q = d dx . Abe [21] recently remarked that This property provides some insight into the generalized entropic form S q . Indeed, the inspiration for its use in order to generalize the usual thermal statistics came [1] from multifractals, and its applications concern, in one way or another, systems which exhibit scale invariance. Therefore, its connection with Jackson's differential operator appears to be rather natural. Indeed, this operator "tests" the function f (x) under dilatation of x, in contrast to the usual derivative, which "tests" it under translation of x.
Another property which no doubt must be mentioned in the present introduction is that S q is consistent with Laplace's maximum ignorance principle, i.e., it is extremum at equiprobability (p i = 1/W ∀i). This extremum is given by which, in the limit q → 1, reproduces Boltzmann's celebrated formula S = k ln W (carved on his marble grave in the Central Cemetery of Vienna). In the limit W → ∞, S q diverges if q ≤ 1, and saturates at k/(q − 1) if q > 1.
Finally, let us close the present set of properties by reminding that S q has, with regard to {p i }, a definite concavity for all values of q (S q is always concave for q > 0 and always convex for q < 0). In this sense, it contrasts with Renyi's entropy , which does not have this property for all values of q .
Before addressing other relevant quantities, let us introduce the following convenient functions [22]: (hence, e x 1 = e x ) with the definition supplement, for q < 1, that e x q = 0 if 1 + (1 − q) x ≤ 0, (and analogously, for q > 1, e x q diverges at x = 1/(q − 1)) and (hence, ln 1 x = ln x). We can verify easily that For instance, Eq. (10) can be rewritten in the following Boltzmann-like form: Let us also introduce the following unnormalized q-expectation value: hence A 1 corresponds to the standard mean value of a physical quantity A.
If our system is a generic quantum one, its probabilistic description is given by the density operator ρ, whose eigenvalues are the {p i }. Then, the generalized entropy is given by and the unnormalized q-expectation value of an observable A which does not necessarily commute with ρ is given by Eq. (16) can be rewritten as follows: and also as If our system is a generic classical one, the relevant variables are typically continuous variables, and its probabilistic description is given by a distribution of probabilities p(r), where r is a dimensionless variable in say a many-body phase space. Then, the generalized entropy is given by and the unnormalized q-expectation value of an observable A(r) is given by Although we shall, in what follows, be illustrating the present formalism with the case of W discrete microscopic possibilities, the generic quantum and classical discussions follow along the same lines, mutatis mutandis.

B. CANONICAL ENSEMBLE
Once we have a generalized entropic form, say that given in Eq. (1) (or an even more general one, or a different one), we can use it in a variety of ways. For instance, if we are interested in information theory, some optimization algorithms, image processing, among others, we can take advantage of a particular form in different ways. See, for instance, [17,19,23,24] and references therein, where it can be verified that not less than 25 (!) different entropic forms have received, along the years, a great variety of technological and mathematical applications. For instance, the Renyi entropy mentioned above has been quite useful in the geometrical characterization of strange attractors and similar multifractal structures (see [25] and references therein).
However, if our primary interest is Physics, this is to say the (qualitative and quantitative) description and possible understanding of phenomena occurring in Nature, then we are naturally led to use the available generalized entropy in order to generalize statistical mechanics itself and, if unavoidable, even thermodynamics. It is along this line that we shall proceed from now on (see also [26]). To do so, the first nontrivial (and quite ubiquitous) physical situation is that in which a given system is in contact with a thermostat at temperature T . To study this, we shall follow along Gibbs' path and focus the so called canonical ensemble. More precisely, to obtain the thermal equilibrium distribution associated with a conservative (Hamiltonian) physical system in contact with the thermostat we shall extremize S q under appropriate constraints. These constraints are [15] and where {ǫ i } are the eigenvalues of the Hamiltonian of the system. We shall refer to ... q as the normalized q-expectation value and to U q as the generalized internal energy (assumed finite and fixed). It is clear that, in the q → 1 limit, these quantities recover the standard mean value and internal energy respectively. We immediately verify that, for any observable, The outcome of this optimization procedure is given by It can be shown that, for the case q < 1, the expression of the equilibrium distribution is complemented by the auxiliary condition that p i = 0 whenever the argument of the function becomes negative (cut-off condition). Also, it can be shown [15] that Furthermore, it is important to notice that, if we add a constant ǫ 0 to all {ǫ i }, we have (as it can be self-consistently proved) that U q becomes U q +ǫ 0 , which leaves invariant the differences {ǫ i − U q }, which, in turn, (self-consistently) leaves invariant the set of probabilities {p i }, hence all the thermostatistical quantities. It is also trivial to show that, for the independent systems A and B mentioned previously, U q (A + B) = U q (A) + U q (B), thus recovering the same form of the standard (q = 1) thermodynamics. It can be shown that the following relations hold: and where Let us now make an important remark. If we take out as factors, in both numerator and denominator of Eq. (25), the quantity 1 + (1 − q)βU q / W j=1 (p j ) q , and then cancel them, we obtain where β ′ is an increasing function of β [27].
Let us now address the all important question of the connection between experimental numbers (those provided by measurements), and the quantities that appear in the theory. The definition of U q suggests the following normalized q-expectation values where O is any observable which commutes with the Hamiltonian, hence with ρ. If it does not commute, Eq. (34) is generalized into Consistently, O q is the mathematical object to be identified with the numerical value provided by the experimental measure. Later on, we come back onto this crucial point. At this point let us make some observations about the set of escort probabilities from which follows the dual relation The W = 2 illustration of P (q) i is shown in Fig. 1. As anticipated, q < 1 (q > 1) privileges the rare (frequent) events.
Let us first comment that Eqs. (36) and (37) have, within the present formalism, a role somehow analogous to the direct and inverse Lorentz transformations in Special Relativity (see [29] and references therein). Second, we notice that O q becomes an usual mean value when expressed in terms of the probabilities {P The final equilibrium distribution reads If the energy spectrum {ǫ i } is associated with the set of degeneracies {g i }, then the above probability leads to the following one (associated with the level ǫ i and not the state i) If the energy spectrum {ǫ i } is so dense that can practically be considered as a continuum, then the discrete degeneracies yield the function density of states g(ǫ), hence The density of states is of course to be calculated for every specific Hamiltonian (given the boundary conditions). For instance, for a d-dimensional ideal gas of particles or quasiparticles, it is given [30] by g(ǫ) ∝ ǫ d δ −1 , where δ is the exponent characterizing the energy spectrum ǫ ∝ K δ where K is the wavevector (e.g., δ = 1 corresponds to the harmonic oscillator, δ = 2 corresponds to a nonrelativistic particle in an infinitely high square well, etc). In Figs. 2 and 3 we see typical energy distributions for the particular case of a constant density of states. Of course, the q = 1 case reproduces the celebrated Boltzmann factor. Notice the cut-off for q < 1 and the long algebraic tail for q > 1.
All the above considerations refer, strictly speaking, to thermodynamic equilibrium. The word thermodynamic makes allusion to "very large" (N → ∞, where N is the number of microscopic particles of the physical system). The word equilibrium makes allusion to asymptotically large times (t → ∞ limit) (assuming a stationary state is eventually achieved). The question arises: which of them first? Indeed, although both possibilities clearly deserve the denomination "thermodynamic equiulibrium", nonuniform convergences might be involved in such a way that lim N →∞ lim t→∞ could differ from lim t→∞ lim N →∞ . To illustrate this situation, let us imagine a classical Hamiltonian system including two-body interactions decaying at long distances as 1/r α in a d-dimensional space, with α ≥ 0. If α > d the interactions are essentially short-ranged, the two limits just mentioned are basically interchangeable, and the prescriptions of standard statistical mechanics and thermodynamics are valid, thus yielding finite values for all the physically relevant quantities. In particular, the Boltzmann factor certainly describes reality, as very well known. But, if 0 ≤ α ≤ d, nonextensivity is expected to emerge, the order of the above limits becomes important because of nonuniform convergences, and the situation is certainly expected to be more subtle. More precisely, a crossover (between q = 1 and q = 1 behaviors) is expected to occur at t = τ (N). If lim N →∞ τ (N) = ∞, then we would indeed have two (or even more) different and equally legitimate states of thermodynamic equilibrium, instead of the familiar unique state. The conjecture is illustrated in Fig. 4.
A remark is necessary with regard to both sets just mentioned. Indeed, these properties have in fact been studied, whenever applicable, within unnormalized q-expectation values for the constraints, rather than within the normalized ones that we are using herein. Nevertheless, they still hold because they have been established for fixed β, which, through Eq. (33), implies fixed β ′ .

III. THEORETICAL EVIDENCES AND CONNECTIONS
A. Levy-type anomalous diffusion An enormous amount of phenomena in Nature follow the Gaussian distribution: measurement error distributions, height and weight distributions in biological individuals of given species, Brownian motion of particles in fluids, Maxwell-Boltzmann distribution of particle velocities in a variety of systems, noise distribution in uncountable electronic devices, energy fluctuations at thermal equilibrium of many systems, to only mention a few. Why is it so? Or, equivalently, what is their (thermo)statistical foundation? This fundamental problem has already been addressed, particularly by Montroll, and satisfactorily answered (see [7] and references therein). The answer basically relies onto two pillars, namely the BG entropy and the standard central limit theorem. However, the Gaussian is not the only ubiquitous distribution: we also similarly observe Levy distributions (in micelles [69], supercooled laser [70], fluid motion [71], wandering albatrosses [72], heart beating [73], financial data [74][75][76], among many others). So, once again, what is the (thermo)statistical foundation of their ubiquity? This relevant question has also been addressed, once again by Montroll and collaborators [7] among others. In this case however, a satisfactory answer has been missing for a long time. The first successful step toward (what we believe to be) the solution was performed in 1994 by Alemany and Zanette [77], who showed that the generalized entropic form S q was able to provide a power-law (instead of the exponential-law associated with Gaussians) decrease at long distances. Many other works followed along the same lines [78,79]. In [79] it was exhibited how the Levy-Gnedenko central limit theorem also plays a crucial role by transforming, through successive iterations of the jumps, the power-law obtained from optimization of S q into the specific power-law appearing in Levy distributions. Summarizing, in complete analogy with the above mentioned Gaussian case (and which is recovered in the more powerful present formalism as the q = 1 particular case), the answer once again relies onto two pillars, which now are the generalized entropy S q and the Levy-Gnedenko central limit theorem.
The arguments have been very recently re-worked [80] out on the basis of the normalized q-expectation values introduced in [15]. These are the results that we briefly recall here.
Let us write S q as follows: where x is the distance of one jump, and σ > 0 is the characteristic length of the problem. We optimize (maximize if q > 0, and minimize if q < 0) S q with the norm constraint ∞ −∞ dx p(x) = 1, as well as with the constraint We straightforwardly obtain the following one-jump distribution. If q > 1: If q = 1: If q < 1: We see that the support of p q (x) is compact if q ∈ (−∞, 1), an exponential behavior is obtained if q = 1, and a power-law tail is obtained if q > 1 (with p q (x) ∝ (σ/x) 2/(q−1) in the limit |x|/σ → ∞). Also, we can check that is finite if q < 5/3 and diverges if 5/3 ≤ q ≤ 3 (the norm constraint cannot be satisfied if q ≥ 3). Finally, let us mention that the Gaussian (q = 1) solution is recovered in both limits q → 1 + 0 and q → 1 − 0 by using the q > 1 and the q < 1 solutions respectively. This family of solutions is illustrated in Fig. 5.
We focus now the N-jump distribution p q (x, N) = p q (x) * p q (x) * ... * p q (x) (N-folded convolution product). If q < 5/3, the standard central limit theorem applies, hence, in the limit N → ∞, we have i.e., the attractor in the distribution space is a Gaussian, consequently we have normal diffusion. If, however, q > 5/3, then what applies is the Levy-Gnedenko central limit theorem, hence, in the limit N → ∞, we have where L γ is the Levy distribution with index γ < 2 given by Through the Fourier transforms of both Eq. (48) and (49), we can characterize the width and These results are depicted in Fig. 6. This result should be measurable in specifically devised experiments. More details can be found in [80] and references therein. What we wish to retain in this short review is that the present formalism is capable of (thermo)statistically founding, in an unified and simple manner, both Gaussian and Levy behaviors, very ubiquitous in Nature (respectively associated with normal diffusion and a certain type of anomalous superdiffusion).

B. Correlated-type anomalous diffusion
There are some phenomena exhibiting anomalous (super and sub) diffusion of a type which differs from the one discussed in the previous subsection. We refer to the so called correlated-type of diffusion. We consider here a quite large class of them, namely those associated with the following generalized, Fokker-Planck-like equation: where (µ, ν) ∈ R 2 , D is a dimensionless diffusion-like constant, F (x) ≡ −dV /dx is a dimensionless external force (drift) associated with a potential V (x), and (x, t) is a dimensionless 1 + 1 space-time. If µ = 1, we can interpret p(x, t) as a probability distribution since dx p(x, t) = 1, ∀t can be satisfied. If µ = 1, then p(x, t) must be seen as a density function. The word "correlated" is frequently used in this context due to the fact , an effective diffusion emerges, for ν = 1, which depends on p(x, t) itself, a feature which is natural in the presence of correlations. The µ = 1 particular case of this nonlinear equation is commonly denominated "Porous medium equation", and corresponds to a variety of physical situations (see [46] and references therein for several examples). The first connection of Eq. (53) with the present nonextensive statistical mechanics was established in 1995 by Plastino and Plastino [45]. They considered a particular case, namely µ = 1 and F (x) = −k 2 x with k 2 > 0 (so called Uhlenbeck-Ornstein processes), and found an exact solution which has the form of Eq. (43)(44)(45). Their work was generalized in [46] where arbitrary µ and F (x) = k 1 − k 2 x were considered. The explicit exact solution of Eq. (53), for all values of (x, t), was once again found by proposing an Ansatz of the form of Eqs. (45)(46)(47), i.e., the form which optimizes S q with the associated simple constraints. This form eventually turns out to be the Barenblatt one, useful in related problems. Here, let us restrict ourselves to just reproduce the exact solution of Eq. (53) assuming that p(x, 0) = δ(x), this is to say, a Dirac delta distribution. We obtain [46] where and and Summarizing, by using the form which optimizes S q , it has been possible to find the physically relevant solution of a nonlinear equation in partial derivatives with integer derivatives. It can be shown [81] that the problem that was solved in the previous subsection corresponds to a linear equation in partial derivatives but with fractional derivatives. We believe that we are allowed to say that an unusual mathematical versatility has been observed, within the present nonextensive formalism, in this couple of nontrivial examples of anomalous diffusion.

C. Stellar polytropes and other self-gravitating systems
The present formalism has been applied to a variety of astrophysical [42,43] and cosmological [82] self-gravitating systems. In some sense, this is something natural to do given the long range of the gravitational interaction. This was, in fact, the first physical application of nonextensive statistics. We do not intend here to reproduce details. Our present aim is to remind that it is well known in astrophysics that, within the standard thermodynamical approaches, it is not possible to simultaneously have finite values for the total energy, entropy and mass of a self-gravitating system. Plastino and Plastino were the first to show, in 1993, that this physically desirable situation can be achieved if we allow q to sufficiently differ from unity ! In fact, it can be shown (by considering the Vlasov equation in D-dimensional Schuster spheres) that the problem becomes a mathematically well posed one if q < q * , where the critical value q * is given by .
For D = 3 we recover the 7/9 relatively known value. Also, we notice that D = 2 implies q * = 1, which is very satisfactory since it is known that D < 2 gravitation is tractable within standard thermodynamics.

D. Zipf-Mandelbrot law
The problem we focus here first appeared in Linguistics. However, its relevance is quite broad, as it will soon become clear. Suppose we take a given text, say Cervante's Don Quijote, and order all of its words from the most to the less frequent; we refer to the ordered position of a given word as its rank R (low rank means high frequency ω of appearance in the text, and high rank means low frequency). Zipf [83] discovered that, in this as well as in a variety of similar problems, the following law is satisfied: where A > 0 and ξ > 0 are constants. Later on, Mandelbrot [74] suggested that such behavior was reflecting a kind of fractality hidden in the problem; moreover, he suggested how the Zipf law could be numerically improved: This expression has been useful in a variety of analysis. The connection we wish to mention here is that in 1997 Denisov [84] showed that, by extending (to arbitrary q) the well known Sinai-Bowen-Ruelle thermodynamical formalism of symbolic dynamics (i.e., by considering S q instead of S 1 ), the Zipf-Mandelbrot law can be deduced. He obtained and D = d/(q − 1) (d being a positive constant), i.e., Clearly, to make the discussion complete, a model would be welcome, which would provide quantities such as q and d. Nevertheless, Denisov's arguments have the deep interest of explicitly exhibiting that the Zipf-Mandelbrot law can be seen as having a nonextensive foundation. Fittings with experimental data will be shown later on in connection with the citations of scientific papers.

E. Theory of financial decisions; Risk aversion
An important problem in the theory of financial decisions is how to take into account extremely relevant phenomena such as the risk aversion human beings (hence financial operators) quite frequently feel. This kind of problem has, since long, been extensively studied by Tversky [85] and co-workers. The situation can be illustrated as follows. What do you prefer, to earn 85,000 dollars or to play a game in which you have 0.15 probability of earning nothing and 0.85 probability of earning 100,000 dollars ? In fact, most people prefer take the money. The problem of course is the fact that the expectation value for the gain is one and the same (more precisely 85,000 dollars) for both choices, and therefore this mathematical tool does not reflect reality ! The same problem appears if one expects to loose 85,000 and the chance is given for playing a game in which, if you win, you pay nothing, but, if you loose, you pay 100,000 dollars. In this case, most people choose to play. So, the experimental facts are that most human beings are risk-averse when they expect to gain, and risk-seeking when they expect to loose ! The problem is how to put this into mathematics.
Following [86], let us introduce, for the above gain problem, normalized q-expectation values as follows: gain take the money q = 85, 000 (65) and gain play the game Since most people would prefer the money, this means that most people have q < 1 for this particular decision problem. For the loss probem we have: and gain play the game Since in this case most people would prefer to play, this means that, consistently with the previous result, most people have q < 1 for the particular decision problem we are considering now. In some sense, we have some epistemological progress ! Indeed, the statement "most people have (for this type of amount of money) q < 1", unifies the previous two separate statements concerning expectation to gain and expectation to loose. Let us address now the following question: how can we measure the value of q associated with a particular individual ? We illustrate this interesting point with the example of the gain. The person is asked to choose between having V dollars or playing a game in which, if the person wins, the prize will be 100, 000 dollars and, if the person looses, he (she) will receive nothing. As before, the person is informed that his (her) probability of winning is 0.85 (hence, the probability of loosing is 0.15). Then we keep gradually changing the value V and asking what is the preference. At a certain critical value, noted V c , the person will change his (her) mind. Then, the value of q to be associated with that person, for that problem, is given by the following equality (See Fig. 7). The ideally rational operator corresponds to q = 1. For this gain problem, the risk-averse operators correspond to q < 1, and the risk-seeking ones to q > 1. It is clear that models for stock exchange can be formulated by using these remarks. Such an effort is presently in progress [87].

F. Physiology of vision
Physiological perceptions such as the visual perception are since long known to focus upon rare events (e.g., a red spot on a white wall). Barlow [88], among others, has recurrently stressed our attention on the fact that, at the action decision level, the various possibilities should enter with a weight proportional to − ln p i , and not proportional to p i , p i being the a priori probability of occurrence of that particular event; indeed, − ln p i diverges when p i → 0. He even argues that evolutionary arguments hold very well together with such hypothesis. To privilege rare events is precisely what happens, in the present formalism, whenever q < 1. Let us be more specific: if we consider the 0 < q << 1 limit, we obtain [89] and where O is an arbitrary observable. Leaving aside several constant quantities that appear above, we immediately observe the prominent role which − ln p i plays in these expressions. Consistently, the q → 0 limit of the present formalism could well be of some utility in the theoretical analysis of the physiological phenomena focused here.

IV. EXPERIMENTAL EVIDENCES AND CONNECTIONS
A. D = 2 turbulence in pure-electron plasma A few years ago, in 1994, Huang and Driscoll [6] exhibited some quite interesting nonneutral plasma experimental results obtained in pure-electron plasma (confined in a 20 cm long and 6 cm wide metallic cylindric Penning trap with a 10 −10 torr vacuum in its interior) in the presence of an external axial magnetic field (507 Gauss). In the interval 2-100 ms after every single electric shot (generating the electron plasma), it was observed a turbulent axisymmetric metaequilibrium state, the electronic density radial distribution of which was measured. Its average (over typically 100 shots) monotonically decreased with the radial distance, disappearing at some radius sensibly smaller than the radius of the container (i.e., a cut-off was observed). The experiment was recently redone [90] under slightly modified experimental conditions (a slow external rotation was imposed in such a way as to compensate the small energy dissipation existing in the plasma), and essentially the same metaequilibrium state was observed during lapses of time as long as 27 hours, or even longer ! In addition to the 1994 experiment, the authors also proposed [6] a phenomenological theory trying to reproduce the experimentally observed profile. Their proposal consisted on the optimization, for a given model, of a functional of the electron density ρ(r) under constraints, namely conservation of total mass, angular momentum and energy. They presented four different attempts. The first one (Point Vortex Maximum Entropy) consisted in optimizing, for a point vortex representation of the plasma, the BG entropy: it failed in reproducing the experimental data. The second attempt (Fluid Maximum Entropy) was essentially the same as the previous one, but using a fluid model for the plasma: the failure was even bigger. They assumed next that the problem possibly relied, not so much in the particular plasma model, but rather in the chosen functional to be optimized. In their third attempt (Global Minimum Enstrophy), they turned back to the point vortex model, but optimized the enstrophy instead of the BG entropy. The result was better than the two first attempts, but had the physically unacceptable feature of producing a negative electron density at sufficiently high radius. They then addressed their fourth attempt (Restricted Minimum Enstrophy), whose only difference with the third one was the fact of introducing an out-of-the pocket cut-off of the electron density at the proper value of the radius. This procedure was, finally, successful, and a very good first-approximation fitting was obtained ! The effort done by Huang and Driscoll was, on top of the high merit of a remarkable experiment, extremely pedagogical and elucidating: the main theoretical problem was not the model, but rather the choice of the functional to be optimized, i.e., the statistics.
The next important step in this story was done by Boghosian. He realized in 1995 and published [43] in 1996 that the Huang and Driscoll fourth, successful attempt precisely corresponds to the optimization of S q with q = 1/2 ! Indeed, by following the recipes of the present generalized thermostatistics, he re-obtained, for the electron density profile, the same differential equation produced within the Restricted Minimum Enstrophy phenomenological theory, with the supplementary bonus of not having to introduce in an ad hoc manner the necessary cut-off. Indeed, as already argued, all q < 1 cases exhibit a cut-off intrinsic to the formalism, and the radial position of that cut-off nicely fits the experimental value. The next step was performed in 1997 by Anteneodo and myself [91] (in fact, after related remarks by Boghosian himself). The Restricted Minimum Enstrophy theory is based on the enstrophy functional, which belongs to the general discussion of Casimir invariants; its form is in fact that of the order 2 Casimir invariant. Consequently, an epistemologically conservative theoretical viewpoint is to appreciate Boghosian's effort as just a formal interesting remark, with no real physical necessity. It happens, however, that, for r → r c − 0, (r c ≡ cut-off radius) the enstrophy theory yields ρ(r) ∝ (r c − r) whereas the experimental data fit much better a vanishing derivative at r c ! We followed [91] along Boghosian's lines and generalized his theory for arbitrary q. We obtained the generalized differential equation for ρ(r) and showed that ρ(r) ∝ (r c − r) q/(1−q) . Consequently, the experimental data fit better for q slightly above 1/2. This, together with the numerical solution of the differential equation, advanced q ≃ 0.55 as a better value for satisfactory overall fitting. (Better fittings would probably demand for a model more sophisticated than the point vortex one used here). The conceptually important point of this discussion is that Casimir invariants are characterized by integer exponents (in ρ(r)), hence none of them can be related to a value of q close to 0.55. From this standpoint, the present formalism appears as the only satisfactory phenomenological theoretical approach available in the literature at the present time.
The last step of this analysis was performed very recently by Anteneodo [92]. Indeed, the calculations above recalled [43,91] were done by using unnormalized q-expectation values. However, as already mentioned and used in the present review, it has been recently argued [15] that normalized q-expectation values should be used instead. It is therefore important to check that the present discussion and results for turbulence remain essentially invariant. This is now done [92], and it is this theory that we present in what follows.
The generalized entropy and associated constraints are given by gd 2 r = 1 (mass conservation) r 2 g q d 2 r g q d 2 r = L q ≡ L (angular momentum conservation) where g(r) is the probability distribution. Moreover, the scaled electrostatic potential satisfies 19 The constrained optimization of S q [g] (δ(S q − α gd 2 r − λL q − βU q )) now yields or, taking the Laplacian of both sides, which can be rewritten as where q . This equation precisely is the one appearing in [91], which, for q = 1/2, recovers that of [43]. For any chosen q, the values of the parameters (A, B) are obtained by imposing the experimental values of total angular momentum and energy. This phenomenological theory has, therefore, only one fitting parameter (q). As said before, q = 1/2 exactly reproduces the Huang and Driscoll's Restricted Minimum Enstrophy profile. The best overall fitting is, however, obtained for a value of q slightly above 1/2.

B. Solar neutrino problem
As easily conceivable, the core of the Sun is a very complex and turbulent plasma, within which an enormous amount of nuclear reactions take place. Many of them constitute chains of nuclear reactions in which neutrinos are produced. For instance, the p-p chain is described in [93]. Through a quite complete analysis of the production of neutrinos within the so called Solar Standard Model (SSM), it is possible to predict the neutrino flux onto the Earth. However, the actual flux measured in a variety of underground laboratories (Gallex, Sage, Kamiokande, Super-Kamiokande, Homestake) roughly amounts to only half of the predicted value ! This problem is currently referred to as the "solar neutrino problem". Two nonexclusive sources of explanation of this enigmatic discrepancy are: (i) the possible neutrino oscillations, which would make that only part of the predicted value would be detectable on the Earth; (ii) the current use of the SSM might be incorrect because it uses BG thermal statistics, which could be inappropriate for the solar plasma. Clayton [10] was the first to address the second possibility, as far as 25 years ago ! Indeed, he assumed an hypothetic distribution of energies essentially given by The particular value δ = 0 obviously recovers BG statistics. Clayton showed that a small value of δ (δ ≃ 0.01) was enough to make the theory consistent with the experimental data that were available at that time. Quarati and co-workers remarked (preliminarily in 1996 [94], and in more refined calculations since then [95]) that, since the needed δ is very small, the ansatz distribution could as well be the power-law one which appears in the present formalism. By identifying the first corrections (to BG) of both distributions, they obtained Consequently, values of q quite close to unity are enough to fit the solar neutrino discrepancy. Once again, we verify the extreme efficiency that modifications of the statistics can have.

C. Peculiar velocities in Sc galaxies
From the data obtained by the Cosmic Background Explorer (COBE), it has been possible to infer the distribution of peculiar velocities of certain groups of spiral (Sc) galaxies (we recall that by peculiar velocity we mean the residual velocity after the global universe expansion velocity has been substracted). Bahcall and Oh [11]  In fact, all the attempts exhibit a long tail towards high velocities, whereas the experimental data show a pronounced cut-off at about 500 Km s −1 . It is relevant to mention that all the models that were used had several fitting parameters, and nevertheless could not get rid of the tail. A fitting was then advanced [96] using the present formalism with only two free parameters, one of them being q and the other one a characteristic velocity. The function that was used was the q-generalized Maxwell distribution, essentially corresponding to an ideal classical gas. The quality of the fitting is quite remarkable, far better than those corresponding to the already mentioned four attempts. Once again, one sees that modifications of the statistics can be sensibly more efficient than modifications of the model. A famous example along this line is provided by the completely different physics associated with a gas of free fermions or of free bosons, i.e., a Fermi-Dirac ideal gas or a Bose-Einstein ideal gas (same model but different statistics).

D. Nonlinear inverse bremsstrahlung absorption in low pressure argon plasma
Liu et al [12] provided in 1994 strong evidence of the existence of non-Maxwellian velocity distributions in a specific plasma experiment, where low pressure argon is exposed to pulsed discharges. During the afterglow, measurements of the inverse bremsstrahlung of intense microwaves is performed. The experimental setting is such that Coulombian collisions are dominant. The experimental data were fitted with the following flat-topped distribution: with m ≥ 2. Souza and myself [97] showed in 1997 that the same data can equally well be fitted with with q ≥ 1. Furthermore, if we expand both fitting functions in the neighborhood of the Gaussian case, we obtain that In both fittings, the exponents m and q depend on the microwave power. In order to discriminate between the two fitting functions, quite precise and systematic experiments would be needed, in particular exploring the actual dependence of the results on the power.

E. Cosmic background microwave radiation
The most accurate data concerning the cosmic microwave background radiation have been obtained with the FIRAS (Far-infrared absolute spectrophotometer) instrument in the COBE (Cosmic background explorer) satellite [98]. These data are known to follow, in the 2 − 20 cm −1 region, Planck's black-body law. In 1995, Sa Barreto, Loh and myself [99], as well as Plastino, Plastino and Vucetich [100] (and several others since then), analyzed within what precision one is allowed to assume q = 1. The result that has systematically come out from these analyses is |q − 1| < 10 −4 . If new observations were performed in the future which would be say 10 times more precise than the available ones [98], this bound would be attained. Consequently, we would know better within what degree of confidence extensive thermostatistics can be used for this cosmological problem. If q = 1 turns out to be clearly confirmed, it is not excluded that we would have to revise our notions about the structure of space-time at the appropriate scales (possibly, Planck's length). It might come out that the physics at that level are better described by finite-difference equations than by differential equations !

F. Electron-positron collisions
The electron-positron annihilation into a virtual photon and the subsequent creation of a quark-antiquark pair provides the cleanest environment for the hadroproduction. Each of the two initial partons begins a complex cascade related to the strong-coupling longdistance regime of Quantum Chromodynamics. A partially successful global description of the hadroproduction has been provided through a thermodynamical equilibrium approach, mainly that of Hagedorn in 1965 [101]. This theory provides the following prediction: where σ is the distribution of the transverse momenta p T , T 0 is a characteristic temperature which Hagedorn predicts to be independent from the electron-positron collision energy W in the mass center referential, and c is a constant. This theory fits the data quite well for small W , say W < 10 Gev, but exhibits a pronounced failure for W increasing up to say 160 Gev. Very recently, Bediaga, Curado and Miranda [102] have used, along Hagedorn's lines, the present generalized statistics. The results are indicated in Figs. 8 and 9. Remark that (i) q varies smoothly and monotonically with varying W (Hagedorn's theory is recovered in the W → 0 limit), and (ii) T 0 ≃ 0.11 Gev and practically independs from W as desirable from Hagedorn's arguments. These results can be considered as a strong evidence of the applicability of the nonextensive thermostatistics to specific anomalous systems.

G. Emulsion chamber observation of cosmic rays
Cosmic rays can be observed by using a variety of detectors (such as Pb detectors; see [103] and references therein). Typically, showers of (clustered or individual) elementary particles appear which start at the so called vertex. These vertex are localized at various depths. The distribution of their depths can be measured (see [104] and references therein for the measurements done at the Mount Pamir lead chambers) . This distribution was recently fitted by Wilk and Wlodarcsyc [105] with the q = 1.3 function which emerges within the present formalism.

H. Reassociation of heme-ligands in folded proteins
In the folded conformational state, proteins might exhibit fractal effects. One such case might be the time evolution of the re-association of molecules that have been taken away from their natural positions. For instance, if O 2 molecules are dissociated, through light flashes, from their natural F e positions in a heme protein and reach positions outside the heme pocket, they tend to start rebinding, and, for so doing, they might have to follow a fractal path, or be under the dynamical influence of fractal excitations (e.g., fractons). Anyhow, this re-association phenomenon has been lengthily studied by Frauenfelder et al [106]. If we define ξ ≡ N(t)/N(0) where N(t) is the number of molecules that have not yet re-associated at time t, the ξ(t) monotonically vanishes with t. The results obtained by photo-dissociating CO molecules from Sigma Type 2 sperm whale Myoglobin (Mb) dissolved in a glycerol-water solution are shown in Fig. 10. For times not too long, the experimental data have been fitted by Frauenfelder et al [106] with where t 0 and n smoothly depend on the temperature T . Bemski, Mendes and myself [107] have argued that, within the generalized formalism, the following equation naturally appears: Its solution is given by This expression recovers, for q = 1, the usual exponential relaxation, and reproduces the Frauenfelder form through the identifications 1/(q − 1) ≡ n and 1/[(q − 1)λ q ] ≡ t 0 . Besides reobtaining the Frauenfelder empiric law, the present scheme allows for a better approximation if a crossover is admitted. More precisely, the above differential equation can be generalized as follows: The general solution involves [107] a hypergeometric function. The fitting is shown in Figs. 10 and 11. A detailed model which would justify the above phenomenological differential equation would be welcome.

I. Diffusion of Hydra Vulgaris
Upadhyaya et al [108] are presently performing interesting experiments on Hydra Vulgaris (a cylindrical body column with inner and outer cells, respectively referred to as endodermal and ectodermal respectively) in physiological solution. The endodermal cells are more adhesive than the ectodermal ones. The authors have measured the velocity distribution P (|V y |) of the "vertical" component of the velocity during the diffusion of endodermal Hydra cells in an ectodermal aggregate. The results are presented in Fig. 12, where the velocity unit is 10 −6 m/hour and the probability is represented by the histogram of the number of counts. These results were fitted with with the values of (a, b, c) indicated in the figure. Through the identification we precisely have the law which emerges within the present formalism, namely with q = 1.53. The next desirable step of course is to formulate a specific model for Hydra which would lead to this law, but this remains to be done.

J. Citation of scientific papers
An interesting study was recently done by Redner [109], in which the statistics of citations of scientific papers is focused. He exhibited the number N(x) of papers which have been cited x times for two long series, namely one (6 716 198 citations of 783 339 papers) from the Institute of Scientific Information (ISI) and another one (351 872 citations of 24 296 papers) from the Physical Review D (PRD). As expected, in both examples, N(x) monotonically decreases with x. Redner fitted the (relatively) low-x data with a stretched exponential of the form with β = 0.44 and 0.39 for the series ISI and PRD respectively. Also, he remarked that the large-x data exhibit a power law, namely close to ∝ 1/x 3 . He argues that this different functional behavior for low and large values of x must reflect different phenomenologies in these two regimes. In contrast with this viewpoint, Albuquerque and myself [110] argue that this is not necessarily so since the data can be quite satisfactorily fitted with a single function, namely with q = 1.53 and 1.64 for the series ISI and PRD respectively: see Fig.13 The satisfactory quality of the fittings is, after all, not so surprising, since we have mentioned earlier in this paper the connection [84] of this formalism with the Zipf law.

K. Electroencephalographic signals of epilepsy
It is since long known that the analysis of signals can be done within formalisms which use entropic forms. One such application has been recently done on EEG records of epilleptic humans and turtles [111]. The simultaneous use of wavelet-based multiresolution analysis including the nonextensive entropy S q leads to signals whose interpretation can be clinically neat and pharmacologically convenient. The authors of this novel processing suggest perspectives for building up automatic detection devices.

L. Cognitive psychology
The development of artificial neural networks and their connections with statistical mechanics (e.g., the Hopfield model for associative memory) makes quite natural the approach of cognitive problems with the present nonextensive formalism. Within this philosophy, we performed [112] an experiment of learning/memorization (of 5 × 5 and 7 × 7 square matrix having circles and crosses randomly distributed once for ever) with students of the University level; 150 students were interviewed, the first 30 in order to optimize the experimental protocole, and the other 120 to make the measurements of the time-evolution of the total amounts of errors when the original matrix was successively shown and hidden. The average results were then fitted with those obtained, for the same task, with a learning machine [113] having a perceptron architecture and an internal dynamics based on the Langevin equation [44] generalized by Stariolo to arbitrary q. The (average) learning time of the machine turned out to monotonically increase with q, exhibiting a practically divergent derivative at q = 1. The best human-machine fit occurred for q slightly above unity. More experiments and comparisons along these lines would be very welcome. Indeed, they would help better understanding some cognitive phenomena, on one hand, and could generate efficient machines for specific tasks, on the other.

M. Turbulence and time evolution of financial data
In 1996 Ghashghaie et al [114] compared financial data with those obtained from turbulent behavior and showed very similar behaviors when appropriate scalings are used. Ramos et al [115] have recently shown that all these data can be satisfactorily fitted with the functional form which emerges from the present formalism. Olsen and Associates data containing bid-ask quotes for US dollar-German mark exchange rates (1,472,241 records) are presented in Fig. 14 (probability density P ∆t (∆π) of price changes; ∆π = π(t) − π(t + ∆t) with ∆t = 640s, 5120s, 40960s, 163840s from top to bottom in the figure). The turbulent flow data [116] are presented in Fig. 15 (probability density P ∆r (∆v) of velocity differences; ∆v = v(r) − v(r + ∆r) for spatial scale delays ∆r = 3.3η, 18.5η, 138η, 325η from top to bottom in the figure, where η is the Kolmogorov scale, i.e., the critical limit for occurence of viscous dissipation). All these data exhibit a slight left-right assymmetry, which has been taken into account by Ramos et al: they used the same q for both sides but different widths. Needless to say that specific models leading to these fitting functions are very welcome.

V. COMPUTATIONAL EVIDENCES AND CONNECTIONS
A. Thermalization of a hot gas penetrating in a cold gas In 1991, Waldeer and Urbassek [117] made, assuming d = 3 Newtonian mechanics, a computational simulation in which a certain amount of high energy particles penetrate into a cold gas and are thermalized through the interactions between molecules. The cold gas is initially put at BG thermal equilibrium at temperature T C . The high energy particles at time t = 0 are randomly distributed in energy at a quite high energy per particle. The interaction potential was assumed to be hard sphere at short distances and decreasing, at long distances, like r −α . They analyzed three typical situations, one with α → ∞, hence well above d (i.e., very short range interactions), the second one with α = 4 (i.e., short range interactions), and the last one with α = 8/3, which is below d (i.e., long range interactions). In their simulation, they follow the time evolution of the energy distribution of the hot particles. After a transient, this distribution evolves with a regular pattern. For α > d, this pattern basically is the BG distribution with a temperature T (t) which gradually approaches T C from above (with lim t→∞ T (t) = T C ), in other words, through curves which approximatively are straight lines in a log-linear plot. For α < d, this approximation occurs through curves which are close to straight lines... in a log-log plot ! (Notice that the curvature in log-log plots tends to be upwards for α < d, whereas it is downwards for α > d; see Figs. 1, 2 and 3 of [117]). This power-law behavior is typical of q > 1. This pecualiarity was invoked by Koponen [9] in 1997 as a justification for using the present generalized formalism to discuss electron-phonon relaxation in ion-bombarded solids if the interactions are long-ranged. A study like that of Waldeer and Urbassek [117] which would systematically address the details of that thermalization by gradually varying α across d is missing and would certainly be very welcome.

B. Long-range classical Hamiltonian systems: Static properties
Let us focus here on what we refer to as weak violation of BG statistics. We use this expression to distinguish it from what we call strong violation of BG statistics. Both of them lead to nonextensive quantities, but, whereas the strong violation concerns q = 1, the weak one concerns q = 1 calculations. To make all this explicit we shall here focus on classical systems, i.e., all observables are assumed to commute. Let us consider the following paradigmatic Hamiltonian: where m is a microscopic mass, {p i , r i } are the d-dimensional linear momenta and positions associated with N particles, and r ij ≡ r j − r i . A typical situation is that of a finite confined system but, if some care is taken, the system could as well be thought of as having periodic boundary conditions. To be specific, let us assume where, in order to avoid any singularity at the origin (for any dimension d not exceedingly high), we have assumed, for the repulsive term, the Lennard-Jones exponent 12. What we desire to focus on in the present discussion is possible singularities associated with infinite distances, i.e., the effects of long-range (attractive) interactions. The case (α, d) = (6, 3) precisely recovers the standard Lennard-Jones fluid; the case (α, d) = (1, 3) is asymptotically equivalent to Newtonian gravitation; the case (α, d) = (d − 2, d) is asymptotically equivalent to d-dimensional gravitation (i.e., the one associated with the solutions of the d-dimensional Poisson equation); the case (α, d) = (3, 3) basically reproduces the distance dependance of permanent dipole-dipole interaction. The range of the (attractive) interaction increases when α decreases; α → 12 corresponds to very short-ranged interactions, whereas α = 0 corresponds to the situation of the Mean Field Approximation, where every particle (attractively) interacts with every other with the same strength, in all occasions. A typical quantity to be calculated within BG statistics is the following one (basically related to the T = 0 internal energy per particle): where the distances r have been expressed in units of a characteristic length of the problem. We immediately verify that this integral converges if α > d, and diverges if 0 ≤ α ≤ d.
Consequently, thermodynamic calculations in the 0 ≤ α ≤ d case have to be done with some care, and not blindly following the standard rules associated with BG statistics (i.e., q = 1). It is in this sense that we use the expression "weak" violation of BG statistics. The care to which we refer is the fact that we have to strictly consider the finite size of the physical system. Consistently, a relevant quantity that emerges naturally is We can check that, in the N → ∞ limit, we have As it will become transparent later on, what these regimes imply is that the system is extensive for α/d > 1 (hence standard thermodynamics apply), whereas it is nonextensive for 0 ≤ α/d ≤ 1, and special scalings become necessary [118,119] in order to have both a mathematically well posed problem, and a physical unfolding (or qualification) of the nonextensive region. The α/d > 1 regime has since long been analyzed [120], and it is well known that extensivity (or, stability, as also referred to) is lost for α/d ≤ 1. However, to the best of our knowledge, the scalings associated with the quantity N * , as well as its numerically efficient collapsing properties, were introduced for the first time by Jund et al [119] in 1995.
A quantity related to N * , namelyÑ turns out to be even more convenient. It is defined through [121] In the N → ∞ limit, we havẽ In the limit α/d → ∞,Ñ → 1; in the limit α/d → 1 + 0,Ñ ∼ 1/(α/d − 1); in the limit α/d → 0,Ñ ∼ N. Roughly speaking,Ñ characterizes the effective number of neighbors that can be associated with a given particle. This is the convenience to which we referred above.
We are ready now to present the kind of size-scalings we expect to be necessary for thermodynamically describing a generic classical Hamiltonian system with the type of interactions above mentioned. Let us focus on a simple fluid, and start with the standard case, i.e., α > d. Its Gibbs energy G (T, p, N) is given by where U, S, V, N, T and p respectively are the total internal energy, total entropy, total volume, total number of particles, temperature and pressure. In the N → ∞ limit, we obtain where the corresponding densitary variables have been introduced.
In contrast with the above, if we have 0 ≤ α ≤ d, the scalings are different, namely Consistently we have These equations recover the previous ones, i.e., those associated with the α > d, as a particular case. Indeed, for α > d,Ñ becomes a constant. So, we see that long range interactions have important thermodynamical consequences, namely (i) the energy quantities (G, U, which normally appear alone) that were extensive for α > d loose their extensivity; (ii) the non-energy quantities (S, V , which normally appear in canonical pairs with intensive quantities) that were extensive for α > d preserve their extensivity; (iii) the control parameters (T, p) that were intensive for α > d loose their intensivity. Consistently, to have mathematically well defined and physically useful equations of states and related quantities, everything must refer to finite quantities, hence, we must express all relations with the above rescaled variables. This does not imply that thermal equilibrium occurs through sharing equal values ofT ,p, etc. The zero-th principle of thermodynamics appears to hold in the usual way, even if we have long range interactions in the system. Although we have illustrated these features on a fluid, it is clear that the same considerations hold for all types of thermodynamical systems (magnets, dielectric substances, elastic solids, etc).
Two particular remarks must be made at this point: (i) When every element of the system equally interacts with each other (i.e., α = 0), N = N, and consequentlyT = T /N. In what concerns the thermostatistical approach of a system, this is equivalent to dividing the microscopic coupling constants by N, a familiar feature that is artificially imposed in all Mean Field calculations. We have used the word "artificial" because, whenever α ≤ d, the Hamiltonian which includes the microsopic interactions indeed is nonextensive! (and so is U). To divide the two-body coupling constants by N when α = 0 (or, byÑ, when α > 0) certainly is an artificial manner of forcing to be extensive a Hamiltonian which physically is not. This practice is traditionally frequent among magneticiens (who divide J by N), but certainly not among astronomers, who normally do not even think about what would be a very strange way of renormalizing the gravitational constant G ! (ii) If a singularity (for example, a critical phenomenon) occurs under particular physical conditions, it must generically occur at finite values ofT ,p, etc, and not at finite values of T, p, etc. Let us illustrate this for the simple case of a critical temperature:T c must be finite, hence, if α/d > 1, it must be T c (α, d) ∝ 1/(α/d − 1) This implies that T c must generically diverge for all classical systems at the extensive-nonextensive frontier. More precisely where A(d) is a system-dependent finite constant. In fact, let us anticipate that this precise behavior has been observed in the systems available in the literature [119] (d = 2 and d = 3 Lennard-Jones-like fluids, d = 1 and d = 2 Ising and Potts ferro-and antiferromagnets, etc), with no exception. To illustrate the connection between fluid models like the extended Lennard-Jones one above considered, and localized spin systems, let us briefly focus on the Ising ferromagnet. The simplest long-range N-spin cubic-lattice Hamiltonian of this kind is given by [122] where, for d = 1, r ij =1, 2, 3, ...; for d = 2, r ij = 1, 2 1/2 , 2, ...; for d = 3, r ij = 1, 2 1/2 , 2 1/3 , 2, ... ; and so on. Clearly, the limit α → ∞ recovers the first-neighbor d-dimensional spin 1/2 ferromagnet, whereas α = 0 corresponds to the Mean Field Approximation. For this model, kT c(α, d)/J diverges for 0 ≤ α/d < 1, decreases for α/d increasing above unity, and approaches the first-neighbor value (0 for d = 1, 2.269... for d = 2, etc) for α/d → ∞. Also, kT c /J ∼ A(d)/[(α/d) − 1] in the α/d → 1 + 0 limit. The introduction of T nicely enables the unfolding of the region where T c diverges. Indeed, kT c /J ≡ kT c /(JÑ) is finite in both extensive (α/d > 1) and nonextensive (0 ≤ α/d ≤ 1) regions, thus providing an enlightening unification. Finally, let us mention that, it seems that all equations of states (e.g., lim N →∞ M(T, N)/N can be, for all α/d ≤ 1, mapped into that associated with the Molecular Field Approximation. This simplifying feature appears to hold only for the static thermodynamic properties, and not for the dynamical ones, as will become clear later on (in Subsection V.G).
To close this subsection, let us emphasize that what we have been focusing on here is what we refer to as the weak violation of BG statistical mechanics (see Fig. 4). These are analytical or Monte Carlo q = 1 calculations (i.e., the energy distribution obeys the Boltzmann factor), but the variables must be scaled withÑ , which is not at all necessary for the standard, short-ranged interacting systems.

C. Long-range tight-binding systems
In the previous subsection we addressed classical systems. It is clear, however, that similar nonextensivity is expected to emerge in quantum systems if long-range interactions are present. One such Hamiltonian is the tight-binding-like which follows [123,124]: where c + i and c i are the creation-annihilation operators associated with electrons on site i, the {ǫ i } are the on-site energies, and V is the inter-site energy. The T = 0 electron diffusive properties corresponding to this Hamiltonian exhibit a variety of anomalies intimately related toÑ , as preliminary shown by Nazareno and Brito [123] and studied with more details in [124].

D. Granular systems
In 1995, Taguchi and Takayasu [8] simulated a vertically vibrated bed of powder with inelastic collisions and studied the distribution of horizontal velocities. In the lower layers (so called solid phase) they observed a standard Maxwellian (Gaussian) distribution of velocities. The situation was sensibly different in the upper layers (so called fluidized phase). Indeed, there the distribution was a Student's t-distribution, precisely the one appearing in Eq. (42) with q = 3, and assuming an energy proportional to the (velocity) 2 (together with d = r = 2, hence a constant density of states). This anomaly must be related to the fractallike granular clusterization which occurs in real space [125] but a deep analysis would be welcome. Also would further simulations, for instance of the cooling type. Studies of such computational models, either externally forced or just left to their own isolated evolution, can provide important physical insights, especially if quantities like the energy distribution, the Lyapunov spectrum (or at least its maximum value) or possible multifractality are focused on.

E. d = 1 dissipative systems
One-dimensional maps constitute the simplest systems which might present chaos. Basically they consist of the following recurrent equation: where h(x; a) is a rather simple nonlinear function, and a is a control parameter. Typically, both x and a are real numbers (but higher-dimensional situations are of course possible, and frequently studied). The logistic map, for instance, exhibits this structure. Typically, for a < a c , the system exhibits simple orbits, the attractor being a cycle whose number of elements is finite. For a > a c , the system can exhibit attractors with an infinite number of elements. The value a c is the critical one, usually referred to as the chaos threshold; the associated attractor typically constitutes a multifractal characterized by the so called probabilities corresponding to the ratios of numbers of points belonging to each window) the t = 0 value of the BG entropy S 1 (0), which is going to be very close to zero (strictly zero in the (M, N) → (∞, ∞) limit). We then allow each of the N points to evolve according to the logistic map until an attractor is achieved. The entropy S 1 (t) will grow with t until arrival to a saturation value S 1 (∞) which depends on (M, N) (necessarily lim N →∞ S 1 (∞) < ln M). In the M → ∞ limit, the growth of S 1 (t) is in fact linear (see, for instance, [128], which enables the following characterization of the so called Kolmogorov-Sinai entropy: Quite generically, the Pesin inequality holds, which states (for one-dimensional nonlinear dynamical systems) that and K 1 = 0 otherwise. Since we are only presenting a sketch of the situation, let us from now on address those particular systems for which the equality holds [129]. For those it is and K 1 = 0 otherwise. This type of analysis is convenient either if we have simple orbits (i.e., strong insensitivity to the initial conditions, i.e., for λ 1 < 0) or if we have strong chaos (i.e., strong sensitivity to the initial conditions, i.e., λ 1 > 0). But this analysis is a very poor one at say the edge of chaos, where λ 1 = 0 and we have weak sensitivity to the initial conditions. It is to unfold this type of situation that S q becomes extremely useful. Let us show how. At the chaos threshold we have K 1 = λ 1 = 0. But if we follow the same procedure we just described for calculating K 1 , but using instead S q (t), an interesting phenomenon can be revealed, which we describe now. Let us first define the following generalized Kolmogorov-Sinai entropy: A value q * is generically expected to exist [130] such that, for q > q * (q < q * ), K q = 0 (K q diverges), and K q * is finite ! Furthermore, it can be argued [38] that the above Pesin equality can be generalized as follows: and K q = 0 otherwise. It is through this important type of (in)equality that the connection emerges between S q and the power-law time-dependence of the sensitivity to the initial conditions. The particular value q * above described is what was numerically calculated in [38,126,127], and satisfies the scaling (115) Some of the above statements can be trivially checked with the logistic map at its chaotic region (i.e., for λ 1 > 0). We know in that case that S 1 (t) ∝ t, hence (assuming the simple case of equiprobability) the total number of possibilities W (t) grows exponentially with t. For any q > 1, S q (t) is always bounded, then K q necessarily vanishes. For any q < 1, S q (t) grows like the 1/(1 − q) power of W (t), which in turn, as said before, grows exponentially with t, hence necessarily K q → ∞. We conclude that q * = 1 ! The same picture is expected to hold for weak chaos.

F. Self-organized criticality
In the previous example, fine tuning (e.g., a = a c ) is necessary to observe the anomalous (q = 1) behavior. Let us address dissipative systems with many degrees of freedom, very particularly those which do not need fine tuning. Would robust systems like those exhibiting self-organized criticality [131] (SOC) also present q = 1 behavior? The answer is yes, as it has been clearly exhibited in at least three computational systems, namely the Bak-Sneppen model for biological evolution, the Suzuki-Kaneko model for imitation games and the Bak-Tang-Wiesenfeld model for sandpiles [132]. In these systems, the Hamming distance plays the role played by ξ in the previous ones. Also, the relevance of the order of the t → ∞ and N → ∞ limits has been exhibited. Like in the conjectural Fig. 4, the q = 1 behavior is observed only in the lim t→∞ lim N →∞ order. On what model ingredients does q depend? Is a taxonomy in universality classes analogous to that of standard critical phenomena possible? Is a multifractal f (α L ) function hidden somewhere? Does a scaling law like that of Eq. (115) still hold? Again, additional computational effort is very welcome.

G. Long-range classical Hamiltonian systems: Dynamic properties
Let us finally focus on the "heart" of statistical mechanics, the dynamics of the systems on which Boltzmann himself was meditating, namely the Hamiltonian systems with many degrees of freedom. Although lots of interesting quantum nonextensive phenomena must exist, here we shall restrict ourselves to the classical ones. We expect them to be able to provide nonextensive anomalies in a kind of pure, or simpler manner. Since a classical canonical Hamiltonian must satisfy the Liouville theorem, the Lyapunov spectrum must be symmetric with respect to zero, the corresponding eigenvalues being necessarily coupled in pairs of positive and negative values with the same absolute value. Consequently, to study the sensitivity to the initial conditions it suffices to study the maximum Lyapunov exponent. If it is positive, the system will generically be strongly chaotic, and will therefore easily satisfy the ergodic/mixing hypothesis (equality of time and ensemble averages). If, however, the maximum Lyapunov exponent vanishes, the entire spectrum will necessarily vanish, hence the system will be, at most, weakly chaotic, and will therefore have difficulties in satisfying the ergodic/mixing hypothesis, at least at not extremely large times (reflecting the macroscopic size of the system). The d = 1 coupled planar rotators N-body model with a two-body coupling constant proportional to 1/r α (r ≡ distance between two given rotators) has been recently studied (for α = 0 in [133], and, for α ≥ 0 in [134]) in the microcanonical ensemble. It has been established that, above a critical (conveniently normalized) total energy, the maximum (conveniently normalized) Lyapunov exponent is, in the N → ∞ limit, positive (zero) for α > 1 (0 ≤ α ≤ 1). More precisely, this maximum Lyapunov exponent is proportional to 1/N κ where κ(α) appears to be a monotonic function which decreases from κ(0) to zero while α increases from 0 to 1, and remains zero for all α ≥ 1). It must be recalled that it is only for α > 1 that the standard BG prescriptions provide finite integrals in the relevant calculations. If we were to discuss the d-dimensional version of the same model, we would certainly have κ(α, d). It is certainly possible that it is κ(α/d), and it would be interesting to check such a hypothesis. If we were to consider not planar (like the XY ferromagnet) but rather three-dimensional (like the Heisenberg ferromagnet) rotators, would κ be insensitive to that, or it would depend on the specific model? What would happen if, instead of localized rotators, we were to consider a long-ranged Lenard-Jones-like fluid, or d-dimensional gravitation? All these questions are certainly interesting, and worthy of further computational efforts.
In a recent paper, an essentially α = 0 model was considered [135], and, under certain circumstances, a crossover was found between anomalous (at times smaller than τ (N)) and normal (at times larger than τ (N)) diffusion, with τ (N) ∝ N. What happens if α > 0?
What happens for other models? The behavior observed for the particular model that was studied is consistent with Fig. 4. But is it exactly that conjecture that is going on? Only the study of the energy distributions (of single particles or of relatively large subsystems of the N-body system) themselves can provide the answer. What about the distributions of velocities? Are they Maxwellian (i.e., Gaussian) for α/d > 1 and non-Maxwellian otherwise? Are they Levy's or Student's t-distributions for α/d ≤ 1? If so, what is the dependence of q(α, d) ? Maybe q(α/d) ? Are the associated fluctuations anomalously time-correlated? Can nonmarkovian processes be present when the system is nonextensive (i.e., when 0 ≤ α/d < 1)? Plenty of intriguing questions that, sooner or later, will have to be answered, mainly through computational work (at least the first approaches). Better sooner than later!

H. Optimization techniques; Simulated annealing
The so called Optimization problem consists basically in determining the global minimum (or minima, if degeneracy is present) of a given cost funtion E(x), where x is a discrete or continuous d-dimensional variable. This problem can become extremely complex depending on the cost function having a large number of local minima, and on the dimension d being high. For the ubiquitous cases (in physics, chemistry, neural networks, engineering, finances, etc) for which analytic discussion is not tractable, a variety of computational algorithms have been developed. A special place among these is occupied, because of its efficiency and paradigmatic value, by the Simulated Annealing (SA) introduced in 1983 by Kirkpatrick et al [136]. Its denomination comes from its total analogy with the well known annealing technique, frequently used in Metallurgy for making a molten metal to reach its crystalline state (global minimum of the relevant thermodynamic energy). In SA, one or more artificial temperatures are introduced and gradually cooled, acting as a source of stochasticity, extremely convenient for eventually detrapping from local minima. Near the end of the process, the system hopefully is in the attractive basin of one of the global minima, the temperature is practically zero, and the algorithm asymptotically becomes a steepest descent one. The challenge is to cool the temperature the quickest we can but still having the guarantee that no definitive trapping in any local minimum will occur. More precisely speaking, we search for the quickest annealing (i.e., in some sense approaching a quenching) which preserves the probability of ending in a global minimum being equal to one. SA strictly follows a BG scheme. Let us illustrate for continuous x. The system "tries" to visit, according to a visiting distribution assumed to be Gaussian in the neighborhood of its actual state. The jump is always accepted if it is "downhill", i.e., if the cost function decreases. If it is "uphill", the jump might be accepted with a probability given by the Boltzmann factor corresponding to that cost function. It has been shown that the probability of ending on the global minimum equals unity if T (t) decreases logarithmically with time t. This algorithm is sometimes referred to as Classical Simulated Annealing (CSA) or Boltzmann machine. We easily recognize that, if instead of decreasing, the temperature was maintained fixed, this procedure precisely would be the well known Metropolis et al one for simulating BG thermostatistical equilibrium.
This optimization machine has been generalized within the present statistics as follows [58]. The visiting ditribution is generalized to be a q V -Gaussian, and the acceptance Boltzmann factor is generalized to be a q A -generalized factor, where q V and q A respectively are the visiting and acceptance entropic index. The cooling schedule is generalized as follows: This is the Generalized Simulated Annealing. This machine is characterized by (q V , q A ). The choice (1, 1) corresponds to CSA, and the choice (2, 1) corresponds to the so called Fast Simulated Annealing (FSA). The CSA corresponds to a cooling given by The FSA corresponds to a faster cooling given by The limiting case q V = 3 corresponds to These particular cases illustrate the great computational advantage that can be obtained by speeding up the algorithm by conveniently choosing q A (see also [68]). In practice, a convenient choice for q A is slightly below 3. The choice of q V seems to be more modeldependent. Details can be seen in a by now vast literature [59][60][61][62][63][64][65][66], in which applications have been done and variations have been performed concerning a variety of classical and quantum physical problems, the Traveling Salesman Problem, and many others. The first application [67] in quantum chemistry concerned simple molecules of the series CH 3 − R and some others, including the H 2 0 3 one, by using the MOPAC program package. Nowadays, Straub (in Boston), Okamoto (in Okazaki), and Ellis-Mundim-Bisch (respectively in Chicago, Salvador and Rio de Janeiro) and co-workers are currently improving and applying these techniques to complex molecules such as polypeptides, in particular with the aim of studying the important, though hard, protein folding problem.

VI. FINAL REMARKS
Boltzmann-Gibbs statistical mechanics and standard thermodynamics do not seem to be universal. They have domains of applicability quite poorly known nowadays. The precise knowledge of the restrictions for their validity is conceptually and practically very important. A nonextensive generalization of these formalisms is now available [1,14]; see Table I. It has been developed to cover at least some of these difficulties. Several important types of systems have been focused on in the present paper, which should substantially clarify the situation. These efforts span a wide epistemological variety, which goes from clear-cut theories to phenomenological ones, to quite well or less well understood fittings and connections. As exhibited at length here, the areas on which this formalism has been satisfactorily applied includes physics, astronomy, chemistry, mathematics, biology, economics, linguistics, cognitive psychology, etc. However, in spite of all this sensible progress, some inter-related crucial points are still to be understood and established on a neat and transparent basis. These include (i) the zeroth principle of thermodynamics and its connections with the thermodynamic limit, properties which would in principle exhibit the mathematical connection of the weak violation of the BG statistics (i.e., introduction ofÑ within the q = 1 formalism) with the strong one (q = 1); (ii) the functional dependence of q on (α, d) for long-range interacting Hamiltonians (0 ≤ α/d ≤ 1) and its connection with anomalous diffusion; (iii) the physical interpretation of the {p i } distribution and of the escort one {P i }, as well as clear-cut prescriptions for using one or the other when fitting experimental data (in the meanwhile, it appears that {P i } is the one to be used for equilibrium distributions of kinetic, potential or total energies, whereas {p i } is the one to be used for real-space diffusive nonequlibrium phenomena, either Levy or Student's [137] distributions according to whether correlations between jumps are suspected to be absent or present); (iv) what are the generic physical conditions for using the microcanonical, canonical and grand-canonical ensembles, under what exact and fully specified conditions they are expected to be thermodynamically equivalent, and the possible relevance for the so called thermogravitational catastrophe; (v) the clear connection with microscopic dynamic properties such as (partial) lack of ergodicity and mixing, the generalization (to weak chaos) of the Pesin inequality, the complete domain of validity of scaling relations connecting q to multifractality, and the possible relevance for SOC, spin-glasses and similar phenomena; (vi) the clear physical connection with quantum groups and, in general, deformations of relevant Lie algebras, and through these, the possible relevance for quantum gravity and the deep (possibly discrete, multifractal-like) structure of space-time.
On speculative grounds, one might think of two conjectures, to be clarified (i.e., rigorously formulated), confirmed or refuted. The first of these conjectures can somehow (on intuitive grounds) be formulated as follows. Strongly ergodic/mixing phenomena are ubiquitous in Nature; essentially, they are driven by microscopic interactions which are short-ranged in space-time (short-range forces, short-range memory, nonfractal boundary conditions); their basic geometry tends to be continuous, Euclidean-like; their thermodynamics is extensive; their central laws (energy distribution at equilibrium, time-relaxation towards equilibrium) generically are exponentials; and their thermostatistical foundation is Boltzmann-Gibbs statistics (i.e., q = 1). But weakly ergodic/mixing phenomena also are ubiquitous in Nature (e.g., biological, socio-economical, human cognitive phenomena, etc); essentially, they are driven by microscopic interactions which are long-ranged in space time (long-range forces, nonmarkovian memory, fractal boundary conditions); their basic geometry tends to be discrete, multifractal-like; their thermodynamics is nonextensive; their central laws (energy distribution at equilibrium, time-relaxation towards equilibrium) generically are power-laws; and the thermostatistical foundation of (at least some of) them (hopefully !) is the q = 1 statistics. The allowance for nonextensivity, in general, and for a nonextensive entropy, in particular, appears to be the "price" to be payed in order that Boltzmann's "mechanical" (i.e., one system evolving along time) manner of thinking about macroscopic systems coincides, at the level of the concrete mathematical results to be compared with the experimental data, with Gibbs' "ensemble" (i.e., many systems at a fixed time) manner of thinking. This coincidence of results is, since one century, well known and understood for standard sytems. Our aim here is to extend it to a large variety of anomalous systems.
The second of these conjectures is, at the present moment, so hard to rationalize that I dare to mention it here only because, after having been exposed to so many mathematical and physical arguments (that have been included in the present review), the reader might accept to honor me with his (her) indulgence, and have a look at the following few, intuitive lines. I believe that a deep analogy (maybe a kind of isomorphism, through the use of mathematical structures like the co-homology groups) exists between crystallographic structures such as crystal -quasicrystal -fluid, and nonlinear dynamics such as integrable -(weak) chaotic -(strong) chaotic. In some sense, they appear as space and time versions of the same mathematical structures. The first case concerns crystals (i.e., d-dimensional Bravais lattices) and integrable dynamics (i.e., motion on simple orbits), and its essential invariance is the discrete translational one. The third case concerns strongly disordered systems like fluids (liquid, gases) and strongly chaotic dynamics, and its essential invariance is the continuous translational one. Finally, the second, and intermediate, case is by far the most subtle one (and probably this is why it is the one that humanity took the longest time to discover), and concerns quasicrystals (e.g., Penrose tilings, amorphous substances like glasses, spinglasses, and other structures known to have (multi)fractal scalings; probably most of the so called complex spatial phenomena belong to this group) and weakly ergodic dynamics (e.g., edge of chaos, strange attractors, self-organized criticality, probably most of the so called complex time phenomena); its essential invariance is the dilatation one. In the first case we have the (space or time) highest predictability, and statistical methods are out of place. In the third case we have the (space or time) lowest predictability, and statistical methods exhibit their full power. Finally, in the second case, we have an intermediate predictability, and the statistical methods have to be "intrinsically nonlinear" in some sense, in order to be applicable and useful. There will be no surprise for the reader if, at this point, I admit that I believe that the statistical mechanics to be associated with the third present case of course is the BG one, whereas it might be the q = 1 statistical mechanics the one to be associated with the present second case ! (see also [138], where some preliminary, but nevertheless concrete, calculations exhibit this kind of connections). Let us now remind that Wiles' 1995 celebrated proof [139] of Fermat's last theorem was deeply related to quasicrystals since it was based on the proof of the Taniyama-Shimura conjecture about modular elliptic curves and used certain Hecke algebras. Consequently, I hardly dare to explicitly state a simple and unavoidable corollary, namely that, if my present second conjecture turns out to be, in some nontrivial and precise sense, correct, then the q = 1 statistical mechanics must be related to Fermat's last theorem ! Through the complete analysis, in more detailed terms, of the various aspects tackled in the present review, we could learn a lot and, very especially, (precisely) when the celebrated Boltzmann factor is the correct theoretical description of natural systems at thermal equi-librium ! This famous and so useful factor would then become, not a "dogma", as referred to by Takens [140], but a theorem ! We acknowledge extremely fruitful discussions, along the years, with E.M.F. Curado, A. Plastino 42)) of the Boltzmann factor (recovered for q = 1) as function of the energy E at a given renormalized temperature T ′ , assuming a constant density of states. From top to bottom at low energies: q = 0, 1/4, 1/2, 2/3, 1, 3, ∞ (the vertical line at E/T ′ = 1 belongs to the limiting q = 0 distribution; the q → ∞ distribution collapses on the ordinate). All q > 1 curves have a (T ′ /E) q/(q−1) tail; all q < 1 curves have a cut-off at E/T ′ = 1/(1 − q).  Figure 4 -Central conjecture of the present work, assuming a Hamiltonian system which includes two-body (attractive) interactions which, at long distances, decay as r −α . The crossover at t = τ is expected to be slower than indicated in the figure (for space reasons). Figure 5 -The one-jump distributions p q (x) for typical values of q. The q → −∞ distribution is the uniform one in the interval [−1, 1]; q = 1 and q = 2 respectively correspond to Gaussian and Lorentzian distributions; the q → 3 is completely flat. For q < 1 there is a cut-off at |x|/σ = [(3 − q)/(1 − q)] 1/2 . Figure 6 -The q-dependence of the dimensionless diffusion coefficient ∆ q (width of the properly scaled distribution p q (x, N) in the limit N → ∞). In the limits q → 5/3 − 0 and q → 5/3 + 0 we respectively have ∆ q ∼ [4/9]/[(5/3) − q] and ∆ q ∼ [4/(9π 1/2 ]/[q − (5/3)]; also, lim q→3 ∆ q = 2/π 1/2 . Figure 7 -The index q to be associated with a person whose critical value corresponding to Eq. (69) is V c . People with q < 1 (q > 1) tend to avoid (seek) risks for that particular game. The case q = 1 corresponds to an ideally rational agent. Figure 8 -Distribution of the transverse momenta p T obtained in electron-positron frontal collisions of energy W varying from 14 to 161 Gev. The dotted line corresponds to q = 1 (i.e., a Hagedorn type of fitting as given by Eq. (87)) for all values of W . The solid lines correspond to q = 1 fittings. Figure 9 -The values of q and T 0 used in the fittings of Fig. 8. When W approaches zero, q approaches unity, i.e., Hagedorn's theory; T 0 is essentially insensitive to W , as physically desirable. Figure 10 -Time evolution of ξ ≡ N(t)/N(0) associated with MbCO in glycerol-water. Dots: experimental data. Dashed lines: fittings with Frauenfelder's empiric law (Eq. (88) or Eq. (90)). Solid lines: fittings with the solutions of Eq. (91) (see Fig. 11). Figure 11 -Temperature dependences of the parameters used to fit the experimental data of Fig. 10. Figure 12 -Distribution of the "vertical" velocities during diffusion of endodermal Hydra cells in an ectodermal aggregate. The abcissa units are 10 −6 m/hour. The fitting was obtained using q = 1.53 (see the text). Figure 13 -Distribution of ISI and PRD papers having received x citations. (a) and (b) exhibit the fittings in [109]; (c) and (d) exhibit our present fittings (see the text). Figure 14 -Distributions of price changes for US dollar-German mark exchange rates and fittings using assymetric q-distributions (see the text). Figure 15 -Distributions of velocity differences and fitting using assymetric qdistributions (see the text). TABLES TABLE I. Some useful formulae written, through Eqs. (9) and (10) (ln q Z q = ln qZq − βU q ) Internal energy U q = − ∂ ∂β ln q Z q Free energy F q = U q − T S q = − 1 β ln q Z q Anomalous diffusion probability distribution p q (x) =