Introduction

Experiments involving heavy-ion collisions at ultra-relativistic energies, conducted at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), aim at studying a new state of matter known as the quark–gluon plasma (QGP) [1-4]. Viscous hydrodynamics is remarkably successful at describing a wide range of observables and has become the “standard model” for the evolution of ultra-relativistic heavy-ion collisions [5-7]. Abundances of the final state hadrons contain important information about the dynamics of the QGP created in relativistic heavy-ion collisions. While the measured integrated yields are rather well understood within the picture of the thermal particle production, the shape of the transverse momentum distribution, in particular at low transverse momentum (p_T) below 500 MeV/c, is still poorly understood within the state-of-the-art hydrodynamic model calculations [6-9]. The particle production at low transverse momentum is associated with the long-distance scales, which are accessible in heavy-ion collisions and out of reach in hadronic interactions. Its enhancement could indicate transverse momentum and particle species yield redistribution of the thermally produced hadrons due to conventional phenomena [10-15], which is not yet implemented in the current state-of-the-art hydrodynamic models, or new physics phenomena [16-29].

An enhancement at low transverse momentum was first observed at the ISR in high multiplicity pp and α–α collisions when compared to minimum bias pp collisions [30] and in p–A collisions at Fermilab and CERN [31, 32], and later in A–A collisions at the AGS and CERN [21, 33-36]. A less than 10% enhancement at low-p_T was observed for midrapidity pions in both p–Pb and A–A collisions by the NA44 Collaboration [37]. The low-p_T enhancement in A–A collisions showed no target size dependence and was smaller for pions at midrapidity compared to target rapidity [21]. Intriguing explanations proposed at that time included exotic behavior in dense hadronic matter [17-19], the decay of quark matter droplets [38], collective effects [20], baryonic and mesonic resonance decays [21-23], and the possible formation of a transient state with partially restored chiral symmetry in the early stage of the heavy-ion collision [24-26].

The mechanism of low-p_T pion production remains an open question for experiments at the modern heavy-ion colliders. An excess is seen when comparing low-p_T pion yield measured by the ALICE [39-42] at the LHC and the PHENIX and STAR [43-46] at RHIC to the hydrodynamic model calculations [7-9, 47-53]. An indication of the excessive pion yield, though with large uncertainties, is also visible from the comparison of the thermal model fits to measured integrated yields of different particle species [12]. While in most experiments at RHIC and the LHC the pion p_T spectra are measured only above p_T = 0.1–0.2 GeV/c, the PHOBOS experiment at RHIC [54, 55] measured it down to the p_T=30–50 MeV/c. An extrapolation of the blast-wave model [55], fitted to PHOBOS experimental data in the intermediate p_T region, revealed no significant increase in kaon and proton production when compared to low-p_T data. However, the same extrapolation showed a possible enhancement in pion production at very low p_T.

The low-p_T pion excess may arise from physics mechanisms not accounted for in the current hydrodynamic model simulations, like Bose–Einstein condensation [27, 28], increased population of resonances [10], treatment of the finite width of ρ meson [11], or critical chiral fluctuations [29]. Quantification of the low-p_T pion yield excess is important for both, the improvement of fluid dynamic modeling and the search for new particle production mechanisms in heavy-ion collisions. On the experimental side, the proposed next-generation detector ALICE 3 at the LHC [56], which combines excellent particle identification capabilities, a unique pointing resolution, and large rapidity coverage, will allow measurements below p_T ~ 100 MeV/c.

To advance in understanding of the mechanism for low p_T particle production we propose a procedure to systematically quantify the model-to-data differences using modern Bayesian inference analysis techniques, which allows for consistent treatment of the experimental uncertainties. In this paper, we deploy a procedure based on the relativistic fluid framework F_LUIDuM with partial chemical equilibrium (PCE) coupled to T_RUIDuM initial state and F_ASTR_ESO decays to analyse p_T distributions of charged pions, kaons and protons measured in collisions of Pb–Pb at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ [39], Xe–Xe at $\sqrt{s_{\text{NN}}}=5.44\text{ TeV}$ [41], and Au–Au at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [43, 46]. Our results demonstrate the power of the proposed procedure to exploit the precision of the current experimental data in the search for limitations and improvements in the available state-of-the-art hydrodynamic model frameworks.

Modelling of heavy-ion collisions

Our model for simulating high-energy nuclear collisions combines three distinct components. The T_RENT_O model [57] was utilized for the initial conditions, while the F_LUIDuM model with a PCE implementation [58], featuring a mode splitting technique for fast computations, was used for the relativistic fluid dynamic expansion with viscosity. Additionally, the F_ASTR_ESO code [15] was used to take resonance decays into account.

The T_RENTo model involves positioning nucleons with a Gaussian width w using a fluctuating Glauber model, while ensuring a minimum distance d between them. Each nucleon contains m randomly placed constituents with a Gaussian width of v. T_RENTo uses an entropy deposition parameter p that interpolates among qualitatively different physical mechanisms for entropy production [57]. Furthermore, additional multiplicity fluctuations are introduced by multiplying the density of each nucleon by random weights sampled from a gamma distribution with unit mean and shape parameter k. For this study, the T_RENTo parameters are set based on [59]. The inelastic nucleon–nucleon cross sections are taken from the measurements by the ALICE and PHENIX Collaborations [60, 61]. The Pb and Au ions are sampled from a spherically symmetric Woods–Saxon distribution, while the Xe ion comes from a spheroidal Woods–Saxon distribution with deformation parameters β₂=0.21 and β₄=0.0 [62].

The software package F_LUIDuM [58], which utilizes a theoretical framework based on relativistic fluid dynamics with mode expansion [63-65], is used to solve the equations of motion for relativistic fluids. The causal equations of motion are obtained from second-order Israel–Stewart hydrodynamics [66]. As in our previous work [8], we are interested in examining the azimuthally averaged transverse momentum spectra of identified particles at midrapidity. Therefore, we do not consider azimuthal and rapidity-dependent perturbations and only require the background solution to the fluid evolution equations, neglecting terms of quadratic or higher order in perturbation amplitudes.

The Cooper–Frye procedure is used to convert fluid fields to the spectrum of hadron species on a freeze-out surface, which in our work is assumed to be a surface of constant temperature [67]. As in our previous work [8], the hadronic phase, after the chemical freeze-out and before the kinetic freeze-out, i.e. T_kin < T < T_chem, is modeled by a concept of partial chemical equilibrium (PCE), which replaces the need for a hadronic after-burner in the simulation. Our description follows the work described in Refs. [68-70], in which different particle species in a hadronic gas are treated as being in chemical equilibrium with each other, while the overall gas is not. During the PCE, the mean free time for elastic collisions is still smaller than the characteristic expansion time of the expanding fireball, thereby keeping the gas in a state of local kinetic equilibrium. The chemical equilibrium is not maintained if the mean free path of the inelastic collisions exceeds this threshold. On the kinetic freeze-out surface, we take the particle distribution function to be given by the equilibrium Bose–Einstein or Fermi–Dirac distribution (depending on the species), modified by additional corrections due to bulk and shear viscous dissipation [71, 72] and decays of unstable resonances [15]. We use a list of approximately 700 resonances from Refs. [73-75].

As described in Ref. [8], our central framework revolves around certain free parameters: the overall normalization constant Norm, ${\left(\eta /s\right)}_{\text{min}}$ and ${\left(\zeta /s\right)}_{\text{max}}$ in the shear and bulk viscosity to entropy ratio parametrizations, the initial fluid time τ₀, and the two freeze-out temperatures T_kin and T_chem. With our Bayesian inference analysis, we simultaneously determine these six model parameters within predefined intervals (refer to Tab. 1). These intervals are based on physical considerations and knowledge from previous studies [8, 12, 39, 48, 76]. It is worth mentioning that we have confirmed a posteriori that the optimal values fall within these intervals rather than on their boundaries, and in cases where no clear convergence was obtained, larger intervals were employed. Although F_LUIDuM is recognized for its fast execution speeds, the extensive parameter exploration involved in Bayesian analyses necessitates an approach to speed up the simulations. Our approach is based on the usage of an ensemble of artificial neural networks (ANNs) to emulate our model calculations. The training necessitates large datasets to achieve the required accuracy for replacing the simulation outputs. For each collision system, we use the outputs of ten thousand complete model calculations, with parameters distributed within the ranges presented in Tab. 1. The parameter values are generated using Latin hypercube sampling, which ensures a uniform density. With this large population of initial points in the parameter space, the emulator uncertainties result in a few percent. We refer readers to Refs. [8, 77], which provide extensive discussions. The posterior density is inferred from a probabilistic model, we use the numerical Markov-Chain Monte Carlo (MCMC) method [78], which is an efficient approach for exploring the probability space. Without clear guidance on how to precisely handle the degree of correlation in the experimental systematic uncertainties, the Bayesian inference analysis is performed assuming the experimental systematic uncertainties uncorrelated among the different particle species and transverse momentum intervals.

Determination of the optimal fitting range

The pion, kaon, and proton p_T spectra across various collision centrality classes measured by the ALICE Collaboration at the LHC and by the PHENIX and STAR Collaborations at RHIC in different colliding systems and center-of-mass energies, namely Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ [39], Xe–Xe collisions at $\sqrt{s_{\text{NN}}}=5.44\text{ TeV}$ [41], and Au–Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [43, 46] are used in this work. For the RHIC energy, the $\overline{\text{p}}$ spectra are used in the Bayesian inference analysis because in our model the stopping of the baryons from the colliding nuclei (baryon transport at midrapidity) is not included. Differently from our previous work [8], we expanded the Bayesian inference analysis to include multiple centrality intervals covering the range 0–40% for all collision systems. We highlight that our objective in this work is to systematically quantify the low-p_T pion excess and examine its possible dependence on collision centrality, collision energy, and colliding nuclei, rather than constraining physical parameters of the QGP across these three collision systems. For this reason, we have run the full framework separately for each centrality interval and collision system, without attempting to perform a global fit using all available data.

To determine the optimal $p_{\text{T}}^{\pi }$ range for fitting experimental measurements, and consequently to compute the low-p_T pion excess, we performed a Bayesian inference analysis varying each time the p_T interval of the pion spectra, while keeping them for kaon and proton spectra fixed ($p_{\text{T}}^{\text{K}\mathrm{,}p}<2.0$ GeV/c). Because of the large difference in masses, simultaneous inference of pion, kaon, and proton spectra was employed to achieve convergence of the model parameters. Initially, we optimized the starting $p_{\text{T}}^{\pi }$ within the range $x_1<p_{\text{T}}^{\pi }<2.0$ GeV/c, where x₁ ranged from 0.1 to 1.0 GeV/c. This range was chosen to ensure an adequate number of p_T intervals for the Bayesian inference procedure, with the upper limit for x₁ set at 1.0 GeV/c. Subsequently, we optimized the ending p_T by fitting within $0.5<p_{\text{T}}^{\pi }<x_2$ GeV/c, where x₂ ranged from 2.0 to 3.0 GeV/c.

Although the constraint of the QGP physical parameters is not the main focus of this research, it is crucial to monitor their performance and convergence while optimizing the pion p_T range. This ensures that the chosen p_T range in the Bayesian inference procedure leads to convergence. In Fig. 1 the six key parameters for the 0–5% centrality class in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ are shown. The Norm and ${\tau }_0$ parameters are depicted in a ratio format ($\text{Norm}/{\tau }_0$) because in our model the expected entropy density profile is obtained using their ratio [8, 48]. The top panel corresponds to the starting $p_{\text{T}}^{\pi }$ optimization procedure, while the bottom panel focuses on ending $p_{\text{T}}^{\pi }$ optimization. The values reported represent the median of the marginalized Bayesian posterior distributions for each model parameter, while error bars denote the 68% confidence interval.

Fig. 1Full-screen View

Parameter values within different p_T fitting ranges for the 0–5% centrality class in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{TeV}$. The upper panel displays the variation of the starting $p_{\text{T}}^{\pi }$ from 0.1 GeV/c to 1.0 GeV/c. Conversely, the lower panel illustrates the variation of the ending $p_{\text{T}}^{\pi }$ from 2.0 GeV/c to 3.0 GeV/c. The error bars in the figure denote 68% confidence intervals of the marginalized Bayesian posterior distributions for each model parameter

All parameters exhibit varying stability across the fitting ranges. In the top panel, the parameters converge to stable values when x₁ exceeds the p_T threshold of approximately 0.5 GeV/c. This indicates that using the low- p_T pion region (x₁ < 0.5 GeV/c) in the Bayesian inference procedure would introduce instabilities in constraining the physical parameters, demonstrating that a fluid dynamic framework cannot capture the experimentally measured low-p_T pion spectra. In the lower panel of Fig. 1, the parameters start deviating from the converged values again when the Bayesian procedure includes the spectra values for $p_{\text{T}}^{\pi }>2.0$ GeV/c. As we move to higher p_T, it is anticipated that particles are no longer predominantly produced thermally. We attribute this to the limit of the applicability of the fluid dynamics description at high p_T (emerging contribution from hard processes) and concluded that our results for the $p_{\text{T}}<2.5$ GeV/c are stable and not subject to overfitting or overtraining. Instead, contributions from hard partonic scattering processes become more pronounced, and the effects of partonic energy loss begin to dominate the spectral shape. On top of the parameter instabilities, it was observed that even when either the low-p_T or high- p_T spectra are included in the Bayesian procedure, the F_LUIDuM calculations fail to replicate the experimental data accurately. This results in significant discrepancies between the data and the model observed both at low and high p_T. The same study was also conducted for the 30–40% centrality interval in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ to verify the consistency of the findings. Comparable performances were observed across the different centrality intervals analyzed. As a result, the optimal p_T range with respect to a fluid dynamic description was established to be $0.5<p_{\text{T}}^{\pi }<2.0$ GeV/c for all centrality intervals and collision systems.

In Fig. 2 the Bayesian posterior distributions of the model input parameters utilized in this analysis for all centrality classes in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ are reported. It is important to note that we have confirmed the posterior distributions fall within the prior interval specified in Table 1, rather than on its boundaries. This marginal distribution plot illustrates that the parameters exhibit a high degree of consistency across the various centrality classes, lying within one standard deviation. Notably, the parameter T_kin shows a systematic shift in its median value towards more peripheral collisions. This observation aligns with previous findings obtained with the usage of a Blastwave fit [2], and it can be interpreted as a possible indication of a more rapid expansion towards central collisions and with the expectation of a shorter-lived fireball with stronger radial pressure gradients in more peripheral collisions. As discussed in [8], the (η / s)_min remains unconstrained, which we attributed to the limited sensitivity of the current observables to the shear viscosity of the system.

Fig. 2Full-screen View

(Color online) Marginal posterior distributions of the model input parameters for the five analyzed centrality classes in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{TeV}$

This study determined that the optimal p_T range for the pion p_T spectra is $0.5<p_{\text{T}}^{\pi }<2.0$ GeV/c. This range is recommended when similar Bayesian analyses are performed to constrain physical QGP parameters and use p_T differential pion variables. Having established this p_T interval, the one required to quantify the low- p_T pion excess is consequently defined as $0.1<p_{\text{T}}^{\pi }<0.5$ GeV/c for the LHC and $0.2<p_{\text{T}}^{\pi }<0.5$ GeV/c for RHIC.

Low-p_T pion yields

In Fig. 3 the ratios of the experimental spectra over model calculations for pions, kaons, and protons are shown. For those ratios, the model calculations are performed using the Maximum a Posteriori (MAP) estimates of the parameters [57]. The MAP estimate refers to the set of model parameters corresponding to the mode of the posterior distribution, representing the point in parameter space with the highest posterior probability. Given that we use uniform priors in our Bayesian inference, the MAP values are equivalent to those that maximize the likelihood function. The ratios are arranged in rows per particle type and columns per collision system. The bands depict experimental statistical and systematic uncertainties combined in quadrature. F_LUIDuM calculations yield nearly flat data-over-model ratios compatible with unity within one standard deviation across the entire p_T spectrum for pion, kaons, and protons for all centrality intervals and collisions systems in the intervals used in the Bayesian analysis.

Fig. 3Full-screen View

(Color online) Differential yields of pions (π), kaons (K), and protons (p) over model spectra for 0–40% centrality classes in Pb–Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ [39], Xe–Xe collisions at $\sqrt{s_{\text{NN}}}=5.44\text{ TeV}$ [41], and Au–Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [43], respectively. From top to bottom, each row corresponds to particle type, and from left to right, each column showcases one of the three collision systems under consideration. Within every panel, ratios are segmented into five (six) centrality classes ranging from 0% to 40%. The bands represent the statistical and systematic uncertainties of the data, summed in quadrature

We recall that for kaons and protons the p_T interval used for the Bayesian analysis is p_T < 2.0 GeV/c, while for pions it is 0.5 < p_T < 2.0 GeV/c. For p_T > 2.0 GeV/c, the model calculations for all hadrons start to deviate from experimental measurements, suggesting that the higher p_T domain may not be predominantly governed by soft processes, which can typically be described by fluid dynamic calculations. The observed deviations are larger for pions with respect to heavier particles, supporting the idea that hadrons originate from a fluid with a unified velocity field.

In the low-p_T range (p_T < 0.5 GeV/c), unlike kaons and protons, the data-over-model ratios for pions exceed unity across all centrality classes and collision systems, indicating a systematic pion production excess in the experimental measurements with respect to the fluid dynamic production. As discussed in the previous section, even when including the pion spectra in the Bayesian inference analysis, the fluid dynamic calculation is not able to capture this low-p_T interval.

In Fig. 4, the pion excess, computed as the difference between the integral of the experimentally measured pion spectra in the interval p_T < 0.5 GeV/c and the integral of the pions computed within our framework in the same p_T interval, is shown for the three collision systems as a function of centrality. It is important to notice that the excess is computed in two different p_T intervals for LHC and RHIC. At the LHC pion spectra are measured down to p_T = 0.1 GeV/c while at RHIC down to p_T = 0.2 GeV/c. This study focuses on the single-charge π excess, utilizing π⁺ p_T spectra in Pb–Pb at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ [39] and averaging π⁺ and π^- in Xe–Xe at $\sqrt{s_{\text{NN}}}=5.44\text{ TeV}$ [41]. As for the Au–Au system, the π^- p_T spectra from PHENIX [43] are used to compute the excess. For completeness, the pion excess is computed also utilizing the measured pion spectra from the STAR Collaboration (green markers) measured in Au–Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [46]. Due to the limited p_T interval of the STAR measurement (0.2–0.75 GeV/c), it was not possible to perform an independent Bayesian analysis, and the pion excess is computed using the model calculation obtained by the inference analysis of the PHENIX data. The excesses obtained using experimental data from the two Collaborations are compatible within the uncertainties. The STAR measurements are available for the 10–20% interval, while the PHENIX data for the centrality intervals 10–15% and 15–20%. Despite different acceptance in rapidity of the STAR, PHENIX, and PHOBOS experiments at RHIC, all measurements are reported at midrapidity per unit of rapidity and no additional treatment is required when comparing our calculations to these data from different RHIC experiments. Therefore, to calculate the pion excess, the fluid calculations from the two finer centrality classes are averaged into a larger one. The experimental uncertainties from measurements are reported as bars. For consistency with the treatment of the systematic uncertainties in the fit, while computing the pion excess, the experimental systematic uncertainties are propagated as fully uncorrelated across p_T. The total uncertainties represent the quadratic sum of experimental statistical and systematic uncertainties. A decreasing trend in the excess from central to peripheral collisions is observed for all collision systems. The significance of the excess is above 5 for all centrality classes and collision systems, specifically varying from 9.3 to 11.1 across centrality classes at the LHC energies. We estimated the effect of treating the experimental systematic uncertainties as partially or fully correlated and the extracted pion excess remained compatible with our main result reported in the paper. In the future, it will be beneficial to have experimental guidance on the degree of correlation of the uncertainties among p_T intervals and particle species of the experimental observables, which would enable a more thorough treatment of the systematic uncertainties in the analysis.

Fig. 4Full-screen View

(Color online) The integrated absolute single-charge pion excess in the range $0.1(0.2)<p_{\text{T}}<0.5$ GeV/c as a function of centrality in different collision systems. The bars represent the experimental uncertainties and the bands represent the model uncertainties

The uncertainties depicted by the bands in the figures originate from our model reflecting different sources of model and experimental uncertainty. The uncertainties represent the spread in posterior distributions and the extrapolation in the parameter space performed by the neural network (NN) emulator. As illustrated in Fig. 3, the MAP parameters provide a sufficient description of the data, indicating that the predominant source of model uncertainty arises from our NN emulator. Enhancements in the posterior distribution’s precision could be achieved by conducting the calibration with an increased number of design points and a narrower range of parameter values to increase the density of the training points to reduce the interpolation uncertainty.

In Fig. 5, the excess relative to the integral of the experimental data in the interval 0.1 < p_T < 2.0 GeV/c for the LHC and 0.2 < p_T < 2.0 GeV/c for RHIC is presented. When calculating the relative excess the systematic uncertainties between the excess and the integrated pion yields are treated as correlated, which partially cancels them out. For the STAR measurements, due to the limited p_T intervals, the PHENIX measurements are used as the denominator, hence no cancellation in the systematic uncertainties was possible. To obtain the yields for the 10–20% centrality interval from PHENIX, the arithmetic average of the yields for 10–15% and 15–20% was used. The relative excess remains constant as a function of centrality, with a consistent 10–20% excess across different collision systems. The computed relative excess indicates that fluid dynamic calculations account only for 80–90% of the measured pion production in heavy-ion collisions. An excess yield is found in all collision systems and centrality ranges considered.

Fig. 5Full-screen View

(Color online) The excess of pions normalized to the integrated pion yields over $0.1(0.2)<p_{\text{T}}<2.0$ GeV/c for various centrality classes in the different collision systems. The bars represent the experimental uncertainties and the bands represent the model uncertainties

Having performed the Bayesian inference analysis utilizing the available RHIC data at $\sqrt{s_{\text{NN}}}=20\text{ GeV}$, we can further compare our results with the measurements from the PHOBOS experiment at very low-p_T, to determine if a pion enhancement can also be quantified for p_T < 100 MeV [54]. The PHOBOS measurements are reported for the 0–15% centrality interval close to midrapidity (-0.1 < y < 0.4). In Fig. 6 the pion, kaon, and proton invariant yields are reported for the PHOBOS measurement at low-p_T (solid markers) and for the PHENIX data at larger p_T [43] (open markers), in comparison with the MAP model calculation obtained via the Bayesian inference analysis previously described. To obtain the centrality interval 0–15%, the results of our model and the PHENIX data for 0–5%, 5–10%, and 10–15% were averaged. The experimental systematic uncertainties were propagated as correlated across the different centrality intervals. To compute the sum of positively and negatively charged particles in our model, which predicts them in the same amount, the positively charged particle spectra were scaled by the experimentally measured numbers reported in Table IX of Ref. [43]. This correction is significant only for the proton case due to the experimental difference in the measured proton and antiproton spectra. In the bottom panels of Fig. 6 the ratios to the various particle species are reported. No significant enhancement of kaons and protons is observed within the current experimental precision, while a deviation of about 50% is observed for the low-p_T pion. This might indicate that the pion excess below p_T < 0.1 GeV/c saturates and does not keep rising to larger values. However, to compute an integral of the full pion excess it is important to have experimental measurements covering the full p_T without having gaps within the measurement, which is envisioned by the proposed next-generation detector ALICE 3 at the LHC [56].

Fig. 6Full-screen View

The invariant yields of pions, kaons, and protons as functions of p_T measured by the PHOBOS [55] and the PHENIX [43] Collaborations compared to the fluid dynamic calculations

Summary

In summary, we propose a procedure to advance in understanding the mechanism for low p_T particle production in heavy-ion collisions, which allows to systematically quantify the model-to-data differences using modern Bayesian inference analysis techniques and consistently treat the experimental uncertainties. We deploy this procedure using relativistic fluid framework F_LUIDuM with PCE coupled to T_RENTo initial state and F_AST_RESO decays to analyse p_T distribution of charged pions, kaons, and protons measured in heavy-ion collisions at top RHIC and the LHC energies. Despite the limited information about the correlation among systematic uncertainties in the data, our results indicate a systematic excess of pions produced below $p_{\text{T}}\approx 500$ MeV/c for both RHIC and LHC data. Our results demonstrate the power of the proposed procedure to fully exploit the precision of the experimental data and search for limitations and improvements in the available state-of-the-art hydrodynamic model frameworks. Further work, which is beyond the scope of this paper, is required for the interpretation of the observed low-p_T pion excess in terms of transverse momentum and particle species yield redistribution of the thermally produced particles, which are not yet modeled by the current state-of-the-art hydrodynamic models, or as evidence for a new particle production mechanism.