Fingerprints of nonequilibrium stationary distributions in dispersion relations

Distributions different from those predicted by equilibrium statistical mechanics are commonplace in a number of physical situations, such as plasmas and self-gravitating systems. The best strategy for probing these distributions and unavailing their origins consists in combining theoretical knowledge with experiments, involving both direct and indirect measurements, as those associated with dispersion relations. This paper addresses, in a quite general context, the signature of nonequilibrium distributions in dispersion relations. We consider the very general scenario of distributions corresponding to a superposition of equilibrium distributions, that are well-suited for systems exhibiting only local equilibrium, and discuss the general context of systems obeying the combination of the Schrödinger and Poisson equations, while allowing the Planck’s constant to smoothly go to zero, yielding the classical kinetic regime. Examples of media where this approach is applicable are plasmas, gravitational systems, and optical molasses. We analyse in more depth the case of classical dispersion relations for a pair plasma. We also discuss a possible experimental setup, based on spectroscopic methods, to directly observe these classes of distributions.


Nonequilibrium stationary distributions
To set the scene, let us first state precisely the statistical conditions we are considering here, i.e., the circumstances under which Eq. (1) holds. We are dealing, quite generally, with nonequilibrium systems in a stationary state that exhibit only local equilibrium. In this situation, our nonequilibrium system may be virtually divided up into small cells or small regions that remain infinitely close to equilibrium. In each cell, one has local thermal equilibrium and the (local) statistics f 0 (v|β) correspond to those of equilibrium statistical mechanics (in fact, as each cell exhibits local equilibrium, the whole machinery of equilibrium statistical mechanics holds locally). At larger scales however, one has to account for temperature inhomogeneities across different cells, that are encoded into f (β) . Provided that the temperature varies within a time-scale much larger than the local relaxation time (viz., the adiabatic Ansatz 41 ), the long-term distribution arises as a superposition of local equilibrium statistics, averaged over the distribution of the inverse temperature β ≡ 1/k B T (hereafter, k B = 1 ) across the different cells, i.e., Eq. (1). In the statistical mechanics literature, such an approach goes by the name of superstatistics 22 , as it merely consists of a superposition of different statistics. As we will discuss later on, the conditions outlined above are in fact very common and this methodological attitude allows describing a wide range of physical systems under very typical conditions.
In principle, one may construct infinitely many distributions in the form of Eq. (1), by selecting appropriate distributions f (β) . It is known however that three main classes of f (β) emerge as a universal behavior for most experimentally relevant situations. Beside their experimental relevance, these universality classes have a transparent statistical origin that can be understood from the application of the central limit theorem 23,24 or from the maximum entropy principle 42 (see also Ref. 43 ). Assuming a (d-dimensional) MB distribution for f 0 (v|β) , we present in what follows these classes and the corresponding long-term velocity distributions.
(i) χ 2 superstatistics In this case, the inverse temperature β is distributed according to the χ 2 distribution (also known as Gamma distribution) of degree n: where β 0 ≡ �β� is the average of β . The corresponding long-term velocity distribution follows from Eq. (1) as which is equivalent to the Tsallis distribution (q-Gaussian), known in the paradigm of NSM. This can be made more transparent by adopting a slightly different parametrization; by defining an entropic index q := 1 + 2/(n + d) and an effective inverse temperature β := β 0 (n + d)/n ≡ 1/ T , Eq. (3) can be re-expressed in the more familiar form 21 as www.nature.com/scientificreports/ We note in passing that, in the statistics literature, distributions in the form of Eqs. (3)-(4) are known as Student's t-distributions 44 . Note that for large v, these distributions behave asymptotically as power-laws, which make them relevant for many different physical problems.
(ii) Inverseχ 2 superstatistics In this case, instead of β , it is the temperature ( β −1 ) itself that is χ 2 -distributed. Then, β follows an inverse-χ 2 distribution, The corresponding velocity distribution follows as where K α (x) is the modified Bessel function of the second kind. Asymptotically, these distributions exhibit exponential tails in the velocity 45 , which make them of less relevance in the usual scenario of space plasmas and gravitating system that are usually characterized by heavy-tailed distributions. Nonetheless, this type of exponential behavior has been observed for stationary distributions in the case of Vortex glasses and vortex liquids 46 , in fusion plasmas 27 , and in cancer disease-specific mortality probability distributions 28 .
(iii) Log-normal superstatistics In this case, β is distributed according to the log-normal distribution, with an average of β given by β 0 ≡ �β� = µe s 2 /2 . In this case, no closed-form expression for the corresponding velocity distribution B(v) is known to date. Therefore, we will be dealing with this last case numerically. These three classes of distributions cover the rich variety of distributions observed in stationary nonequilibrium systems and, in fact, have substantial empirical evidence; χ 2 superstatistics (or equivalently Tsallis statistics) have been observed in numerous systems, such as dusty plasmas 5 , cold atoms 6,7 , and spin glasses 8 . Experimental evidence for log-normal superstatistics has been reported in the context of Lagrangian and Eulerian turbulence 25,26 , gravitational systems 38 , space plasmas 47 , and air pollution statistics 33 , among other systems 48 , while systems possibly obeying inverse-χ 2 superstatistics have been presented in Refs. [27][28][29] .
So far, we have been considering the case of classical statistics. One may nonetheless follow similar lines of reasoning in the quantum context as well by identifying the local equilibrium distribution f 0 (v|β) to the Bose-Einstein (BE) or Fermi-Dirac (FD) distribution. The latter case is particularly interesting for our analysis, as we will focus our attention on plasma DRs. In this case, the local equilibrium FD velocity distribution reads as where ǫ = mv 2 /2 . At this stage, an important remark is in order. Note that, although phenomena of wave propagation, producing DRs, are intrinsically 1d processes, the distribution is fundamentally a 3d quantity. In the classical case of the MB distribution, this subtlety can be safely disregarded since the 1d MB distribution, obtained by integrating over the two directions perpendicular to the wave propagation, has the same form as the original 3d distribution. In the case of the FD distribution, the reduction of the 3d distribution to its projected 1d version has to be done more carefully, and reads as 49 where β F ≡ 1/T F , with T F being the Fermi temperature and v F the Fermi velocity. In this case, we are not in position to obtain closed-form expressions for the corresponding superstatistical distributions (1). The latter can however be easily computed numerically.
In Fig. 1, we show examples of 1d velocity distributions associated with the three universality classes of superstatistics, i.e., χ 2 [Eq. (2)], inverse χ 2 [Eq. (5)], and log-normal [Eq. (7)]. The top panel corresponds to a local equilibrium MB distribution while the bottom panel shows its generalization to the FD distribution (for T = 0.01T F and µ = 0 ). For a better comparison between the different universality classes, we parametrize the distributions by using a single parameter defined as q := �β 2 �/β 2 0 . The latter can be easily expressed in terms of the parameters of the corresponding distribution f (β) as 2β . www.nature.com/scientificreports/ for the three universality classes. This parametrization is convenient since the limit q = 1 corresponds to a vanishing variance of f (β) , which shrinks into a Dirac delta δ(β − β 0 ) , and equilibrium distributions are recovered in this case. Before closing this section, it may be appropriate to give these distributions a more empirical credit, by confronting them with direct observations, in a context relevant to our analysis. In Fig. 2a, we confront them with independent observations 40 of velocity distributions of a plasma under micro-gravity conditions, obtained through the PK-4 instrument on-board the International Space Station (ISS), that clearly exhibit a non-Boltzmann behavior. We show the best-fits obtained by a nonlinear regression method based on the Levenberg-Marquardt algorithm 50,51 (also known as the damped least-squares method), for the χ 2 and the log-normal universality classes, while we disregard the inverse-χ 2 class, which fails to describe the high-energy part of the observational data, since it exhibits an exponential decay. One may appreciate that both the χ 2 and the log-normal classes nicely fit the data. Note that an even better fit can be reached by considering (as usually done in reports) that the core population is described by a MB distribution, and using the superstatistical distributions to fit only the halo part, i.e., the high energy tails, as shown in Fig. 2b.

Generic dispersion relations: the Schrödinger-Poisson system
To put the discussion in a very general context, we consider here generic dispersion relations that apply to a wide class of systems, formally described by the combination of the Schrödinger equation and the Poisson equation. As a trivial example, one may think of the case of a quantum plasma, where the SP model reads as (10) q := �β 2 � χ 2 www.nature.com/scientificreports/ where C is the Coulombian potential and n 0 is the background ion density (with ions supposed at rest). Another typical example is that of a self-gravitating system, in which case one has where G is the gravitational potential and m is the mass of the particle. Such a formal analogy between plasmas and self-gravitataing systems is well-known and, indeed, is also manifest in the classical regime, i.e., in the Vlasov-Poisson kinetic approach. However, as recently observed by Mendonça 39 , this analogy goes even further and applies to other physical problems as well. This can be made more transparent by rewriting the SP model in the form of an integro-differential equation as where 0 is an arbitrary external potential. Upon choosing the functional form U r − r ′ = 1/ r − r ′ , one covers the case of (i) a quantum plasma ( g = e 2 /4πε 0 and � 0 = −e 2 /ε 0 n 0 ), (ii) a self-gravitating system ( g = −m 2 G and 0 = 0 ), and (iii) ultra-cold atoms in a magneto-optical trap ( g = Q/4π and 0 = −Qn eq , where Q is the effective atomic charge and n eq is the equilibrium density). The model also covers (iv) the case of a BEC, with or without dipolar interactions, for and g = 4π 2 am , where a is the scattering length, C dd is the dipolar interaction strength while θ and ϕ are orientation angles. By employing the quantum kinetic Wigner-Moyal approach 52,53 , one arrives at a generic form of the DR as follows where f 0 (v) is the unperturbed (1d projected) velocity distribution and U(k) is the Fourier transform of U r − r ′ (for plasmas, self-gravitating systems, and ultra-cold atoms, one has U(k) = 4π/k 2 ). For a system, initially in thermodynamic equilibrium, f 0 (v) is usually identified with the MB distribution, or the relevant quantum generalization thereof. For a system exhibiting only local equilibrium, as defined above, the unperturbed distribution  1), where the temperature distribution f (β) is given by one of the three universality classes. While computing DRs, one is most often interested in the case of excitations with large phase velocities, i.e., v ≪ ω/k . In this limit, Eq. (15) can be expanded as follows where ω 0 is the appropriate characteristic frequency, given by For a quantum plasma, it corresponds to the plasma frequency ω 0 ≡ ω p = e 2 n 0 /ε 0 m , while for a selfgravitating system it is given by 4πGmn 0 is the so-called Jeans frequency. Note that, in principle, DRs (16) fully characterize the velocity distribution through the complete (infinite) set of velocity moments v l . In practice however, one is interested in the experimentally relevant limit of excitations with large phase velocities, in which case Eq. (16) is truncated and contains information on the velocity distribution only through its first moments. That is, one cannot reconstruct the distribution function from the DR alone. This turns out to be rather advantageous if one is interested in the universal effect produced by nonequilibrium distributions in DRs. This is because, even in the absence of a closed-form expression for the distribution, as happens in the case of log-normal superstatistics, one can compute the corresponding DR as long as the moments of f (β) are known. In fact, one may observe that, for distributions in the form of Eq. (1), the velocity moments can be expressed as a combination of the moments of the local equilibrium distribution f 0 (v|β) with those of f (β) as follows where �·� f 0 stands for an average over the equilibrium distribution f 0 (v|β) and �·� f (β) is an average over f (β) . For the three universality classes, i.e., Eqs. (2), (5), and (7), one has that can be combined with those of the (1d) MB distribution (l even), to determine all superstatistical velocity moments in an exact fashion. In particular, the first two moments appearing in the DR (16) follow as where the auxiliary functions φ i (q) and γ i (q) read as and Above, i = 1, 2, 3 stand for χ 2 , inverse-χ 2 , and log-normal superstatistics respectively. We can convince ourselves that nonequilibrium distributions have indeed a universal effect on DRs by observing that,by construction, one has q := �β 2 �/β 2 0 ≥ 1 , which implies that φ i (q) ≥ 1 and γ i (q) ≥ 1 , where the equality holds at equilibrium (i.e., for q = 1 ). The effect of nonequilibrium distributions is that of increasing the velocity moments, that are (16) (1 < q < 3/2), www.nature.com/scientificreports/ larger than their equilibrium counterparts; a feature that may be attributed to the heavy tails that are typical of nonequilibrium stationary distributions.
Note that in the case of a local equilibrium distribution given by the FD distribution, the moments are difficult to obtain in an exact form. One may use however approximate expressions 54,55 that are in agreement with the exact forms that are known in the dilute (classical) and ultra-dense (completely degenerate) cases. For instance, for the kinetic energy (i.e., the second-order velocity moment), one usually employs the arithmetic sum of the MB thermal energy and the Fermi energy. Hence, as far as this approximation is valid, the above discussion is applicable to the FD case as well.
Note that the classical limit of the DR (15) is easily obtained by letting → 0 , and reduces to which coincides with the classical DR corresponding to the the Vlasov-Poisson kinetic treatment. In particular, in the case of a plasma, by considering the large phase velocity limit ( v ≪ ω/k ) and neglecting terms beyond the order k 2 , one obtains after making use of Eq. (21), the following DR which is a modified version of the Bohm-Gross DR that have been already discussed in 47 , where D ≡ ε 0 T 0 /n 0 e 2 is the Debye length defined at the mean temperature T 0 . In Ref. 47 , we have shown that DRs in the form of Eq. (25) fit nicely with the data of Van Hoven 18 that have been argued to be a manifestation of Tsallis statistics 19,20 . Here, we compare the DR (25) with the dispersion data of Derfler et al. 56 , obtained using a Langmuir probe. Best-fits obtained using the generic DR (25) are shown in Fig. 3, where we have normalized the frequency to ω 0 and the wavenumber to 1/ D , that is, � ≡ ω/ω 0 and K ≡ D k . One may appreciate that the universal effect induced by nonequilibrium distributions is in a good agreement with the data.
Prior to analyzing more general forms of DRs in the following section, it may be appropriate to comment briefly on the effect of nonequilibrium stationary distributions on the coupling parameter. In a plasma, the latter is defined as the ratio between the average electrostatic potential per particle U and the average kinetic energy K . Independently on the distribution, one can estimate �U� ≈ e 2 /4πε 0 r S where r s is the Wigner-Seitz radius defined as 4πr 3 S n 0 /3 = 1 . On the other hand, for a nonequilibrium distribution in the form of (1), one has �K� = m�v 2 �/2 = 3φ i (q)T 0 /m . The coupling parameter follows as That is, a plasma in a nonequilibrium stationary state tends to be less coupled than its equilibrium counterpart (at the same mean temperature); a feature that may be attributed to the haivy tails of the nonequilibrium distributions, i.e., to the high energy population.

Classical dispersion relations: the case of a pair plasma
In this section, we examine in further details more general forms of classical DR (24), and the corresponding Landau damping, in the case of a plasma. Instead of an electron-ion plasma, we consider the case of a pair plasma, composed of electrons and positrons (or equivalently, two ion species having the same mass 57 ). In fact, because of the mass symmetry and the charge anti-symmetry between the two species, electrons and positrons in general follow (apart from a small asymmetry) the same distribution with the same mean temperature. It is www.nature.com/scientificreports/ therefore meaningful to assign the same superstatistical distribution to both species. The two-species extension of the DR (24) reads as where the subscripts (−) and (+) stand for electrons and positrons, respectively, and ω 0 is the plasma frequency for a pair plasma, ω 0 = 2n 0 e 2 /ε 0 m. Let us first analyze, as previously, the high frequency limit ( v ≪ ω/k ) of the DR (27), i.e., Langmuir modes. In this case, by considering only the real part of the DR for the moment, we recover the modified Gross-Bohm relation (25), where the Debye length is modified accordingly for a pair plasma, i.e., D ≡ ε 0 T 0 /2n 0 e 2 . Note however that, because of to the singularity in the integrand of Eq. (27), the integral is not properly defined in the whole velocity domain. In fact, this singularity induces an imaginary part in the DR, which is responsible for the Landau damping. To study this process, we follow very standard lines (see for instance Ref. 58 ), by setting ω = ω r + iγ , and restricting ourselves to a small imaginary part, i.e., |γ | ≪ ω r . It is convenient to re-express the DR (27) as D(k, ω) = 0 and split it into a real and imaginary parts, i.e., D = D r + iD i , where D r is determined by the Cauchy principal value of the integral, while for D i , one retains the pole contribution to the integral. Then, the imaginary part of the frequency is given as 58 Following these lines, we have for the χ 2 and inverse-χ 2 classes: and whereas for the log-normal class, we are not in position to obtain a closed-form expression for γ . The latter can however be expressed in the following integral form which can be easily solved numerically. In addition to Langmuir modes, corresponding to the high frequency branch, one can also investigate the low frequency branch or ion-acoustic waves (IAWs). In this case, we employ the usual trick 59 that consists of assuming a difference in the (mean) temperatures of the two species, i.e., T + T − and considering a phase velocity such that v th,+ ω/k v th,− , where v th,± are the thermal speeds of the two species. In the end of the process, we let T + → T − = T 0 . In this case, from the real part of the DR (27), we obtain where we use the subscript r to indicate that it corresponds to the real part of the DR. From the imaginary part, we obtain for χ 2 and inverse-χ 2 , the following expressions and respectively, whereas for the class of log-normal superstatistics, we have the following integral form  www.nature.com/scientificreports/ At this stage, it is instructive to observe that the real part of the DR depends on the distribution only through the auxiliary function φ i (q) . This is no longer true for the imaginary part, i.e., for damped modes. In Fig. 4, we show the real part of the DRs (Brillouin diagrams) for the case of Langmuir modes [Eq. (25)] and for IAWs [Eq. (32)], with different values of the auxiliary function φ i (q) . One may see that, indeed, the effect of different nonequilibrium distributions on the DR is qualitatively the same. By inverting Eq. (22), one can deduce, for the different universality classes, the corresponding value of q := �β 2 �/β 2 0 , which measures the degree of temperature inhomogeneities.
The imaginary part of the DR is, on the contrary, more sensitive to the distribution. In Fig. 5, we show the Landau time, i.e., 1/|γ | , normalized to 1/ω r , for Langmuir modes [Eqs. (29)- (31)] and for IAWs [Eqs. (33)- (35)], for the three universality classes of superstatistics. Although one sees the same tendency for the three universality classes, it appears that one can, at least in principle, distinguish between them. In practice however, this may be a highly nontrivial matter, inasmuch as the three universality classes induce qualitatively the same effect.
Before closing this section, it should be noted that beside the two modes discussed above, which are longitudinal modes, one also has transverse electromagnetic waves or light waves. In this case, the classical kinetic treatment yields 59 where c is the speed of light in vacuum. In this case, bearing in mind that light waves are high frequency waves, we may Taylor expand the integrand in Eq. (36) for v ≪ ω/k and make use of Eq. (21), to obtain the following DR while there is no damping in this case. We will not discuss Eq. (37) any further since, because of the term k 2 c 2 , the corresponding DRs are notoriously weakly sensitive to the distribution function.

Thermal Doppler broadening
We discuss in this section a possible experimental setup, based on a spectroscopic method known as thermal Doppler broadening, that may serve to directly identify the three universality classes of superstatistical velocity distributions. The method allows constructing the velocity distribution of a gas given its absorption spectrum and is particularly useful in probing velocity distributions in astrophysical situations, such as interstellar gas clouds 60 , space and astrophysical plasmas 61 , solar flares 62 , etc. www.nature.com/scientificreports/ The mechanism behind this technique is quite simple. It is based on the broadening of the spectral lines caused by the Doppler effect on the velocity distribution of a gas of emitting particles. In fact, the velocity of each particle produces a Doppler shift of spectral wavelengths such that a characteristic wavelength 0 , for particles at rest in the observer's reference frame, is shifted to = 0 (1 + v/c) , where v is the velocity component along the line of sight in the reference frame of the observer. The velocity distribution f(v) of the emitting particles can be fully determined from (and has the same form as) the distribution f ( ) of the spectral shifts, since one has f (v)dv = f ( )d .
For a classical gas in equilibrium, the velocity distribution is a Maxwellian (Gaussian) distribution, and the distribution of the spectral shifts f ( ) has the same form. The same reasoning applies to the superstatistical velocity distributions (1). In this case, the corresponding distribution of the spectral shifts is given as where B(x) is given by Eq. (1), with f (β) corresponding to one of the three universality classes, i.e., Eq. (2), (5), or (7). It is standard practice to characterize the Doppler broadening by the full width at half maximum (FWHM), which in the Maxwellian case, reads as 63 It is a simple exercise to identify the signature of the superstatistical distributions on the FWHM. In the case of χ 2 superstatistics, the distribution of the spectral shifts is merely a power-law and the corresponding FWHM can be easily computed as In the limit q → 1 , the term in the parentheses in Eq. (40) reduces to ln 2 and Eq. (40) reduces to Eq. (39), expressed at the mean temperature T 0 . Note that, if one adopts the parametrization commonly used in Tallis  www.nature.com/scientificreports/ statistics, i.e., q ≡ 1 + (2q − 2)/(q + 1) and β ≡ (q + 1)β 0 /2 , Eq. (40) can be expressed in the language of NSM 64,65 as where is the q-logarithm and T ≡ 1/β . For the two other universality classes, namely the inverse-χ 2 and the log-normal superstatistics, we are not in position to obtain closed-form expressions for the FWHM since this requires the variables to be solved for in an essentially non-algebraic way. One may however estimate the corresponding FWHM numerically. For typical atomic masses and temperatures, one has and � / 0 is of the order of 10 −6 . In this case, for a Maxwellian distribution, one has using Eq. (39), � / 0 ≈ 2.3548 10 −6 . Table 1 shows the estimated values corresponding to the three universality classes of superstatistics for different values of q. One may see that the presence of temperature inhomogeneities tends to decrease the FWHM; a feature that can be attributed to the heavy tails characterizing nonequilibrium distributions. Note that, although the estimated values presented in Table 1 are meant as an illustration, deviations of this order can, in principle, be measured with the present experimental technology that enables low-uncertainty measurements of spectral lines. In fact, the present experimental technology of spectroscopy achieves signal-tonoise ratios that exceed 66,67 10 5 and a resolution of frequency as small as 68,69 10 −12 .
Finally, it should be stressed that, in addition to the thermal broadening discussed here, there exist other mechanisms, of different origins, responsible for the broadening of the spectral line width. One should mention the pressure broadening, due to collisions, and the quantum broadening, arising fundamentally from the uncertainty principle. These effects are known to produce a Lorentzian profile, and are totally independent on the thermal broadening 70 . This independence allows expressing the spectral line profile as a convolution of the two profiles. Using the superstatistical profile (38) for the thermal broadening, one may construct the corresponding Voigt functions, i.e., the joint profiles accounting for both thermal Doppler and Lorentzian broadenings (see for instance Ref. 61 ).

Conclusions
Given the growing empirical evidence in favor of distributions different from those predicted by equilibrium statistical mechanics, it becomes increasingly important to better understand their origins in physical problems. The best methodological attitude to obtain a complete and reliable picture of their origin consists in employing a synergistic approach, combining theoretical knowledge with experimental studies, involving direct and indirect measurements. This paper addresses an example of indirect measurements; those associated with dispersion relations (DRs). We have studied, in a quite general context, the signature of nonequilibrium distributions in DRs, by considering distributions in the form of a superposition of distributions, i.e., superstatistics. We focused our attention on the three universality classes of superstatistics, namely χ 2 , inverse-χ 2 , and log-normal universality classes, that have strong experimental evidence [23][24][25][26][27][29][30][31][32][33] . The extension of our analysis to other forms of distributions is nevertheless straightforward. We discussed the general context of systems obeying a combination of the Schrödinger and Poisson equations and studied more closely the case of classical DRs for a pair plasma. We have also presented a possible experimental setup to directly observe these distributions using spectroscopic methods.
Our analysis sheds light on the universal effect produced by nonequilibrium distribution on DRs. In fact, as DRs generally depend only on the first moments of the distribution, they do not contain the full information about it, and the signature of the latter is merely indicative of the heavy tails that are characteristic of nonequilibrium distributions.
This study may help probing the occurrence of nonequilibrium distributions in a variety of media exhibiting similar elementary excitations such as plasmons, hybrid-phonon modes, or Bogoliubov excitations in Bose-Einstein condensates. Possible new prospects for future research consist in going beyond the linear regime and  Table 1. Estimated values of � / 0 ( ×10 −6 ) corresponding to the three universality classes of supertatistics, for different values of q ≡ �β 2 �/β 2 0 .