Thermofractal Dynamics in the Multisource Thermal Model: A Unified Approach to Multiparticle Production

Introduction

Describing multiparticle production in high-energy proton-proton ($pp$), proton-nucleus ($pA$), and nucleus-nucleus ($AA$) collisions remains a central challenge in relativistic heavy-ion physics [1, 2, 3]. Over decades of experimental exploration at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), two distinct yet successful phenomenological paradigms have emerged to address different aspects of final-state hadron spectra. On one hand, the Multisource Thermal Models account for geometric, flow, and longitudinal expansion effects in relativistic heavy-ion and hadronic interactions. These models treat the emission volume as a spatio-temporal superposition of multiple thermal fireballs, participant partons, or expanding cylindrical sources [1, 4]. While visually intuitive and effective at capturing transverse expansion and longitudinal rapidity profiles, standard multisource formulations rely on tuning several independent source temperatures ($T_1,T_2,T_3,\ldots$) and relative normalisation weights to reproduce low-, intermediate-, and high-$p_T$ regions [1].

On the other hand, the Thermofractals and Tsallis Statistics, grounded in the non-extensive entropy proposed by Tsallis [5], the thermofractal framework [3, 6, 7] provides a thermodynamic mechanism for $q$-exponential energy distributions based on scale invariance. Systematic empirical analyses across RHIC and LHC energies [8, 9, 10, 11] have demonstrated that energy fluctuations in self-similar QCD phase space naturally generate $q$-exponential distributions in $p_T$ [3, 8]. Significantly, thermofractal theory predicts an entropic index $q = 8/7 \approx1.143$ derived from $SU(3)$ phase-space effective degrees of freedom [7]. The incorporation of $q$-statistics and thermofractal concepts in expanding geometry is grounded in fundamental non-equilibrium dynamics [3, 12, 13, 14, 15, 16], and studies on local temperature fluctuations demonstrated that non-extensive statistics emerge naturally from superstatistical thermal ensembles [12, 13, 14], operating in close connection with Beck-Cohen superstatistics [16]. Modern kinetic theory shows that non-extensive statistics also arise from non-local correlations in quantum field theory and non-linear transport dynamics governed by the Plastino-Plastino equation [17], which has been formally derived from the Boltzmann transport equation endowed with $q$-deformed algebraic correlations in self-similar thermofractal systems [19, 20, 21].

In this paper, we show that a thermofractal source governed by $q$ is mathematically equivalent to a continuous superposition of local Boltzmann sources with inverse-Gamma distributed temperatures, establishing a theoretical connection to the multisource thermal model.

Thermofractals and Multisource Models

Thermofractals and Tsallis Non-Extensive Statistics

Standard Boltzmann-Gibbs thermodynamics relies on additive entropy, $S=-k_B\sum_i p_i\ln p_i$, suitable for short-range interactions and uncorrelated phase-space volumes. In high-energy hadronic collisions, long-range color force correlations, scale invariance, and self-similar parton branching violate statistical independence [3, 5]. Tsallis non-extensive statistics generalizes the entropy through a non-extensive entropic index $q$:

which recovers standard Boltzmann-Gibbs entropy in the limit $q\to1$.

Thermofractals represent a physical realization of non-extensive thermodynamics within quantum chromodynamics [3, 6]. A thermofractal is defined as a thermodynamic system exhibiting self-similarity and scale invariance across its phase space [3]. When high energy is deposited into a collision volume, the system forms excited compound states (fireballs or partonic clusters) that recursively decay into smaller, self-similar sub-systems down through sequential generations [3, 7].

The total internal energy $U$ of a thermofractal system is partitioned between the kinetic energy $K$ of $N_{\mathrm{sub}}$ constituent sub-systems and internal energy $\varepsilon_i$ stored within internal compound states [3]:

where $K=\sum_{i=1}^{N_{\mathrm{sub}}}K_i$. The kinetic state density is given by $\rho(K)=A_0K^{\alpha-1}$ with $\alpha=N_{\mathrm{sub}}\cdot d$, where $d$ denotes the spatial dimensionality of the kinetic phase space [3]. The scale-invariant phase-space density $\Omega(U)$ obeys the recursive convolution integral [3, 6]:

The internal density of states scales as a power law governed by an effective phase-space parameter $D$: $\Omega(\varepsilon_i)\propto\varepsilon_i^{D-1}$ [3, 7]. In the thermodynamic limit, this power-law phase-space density yields the non-extensive energy distribution $P(E)$ [3]:

where $C_q$ is the global normalisation factor, $T_0$ is the scale temperature of the primary system, and $q$ is linked to $D$ by [3, 7]:

Non-extensive statistical mechanics in expanding hadronic systems can also be formulated via non-linear transport dynamics [17]. The distribution function $P(E,t)$ evolves according to the generalized Plastino–Plastino Fokker–Planck equation [17]:

In momentum space, the distribution $P(p,t)$ evolves according to the generalised Plastino–Plastino Fokker–Planck equation, with linear drift $A(p)=-\gamma p$ and constant diffusion coefficient $B(p)=B$. The stationary solution depends on the drift and diffusion coefficients and on the boundary conditions. Its relationship to the thermofractal energy distribution in Eq. (4) must therefore be established separately. The stationary solution of Eq. (6) yields $q$-exponential energy distributions [12, 17], and has the form given in Eq. (4).¹ Nonextensive dynamics has also been used to study rapidity distributions in high-energy collisions [18].

The Multi-Source Thermal Model: Physical and Geometric Foundations

The multi-source thermal model provides a macroscopic hydrodynamic and geometric framework to describe multiparticle momentum and rapidity distributions in high-energy hadronic ($pp$) and relativistic heavy-ion ($AA$) collisions [1, 2]. Rather than assuming that the entire interaction volume forms a single static fireball in global equilibrium, the multi-source model treats the emitting system as an inhomogeneous, spatio-temporal superposition of multiple localised emission sub-sources [1, 2].

In the multi-source model, the total single-particle differential yield across transverse momentum $p_T$ and rapidity $y$ is given by a weighted linear combination of $N_{\mathrm{src}}$ discrete thermal sources:

where $k_i$ is the relative yield contribution of the $i$-th source ($\sum_{i=1}^{N_{\mathrm{src}}} k_i=1$) [1].

For a single thermal source at local temperature $T_i$ and chemical potential $\mu$, particle emission follows relativistic phase-space kinetics. The local invariant thermal spectrum $S_i$ takes the form[2]:

where $C(T_i)$ is the temperature-dependent normalization factor. When superimposing thermal sources, distinction must be made between combining unnormalized Boltzmann exponential factors—which directly maps onto the Euler integral representation—and combining fully normalized thermal yields $S_i$, where $C(T_i)$ introduces additional temperature weighting into the superposition kernel.

When collective transverse radial expansion with velocity $\beta_T$ ($\gamma_T=1/\sqrt{1-\beta_T^2}$) is present, the effective single-particle energy in the rest frame of the expanding fluid cell is Lorentz-boosted:

The longitudinal extension of source rapidities $y_b$ is parameterized by a Gaussian profile centered at mid-rapidity $y_0$ with width $\sigma_y$:

Although multi-temperature parameterizations ($N_{\mathrm{src}}=2$ or $3$, with $T_1\approx100$–$120$ MeV, $T_2\approx180$–$220$ MeV, and $T_3\ge300$ MeV) fit transverse momentum spectra across soft and hard regimes, the conventional multi-source model treats $T_i$ and $k_i$ as independent free parameters, lacking a microscopic thermodynamic principle that governs their discrete values [1].

Decomposition of Thermofractals into Thermal Superpositions

Using Euler’s integral identity for the Gamma function,

we set $x=1+(q-1)\frac{E}{T_0}$ and $a=\frac{q}{q-1}$, where $E=m_T\cosh y-\mu$. Substituting these definitions into the equation above yields:

Defining the fluctuating local temperature variable $T_s=\frac{T_0}{(q-1)t}$, we transform the integration variable from $t$ to $T_s$:

where $f(T_s;T_0,q)$ is the normalized inverse-Gamma probability density function [12, 16]:

Eq. (13) demonstrates that a $q$-exponential thermofractal factor maps into a continuous superposition of local Boltzmann exponentials weighted by an inverse-Gamma distribution $f(T_s;T_0,q)$ [12, 14, 16]. Note that if local Boltzmann sources are normalised individually prior to superposition, the temperature-dependent normalisation factor $C(T_s)$ is absorbed into an effective weight function $\tilde{f}(T_s) \propto C(T_s)f(T_s;T_0,q)$, shifting the effective quadrature nodes accordingly.

Defining the inverse temperature variable $\beta\equiv1/T_s$, the distribution $f(T_s;T_0,q)dT_s$ converts into a Gamma distribution $g(\beta)d\beta$:

The expectation value of $\beta$ is $\langle\beta\rangle=\beta_0=q/T_0$. Within Beck-Cohen Superstatistics [12, 16], the relative variance of inverse-temperature fluctuations across local sources is given by:

For $q=8/7$, the relative fluctuation variance is fixed to $(q-1)/q=1/8=12.5\%$.

To evaluate or discretise the exact continuous temperature field given by Eq. (14) for practical applications in multisource model fits, one can apply a discrete quadrature approximation. This method provides an exact transition from the continuous thermofractal temperature field $f(T_s;T_0,q)$ to discrete $N_{\mathrm{src}}$-component multisource models.

Substituting $x=\frac{T_0}{(q-1)T_s}$ into Eq. (13) maps the thermal superposition integral directly onto the standard Gauss-Laguerre weight function:

where the shape exponent is uniquely fixed by the entropic index:

For $q=8/7$, the exponent evaluates to $\alpha=7$. The integral in Eq. (17) can be evaluated numerically with high precision using an $N_{\mathrm{src}}$-point Gauss-Laguerre quadrature:

where $x_i$ are the roots of the generalized Laguerre polynomial $L_{N_{\mathrm{src}}}^{(\alpha)}(x)$, and $w_i$ are the corresponding quadrature weights. The discrete source temperatures $T_i$ extracted in phenomenological multisource fits thus directly correspond to the optimal quadrature nodes:

This mapping highlights the relationship between continuous thermofractal dynamics and discrete multisource models. While the representation of the $q$-exponential Tsallis factor as a continuous superposition of Boltzmann exponentials over a Gamma inverse-temperature distribution in Eq. (15) is exact, the discrete $N_{\mathrm{src}}$-component formulation obtained via Gauss-Laguerre quadrature in Eq. (19) is a finite-degree numerical approximation. For small numbers of sources ($N_{\mathrm{src}}=2$ or $3$), this quadrature provides an optimal discrete approximation to the continuous thermofractal field, rather than an exact algebraic identity.

Up to this point, the thermofractal superposition integral has been formulated in the rest frame of a single, static emission source. In relativistic hadronic and heavy-ion collisions, however, the continuous thermofractal temperature field $f(T_{\mathrm{s}};T_0,q)$ is embedded within a hydrodynamically expanding medium characterised by transverse radial flow and longitudinal expansion along the beam axis.

To account for collective radial expansion with transverse velocity $\beta_{\mathrm{T}}$ (Lorentz factor $\gamma_{\mathrm{T}} = 1/\sqrt{1-\beta_{\mathrm{T}}^2}$), assuming flow aligned with $p_{\mathrm{T}}$, the effective energy of an emitted particle in the rest frame of a local fluid cell is modified to $E_{\mathrm{s}} = \gamma_{\mathrm{T}}[m_{\mathrm{T}}\cosh(y-y_b)-p_{\mathrm{T}}\beta_{\mathrm{T}}]-\mu$, where $y_b$ denotes the longitudinal shift rapidity of the emitting cell. Convolving the single-cell emission yield over both the longitudinal source rapidity distribution $f_y(y_b)$ and the continuous thermofractal temperature distribution $f(T_{\mathrm{s}};T_0,q)$ yields the master differential spectrum:

For a longitudinally expanding system described by $N_{\mathrm{cyl}}$ cylindrical sources centered at rapidities $y_{b,j}$ with spatial widths $\sigma_{y,j}$ [1, 4]:

In the limit of vanishing radial flow ($\beta_{\mathrm{T}}\to0$, $\gamma_{\mathrm{T}}\to1$) and localized rest-frame rapidity emission $f_y(y_b)\to\delta(y_b)$, Eq. (21) reduces smoothly to the single-source Lorentz-invariant spectrum.

To assess the accuracy and range of validity of the discrete multisource representation, we compare, in Figure 1, the exact continuous thermofractal yield [Eq. (13)] with $N_{\mathrm{src}}$-component Gauss–Laguerre quadrature approximations [Eq. (19)] for $q=8/7$ ($\alpha=7$). We emphasise that the exact curve representing the nonextensive distribution provides a good description of the experimental data over 14 orders of magnitude of the measured cross-section, covering the $p_{\mathrm{T}}$ range $0<p_{\mathrm{T}}<1000$ GeV [26]. The multi-source, as well, model has been tested in several works (see [27, 28] and references therein).

For $N_{\mathrm{src}}=2$, the roots of $L_2^{(7)}(x)$ yield discrete temperatures $T_1=1.1667T_0$ and $T_2=0.5833T_0$ with relative normalization weights $k_1=2/3$ and $k_2=1/3$. In Figure 1, we show the invariant differential yield $\mathrm{d}^2N/(\mathrm{d}p_{\mathrm{T}}\mathrm{d}y)$ as a function of $p_{\mathrm{T}}$ for $N_{\mathrm{src}}=2$, 3, and 4 compared against the exact continuous thermofractal result, alongside the corresponding relative deviations $(Y_{\mathrm{quad}}-Y_{\mathrm{exact}})/Y_{\mathrm{exact}}$.

The comparison demonstrates that while a two-source formulation captures the low- and intermediate-$p_{\mathrm{T}}$ spectrum within 5% accuracy up to $p_{\mathrm{T}}/T_0\approx10$, increasing to $N_{\mathrm{src}}=3$ or 4 extends the agreement to higher $p_{\mathrm{T}}$ scales. At ultra-high $p_{\mathrm{T}}$, any finite quadrature eventually underestimates the exact power-law tail, reflecting the physical necessity of the full continuous temperature field $f(T_{\mathrm{s}};T_0,q)$ for high-$p_{\mathrm{T}}$ hard scattering.

Discussion and Conclusions

We have established a formal equivalence connecting thermofractal dynamics, Superstatistics, and the multisource thermal model. By demonstrating that $q$-exponential distributions map into a continuous superposition of Boltzmann sources through an inverse-Gamma temperature density $f(T_s;T_0,q)$, we link microscopic phase-space self-similarity [3, 7] with macroscopic multi-source expansion geometries [1, 4].

Key insights obtained from this synthesis include:

1. Parameter Constraint: Fixing $q = 8/7 \approx1.143$ from $SU(3)$ thermofractal theory constrains the temperature distribution shape $f(T_s)$, replacing unconstrained multi-temperature parameters $(T_1,T_2,T_3,\ldots)$ with a single primary scale parameter $T_0$.

2. Fluctuation Variance: The relative variance of inverse-temperature fluctuations across local sources is fixed to $(q-1)/q = 1/8 = 12.5\%$.

3. High-$p_T$ Asymptotics and pQCD Connection: In the ultrarelativistic regime $(p_T \gg m_0)$, the single-particle differential yield scales asymptotically as:

For $q = 8/7$, this power-law falloff $(p_T^{-6})$ is in direct agreement with the dimensional counting rules and perturbative QCD (pQCD) parton $2 \to2$ hard-scattering expectations observed across RHIC and LHC $pp$ spectra. This demonstrates that thermofractal dynamics smoothly bridges soft thermal particle production at low $p_T$ and hard partonic scattering tails at high $p_T$ without requiring ad hoc model switches. Furthermore, while finite Gauss-Laguerre quadratures $(N_{\mathrm{src}} = 2, 3)$ serve as effective discrete approximations capturing the low- and intermediate-$p_T$ regions, they truncate the spectrum at high momentum. Resolving the exact, continuous power-law asymptotics $(p_T \to\infty)$ highlights the physical necessity of the full, continuous inverse-Gamma distribution $f(T_s;T_0,q)$.

4. Medium Effects and Nuclear Modification Factor ($R_{AA}$): In heavy-ion collisions ($AA$), a dense partonic medium modifies particle production relative to vacuum $pp$ interactions. Within the thermofractal paradigm, medium effects do not simply wash out non-extensive statistics towards BG equilibrium; instead, partonic energy loss, medium-induced jet quenching, and non-local correlations alter the energy structure across subsystem generations. Recent transport analyses based on the Plastino–Plastino Fokker–Planck equation have demonstrated that the nuclear modification factor $R_{AA}(p_T)$ and energy-loss transport parameters (such as the jet quenching coefficient) can be rigorously calculated by modelling the anomalous diffusion of partons traversing the non-additive QGP [24, 25]. The continuous inverse-Gamma distribution $f(T_s;T_0,q)$ established in this work serves as the exact baseline vacuum field, upon which medium-drag interaction channels modulate the effective energy scale $T_0 \to T_0^{\mathrm{med}}(p_T)$, providing a unified physical connection between $pp$ vacuum spectra, transport dynamics in expanding fireballs, and $R_{AA}$ suppression patterns observed across RHIC and LHC energies.

In summary, the unified thermofractal-multisource framework provides a consistent mathematical bridge connecting scale-invariant QCD phase-space dynamics, temperature-fluctuating Superstatistics, and multi-source hadronic expansion models.

Acknowledgments

A. D. thanks the warm hospitality found at the Jan Kochanowski University. This work was partially supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and by the project RID/SP/00015/2024/01.