Hyperstatistics: Hierarchical Averaging

• Hyperstatistics is a hierarchical framework that averages conditional statistics over fluctuating local conditions to yield effective q-exponential distributions.
• It employs a gamma distribution method to derive q-generalized Boltzmann factors, linking the framework with Tsallis nonadditive entropy.
• Applications include thermodynamic oscillators, turbulent flows, and extreme value analyses, providing a unified approach to complex systems.

Hyperstatistics denotes a class of higher-level statistical constructions in which observable distributions are obtained by averaging conditional statistics over fluctuating local conditions. In the 2026 statistical-mechanical usage, it is a framework for complex systems in which Boltzmann–Gibbs statistics breaks down in domains of the system; a gamma distribution of domain Boltzmann factors yields a closed-form $q$-generalized Boltzmann factor $B_q(E)=\exp_q(-\langle\beta\rangle E)$, and the construction is explicitly tied to Tsallis nonadditive entropy. In related literature, the same statistics-of-statistics logic appears as metastatistics of extreme values, where maxima are averaged over fluctuating block sizes and parent-distribution parameters, and as a universal law for relative quantities derived from multiplicative binning of $k/x$-type data. This suggests that “hyperstatistics” is best understood as a family of hierarchical averaging procedures rather than a single universally standardized formalism (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026, Boumali, 5 Jun 2026, Ignaccolo et al., 2012, Kossovsky, 2013).

1. Domain structure and departure from Boltzmann–Gibbs statistics

In the statistical-mechanical formulation, hyperstatistics starts from a heterogeneous system partitioned into sub-domains, domains, and a whole-system level. At the sub-domain level, intensive quantities such as inverse temperature $\beta_i$ still fluctuate. A domain is then characterized by an effective average inverse temperature $\langle\beta_i\rangle$, while the experimentally measured $\langle\beta\rangle$ is an average over these domain-level values. The framework is designed for systems with inherent non-Boltzmann-Gibbsian behaviour and is presented as a way to treat complex systems in which Boltzmann–Gibbs statistics does not hold inside each mesoscopic domain of the system (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026).

This construction is closely related to, but explicitly distinguished from, Beck–Cohen superstatistics. Superstatistics assumes local Boltzmann–Gibbs equilibrium in each domain and averages $e^{-\beta E}$ over a probability density $f(\beta)$. Hyperstatistics instead assumes that the microscopic weights in each domain are already non-Boltzmannian, of Tsallis type, because a distribution of ordinary Boltzmann factors inside the domain has already been averaged. The word “hyper” is used in the literal sense of statistics of statistics: one first replaces ordinary Boltzmann factors by $q$-exponential weights at the domain level, and then considers a distribution of these generalized weights across domains (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026).

A technical nuance emphasized in the oscillator application is that the step

$$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$

is not treated as a general algebraic identity for nonlinear $B_q(E)=\exp_q(-\langle\beta\rangle E)$0. It is instead described as a domain-averaged effective replacement: after an internal transform produces $B_q(E)=\exp_q(-\langle\beta\rangle E)$1 in each domain, one replaces the domain mean by the measured system mean $B_q(E)=\exp_q(-\langle\beta\rangle E)$2. This distinguishes the formal derivation of the domain-level weight from the modeling step that produces the system-wide effective factor (Boumali, 5 Jun 2026).

2. Analytical construction of the $B_q(E)=\exp_q(-\langle\beta\rangle E)$3-generalized Boltzmann factor

The foundational calculation assumes a gamma distribution for the fluctuating inverse temperature within a domain,

$B_q(E)=\exp_q(-\langle\beta\rangle E)$4

with shape parameter $B_q(E)=\exp_q(-\langle\beta\rangle E)$5 and mean $B_q(E)=\exp_q(-\langle\beta\rangle E)$6. Its Laplace transform with kernel $B_q(E)=\exp_q(-\langle\beta\rangle E)$7 yields

$B_q(E)=\exp_q(-\langle\beta\rangle E)$8

With

$B_q(E)=\exp_q(-\langle\beta\rangle E)$9

this becomes the Tsallis $k/x$0-exponential,

$k/x$1

or, in the standard definition adopted by the papers,

$k/x$2

The resulting system-level generalized Boltzmann factor is written

$k/x$3

and, for the investigated Uniform, Gamma, Log-Normal, F, and $k/x$4-Gamma probability distribution functions, the functional form remains of $k/x$5-exponential type (Squillante et al., 23 Apr 2026, Boumali, 5 Jun 2026).

The same framework is formulated in entropy terms through Tsallis nonadditive entropy,

$k/x$6

with the Boltzmann–Gibbs limit recovered as $k/x$7. A central claim of the 2026 proposal is that hyperstatistics preserves the concavity of nonadditive $k/x$8-entropy. The relevant convergence analysis is expressed through a $k/x$9-generalized Gamma function,

$\beta_i$0

which converges only if

$\beta_i$1

This restriction simultaneously controls normalization, the existence of moments, and the admissible range of $\beta_i$2 in the hyperstatistical construction (Squillante et al., 23 Apr 2026).

A recurrent theme is structural independence from the detailed shape of the higher-level density $\beta_i$3. In the formulation used for relativistic oscillators, once domain-level Gamma fluctuations have produced the internal $\beta_i$4-exponential form, subsequent averaging over any normalizable $\beta_i$5 does not change that functional form except for rescaling by $\beta_i$6. The papers therefore present hyperstatistics as an analytically tractable route from fluctuating local Boltzmann factors to a closed $\beta_i$7-generalized Boltzmann factor (Boumali, 5 Jun 2026).

3. Thermodynamic realization in relativistic oscillators

A fully worked thermodynamic realization is given for the one-dimensional Klein–Gordon oscillator (KGO) and Dirac oscillator (DO). The KGO positive-energy spectrum is

$\beta_i$8

with degeneracy $\beta_i$9. For the 1D DO, the full spectrum is

$\langle\beta_i\rangle$0

and after reorganizing the positive-energy states one obtains the same energy levels as the KGO but with degeneracy

$\langle\beta_i\rangle$1

Thermodynamics is built from excitation energies

$\langle\beta_i\rangle$2

so that $\langle\beta_i\rangle$3. This removes the rest-energy shift and enforces third-law behaviour $\langle\beta_i\rangle$4 as $\langle\beta_i\rangle$5 (Boumali, 5 Jun 2026).

With

$\langle\beta_i\rangle$6

the partition function is

$\langle\beta_i\rangle$7

The canonical relations used are

$\langle\beta_i\rangle$8

and the derivatives of the $\langle\beta_i\rangle$9-exponential remain closed: $\langle\beta\rangle$0

$\langle\beta\rangle$1

The resulting formulas for $\langle\beta\rangle$2 and $\langle\beta\rangle$3 are explicit finite sums over $\langle\beta\rangle$4 and $\langle\beta\rangle$5, with KGO and DO differing only in the degeneracy structure (Boumali, 5 Jun 2026).

The reported thermodynamic behaviour is specific. Hyperstatistics successfully reproduces the high-temperature Boltzmann limit $\langle\beta\rangle$6, is structurally independent of $\langle\beta\rangle$7, avoids the unphysical negative regions of Beck’s polynomial bracket, and systematically distinguishes KGO from DO by capturing the enhanced entropy and sharper specific-heat structure caused by spin-induced degeneracy. The two frameworks agree quantitatively for $\langle\beta\rangle$8 and $\langle\beta\rangle$9, but diverge at high temperatures where Beck’s polynomial expansion loses validity and the exact hyperstatistical $e^{-\beta E}$0-exponential remains positive, monotonic, and analytic. For this oscillator problem, because $e^{-\beta E}$1, the sums converge for $e^{-\beta E}$2, and the numerical implementation in the paper uses $e^{-\beta E}$3 (Boumali, 5 Jun 2026).

4. Experimental and phenomenological applications

The 2026 papers present hyperstatistics as a unifying fit architecture across systems as different as RC discharge, cryogenic pumping, relativistic collision data, turbulence, Brownian motion, and dielectric response. In these applications the fluctuating quantity may be an inverse temperature, a relaxation time, an effective variance, or an analogous scale parameter, but the fitted observable is still written in $e^{-\beta E}$4-exponential form after hyperstatistical averaging (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026).

System	Hyperstatistical object	Reported parameterization or outcome
Capacitor discharge	$e^{-\beta E}$5	hyperstatistical fit gives $e^{-\beta E}$6 and $e^{-\beta E}$7
$e^{-\beta E}$8He cryostat pumping	$e^{-\beta E}$9	larger $f(\beta)$0
$f(\beta)$1-Pb midrapidity	$f(\beta)$2	$f(\beta)$3, $f(\beta)$4
Turbulent acceleration	$f(\beta)$5 or $f(\beta)$6-Gaussian form	long tails described without choosing between Gamma or Log-Normal
Brownian VACF	$f(\beta)$7	$f(\beta)$8, $f(\beta)$9, $q$0

For the RC circuit, the physical motivation is a distribution of relaxation times in a real dielectric, replacing the ideal single-$q$1 decay $q$2 by a $q$3-exponential decay. For pumping of $q$4He lines in a closed cycle cryostat, the same mechanism is used for broadly distributed effective pumping times. In midrapidity $q$5-Pb collisions at $q$6 TeV, hyperstatistics is used to fit transverse-momentum spectra with $q$7 and $q$8, consistent with the bound $q$9. For fully developed turbulence, the long-tailed acceleration distribution is fitted directly by a $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$0-exponential form, while for confined Brownian motion the normalized velocity autocorrelation function is fitted by component-dependent $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$1-exponentials before the regime dominated by long-time power-law tails (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026).

A distinct application concerns dielectric response. Using the $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$2-Gamma distribution for relaxation times, the long-$$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$3 tail behaves as

$$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$4

After integrating Debye’s expressions over this distribution, both $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$5 and $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$6 acquire low-frequency power-law behaviour,

$$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$7

which the paper identifies with the universal dielectric response (Squillante et al., 23 Apr 2026).

The Brownian-motion and neuroscience discussions extend the framework beyond static distributions. For the velocity autocorrelation function, hyperstatistics is presented as a compact description of anisotropic memory effects near a wall. For brain dynamics, the proposal is more programmatic: $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$8 and possibly $$\int_0^\infty f(\beta)\,\exp_q(-\beta E)\,d\beta=\exp_q(-\langle\beta\rangle E)$$9 are suggested as phenomenological indices of dynamical regime, with subcritical-like, supercritical-like, and altered states described in hyperstatistical terms. A plausible implication is that the framework is being positioned not only as a fitting formalism but also as a coarse-grained language for complex multilevel dynamics (Squillante et al., 17 Jun 2026).

5. Other hyperstatistical and hierarchical usages

A conceptually similar but mathematically distinct usage appears in the metastatistics of extreme values. There the object of interest is not a generalized Boltzmann factor but the distribution of block maxima conditional on local conditions. If $B_q(E)=\exp_q(-\langle\beta\rangle E)$00 is the number of wet days in year $B_q(E)=\exp_q(-\langle\beta\rangle E)$01 and $B_q(E)=\exp_q(-\langle\beta\rangle E)$02 are the parent-distribution parameters, the conditional cumulative distribution of the annual maximum is $B_q(E)=\exp_q(-\langle\beta\rangle E)$03. The metastatistical or hyperstatistical level is introduced through a factor $B_q(E)=\exp_q(-\langle\beta\rangle E)$04, producing

$B_q(E)=\exp_q(-\langle\beta\rangle E)$05

or, empirically,

$B_q(E)=\exp_q(-\langle\beta\rangle E)$06

The paper states explicitly that this is a textbook hyperstatistical mixture over fluctuating parameters $B_q(E)=\exp_q(-\langle\beta\rangle E)$07, and uses it to model annual rainfall maxima when the number of wet days and Weibull parameters vary from year to year (Ignaccolo et al., 2012).

The same paper contrasts this mixture-based treatment with classical generalized extreme value fitting. Its argument is that ordinary GEV methods fail when convergence is slow and the data are statistically inhomogeneous, whereas the metastatistical construction keeps the finite-$B_q(E)=\exp_q(-\langle\beta\rangle E)$08 structure and averages over inhomogeneous local conditions. Synthetic experiments and the Padova daily precipitation record are used to show that the metastatistics methodology produces correct predictions in cases where the traditional approach based on the generalized extreme value distribution does not (Ignaccolo et al., 2012).

Another higher-level usage concerns relative occurrences of quantities within a data set. Here the central object is the bin-proportion vector produced by multiplicative bin schemes. For an expanding scheme with $B_q(E)=\exp_q(-\langle\beta\rangle E)$09 bins per cycle and inflation factor $B_q(E)=\exp_q(-\langle\beta\rangle E)$10, the infinite-range $B_q(E)=\exp_q(-\langle\beta\rangle E)$11 model yields the General Law of Relative Quantities,

$B_q(E)=\exp_q(-\langle\beta\rangle E)$12

When $B_q(E)=\exp_q(-\langle\beta\rangle E)$13, this reduces to the Benford formula

$B_q(E)=\exp_q(-\langle\beta\rangle E)$14

The paper presents this as a higher-level descriptor of magnitudes shared by many physical and abstract data sets that are logarithmic in the Benford sense; from a hyperstatistics perspective, it is a meta-level law over summary statistics rather than a mixture over temperatures or relaxation times (Kossovsky, 2013).

These constructions share a common architecture: one first specifies a conditional or local statistical law, and then averages it over an additional layer of variability. What differs across the papers is the object being mixed—Boltzmann factors, extreme-value distributions, or bin-proportion structures. This suggests a broad family resemblance among the usages, even when the underlying mathematics is not the same (Ignaccolo et al., 2012, Kossovsky, 2013).

6. Interpretation, limitations, and unsettled terminology

Several interpretive cautions recur in the literature. First, hyperstatistics is not simply synonymous with superstatistics. The 2026 statistical-mechanical papers insist that superstatistics assumes Boltzmann–Gibbs statistics locally, whereas hyperstatistics assumes that a $B_q(E)=\exp_q(-\langle\beta\rangle E)$15-exponential structure has already emerged inside each domain from a distribution of Boltzmann factors. Second, the replacement of a domain average by $B_q(E)=\exp_q(-\langle\beta\rangle E)$16 is treated as an effective modeling step rather than a universal identity. Third, the parameters $B_q(E)=\exp_q(-\langle\beta\rangle E)$17 and $B_q(E)=\exp_q(-\langle\beta\rangle E)$18 may be linked by $B_q(E)=\exp_q(-\langle\beta\rangle E)$19 in the exact Gamma-transform construction, but in data analysis they are sometimes fitted independently and then checked against the convergence condition $B_q(E)=\exp_q(-\langle\beta\rangle E)$20 (Squillante et al., 23 Apr 2026, Squillante et al., 17 Jun 2026, Boumali, 5 Jun 2026).

The framework also has explicit mathematical and empirical limits. The oscillator study notes a convergence bound $B_q(E)=\exp_q(-\langle\beta\rangle E)$21 for sums over $B_q(E)=\exp_q(-\langle\beta\rangle E)$22. The general proposal identifies the choice of hyperdistribution $B_q(E)=\exp_q(-\langle\beta\rangle E)$23 or $B_q(E)=\exp_q(-\langle\beta\rangle E)$24 as a physical modeling question even though the final $B_q(E)=\exp_q(-\langle\beta\rangle E)$25-exponential form is claimed to be universal for the investigated families. The metastatistics of extremes assumes that local yearly parent distributions can be meaningfully parameterized, while the relative-quantities approach presupposes positive data spanning a sufficiently long range in magnitudes. In neuroscience, the application to brain dynamics is presented as potential rather than complete, and the papers explicitly call for systematic comparison across empirical EEG, MEG, and fMRI data and across different states of consciousness (Squillante et al., 17 Jun 2026, Ignaccolo et al., 2012, Kossovsky, 2013).

A common misconception is that positivity of the fitted curve alone guarantees physical admissibility. The comparison with Beck’s asymptotic superstatistics shows why the exact functional form matters: the truncated polynomial bracket can become negative for Gamma and Log-Normal cases or grow unphysically for the F-distribution, whereas the hyperstatistical $B_q(E)=\exp_q(-\langle\beta\rangle E)$26-exponential remains positive and monotonically decreasing for $B_q(E)=\exp_q(-\langle\beta\rangle E)$27 and $B_q(E)=\exp_q(-\langle\beta\rangle E)$28. The claimed advantage is therefore not only empirical fit quality but also analytic control over sign, monotonicity, and convergence (Boumali, 5 Jun 2026).

The terminology itself remains unsettled. One line of work uses “hyperstatistics” for a Tsallis-based two-level statistical mechanics of fluctuating Boltzmann factors; another treats metastatistics of extreme values as a hyperstatistical treatment of inhomogeneous maxima; a third uses a hyperstatistical perspective for universal laws of relative quantities. This suggests that the term currently functions less as a single standardized theory than as a label for hierarchical averaging procedures applied to heterogeneous systems (Squillante et al., 23 Apr 2026, Ignaccolo et al., 2012, Kossovsky, 2013).