Extracting the speed of sound in quark–gluon plasma with ultrarelativistic lead–lead collisions at the LHC

Ultrarelativistic nuclear collisions create a strongly interacting state of hot and dense quark–gluon matter that exhibits a remarkable collective flow behavior with minimal viscous dissipation. To gain deeper insights into its intrinsic nature and fundamental degrees of freedom, we determine the speed of sound in an extended volume of quark–gluon plasma using lead–lead (PbPb) collisions at a center-of-mass energy per nucleon pair of 5.02 TeV. The data were recorded by the CMS experiment at the CERN LHC and correspond to an integrated luminosity of 0.607 nb − 1 . The measurement is performed by studying the multiplicity dependence of the average transverse momentum of charged particles emitted in head-on PbPb collisions. Our findings reveal that the speed of sound in this matter is nearly half the speed of light, with a squared value of 0 . 241 ± 0 . 002 ( stat ) ± 0 . 016 ( syst ) in natural units. The effective medium temperature, estimated using the mean transverse momentum, is 219 ± 8 ( syst ) MeV. The measured squared speed of sound at this temperature aligns precisely with predictions from lattice quantum chromodynamic (QCD) calculations. This result provides a stringent constraint on the equation of state of the created medium and direct evidence for a deconfined QCD phase being attained in relativistic nuclear collisions.


Introduction
When heavy atomic nuclei collide at relativistic speeds, a transformation occurs, giving rise to an exotic state of matter with a temperature above several trillion kelvin and known as the quark-gluon plasma (QGP) [1][2][3][4].In this realm of extreme temperatures, quarks and gluons break free from their confined existence inside hadrons, traversing long distances (e.g., several fm) compared to the size of individual nucleons.The emergence of the QGP represents a fundamental prediction of quantum chromodynamics (QCD) [5,6], the theory that elucidates the nature of the strong force.More remarkably, this strongly interacting QGP matter is found to exhibit the characteristics of an almost "perfect liquid" with little frictional momentum dissipation [7][8][9][10].Its collective dynamics and macroscopic properties are well described by the principles of nearly ideal relativistic hydrodynamics.
The equation of state (EoS) reveals the underlying fundamental degrees of freedom of a substance and is an invaluable tool to infer how the substance will respond to changes in its energy density.In fluid-like environments, the study of sound modes arising from longitudinal compression provides a means to determine the corresponding speed of sound, denoted as c s .This parameter, whose square is defined as the rate of pressure P change in response to variations in energy density ε, c 2 s = dP/dε [11], plays a pivotal role in characterizing the nature of the medium under investigation and in constraining models of corresponding EoS.The exploration of the sound wave propagation in strongly correlated systems, ranging from neutron stars to ultracold atomic gases [12,13], has garnered significant interest in recent years.Various methodologies have been proposed to experimentally extract the speed of sound in a QGP fluid [14][15][16][17][18], offering a direct means to constrain the QCD EoS.Notably, constraints on the speed of sound in hot QCD matter have been inferred through a comparison of relativistic nuclear collision data with theoretical models within a Bayesian framework [15].Recently, an effort to directly extract c 2 s in the QGP phase was made by establishing a connection to an effective static, uniform fluid system [16].That work was based on only two independent measurements of the charged-particle multiplicity density and mean transverse momentum (p T ) in lead-lead (PbPb) collision data from the ALICE experiment at center-of-mass energies per nucleon pair √ s NN = 2.76 and 5.02 TeV, and yielded a value of c 2 s = 0.24 ± 0.04 in natural units at a temperature of 222 ± 9 MeV.This result is in line with lattice QCD predictions, albeit subject to significant experimental uncertainties.
To increase the precision by which the speed of sound can be determined, a new hydrodynamic probe was later proposed in Ref. [17] utilizing the multiplicity dependence of mean p T measurements at a fixed √ s NN .This innovative technique makes use of "ultra-central" collisions in which the ions overlap almost entirely, i.e., collide at a very small impact parameter (b).A conceptual representation of this probe is illustrated in Fig. 1.The impact parameter of a heavy ion collision determines the size of the nuclear overlap region (system size), which is strongly correlated with the energy and entropy deposited in the initial state and the number of emitted charged particles in the final state ("multiplicity", N ch ).As the impact parameter decreases and collisions become increasingly central, both the system size and deposited energy increase, while maintaining a nearly constant initial energy density and temperature.However, this trend reaches its limit when b → 0. In this case, the initial system size is limited by the sizes of the participating nuclei.For symmetric PbPb collisions, this would be the size of a Pb nucleus.More energy and entropy can still be deposited into the fixed volume through fluctuations in the number of interacting partons.By examining the response of the temperature T to the increasing entropy density s at b ∼ 0, the speed of sound can be extracted based on fundamental thermodynamic laws, Here, in terms of experimental observables, s is directly proportional to N ch , while the temperature T relates to the average transverse momentum (⟨p T ⟩) of emitted particles with respect to the beam axis [16].Full hydrodynamic simulations, such as those made possible using the TRAJECTUM model [19], have verified the above relationship, although there are features that are not captured, as will be discussed later.As the c 2 s value depends only on the relative variation in ⟨p T ⟩ and N ch , any global changes to the observables, such as an increase in the system entropy through hadronic resonance decays [20], will not affect the result.
In this paper, we present a precise determination of the speed of sound in QGP using ultracentral PbPb collision data at √ s NN = 5.02 TeV, collected in 2018 by the CMS experiment at the CERN LHC.By achieving a level of precision of several percent, comparable to theoretical uncertainties, our results serve as a robust benchmark for comparison with hydrodynamic simulations and lattice QCD calculations of the EoS.These comparisons provide the most stringent and direct constraints on the degrees of freedom attained by the medium created in these collisions.Tabulated results are provided in the HEPData record for this analysis [21].

The CMS detector
The CMS apparatus [22] is a multipurpose, nearly hermetic detector, designed to trigger on [23,24] and identify electrons, muons, photons, and hadrons [25][26][27].The initial triggering is done with the level-1 system, which uses customized hardware to make the rapid online decision whether or not to accept an event and deliver it to the second system, the High Level Trigger (HLT).The HLT uses a large CPU farm to perform optimized online event reconstruction and characterize an event.A global "particle-flow" algorithm [28] aims to reconstruct all individual particles in an event, combining information provided by the all-silicon pixel and strip tracker, and by the crystal electromagnetic and brass-scintillator hadron calorimeters, operating inside a 3.8 T superconducting solenoid, with data from the gas-ionization muon detectors embedded in the flux-return yoke outside the solenoid.Hadron forward (HF) calorimeters [29], made of steel and quartz fibers, extend the pseudorapidity (η = − ln(tan(θ/2)), where the polar angle θ is defined relative to the counterclockwise beam) coverage provided by the barrel and endcap detectors.Two zero-degree calorimeters (ZDCs) [30], made of quartz-fibers and plates embedded in tungsten absorbers, are used to detect neutrons from nuclear dissociation events.

Data samples, event reconstruction and selection
The data analyzed, before applying the selection described below, consist of 4.27×10 9 minimum bias events, corresponding to an integrated luminosity of 0.607 nb −1 .The minimum bias events are triggered by requiring total energy signals above readout thresholds, which are in the range 6-12 GeV, on both sides of the HF calorimeters [24].Beam-gas interactions and nonhadronic collisions are rejected by requiring the shapes of the clusters in the pixel tracker to be compatible with those expected from particles produced by a PbPb collision [31].The events are also required to have at least one reconstructed primary vertex associated with two or more tracks within a distance of 15 cm from the nominal interaction point along the beam axis.The primary vertex is selected as the one with the highest track multiplicity in the event.Events with concurrent interactions per bunch crossing contribute to about 0.5% of the full data sample and are rejected based on the correlation of total energy deposited in the HF and ZDC detectors, following the procedure used in Ref. [32].The collision centrality in PbPb events, i.e., the degree of overlap or impact parameter of the two colliding nuclei, is commonly determined by the total transverse energy deposit in both HF calorimeters, E HF T, sum [31].As the main focus of this work is on collisions at small impact parameters, we analyzed only the 10% of PbPb events that had the largest E HF T, sum .This class contains the ultra-central collision events of interest.To ease the computational load for high-multiplicity central PbPb collisions, track reconstruction for PbPb events is done in two iterations.The first iteration reconstructs tracks from signals ("hits") in the silicon pixel and strip tracker that are compatible with trajectories of particles with p T > 1.0 GeV, while the second iteration reconstructs tracks compatible with trajectories of particles with 0.3 < p T < 1.0 GeV using solely the pixel detector.In the analysis, the tracks have the additional selection requirement of |η| < 0.5 for the best tracking performance.More details on the track reconstruction and selection can be found in Ref. [33].The tracking efficiency (ε eff ) and misreconstruction rate (ε mis ) are evaluated using the HYDJET [34] event generator, together with a full GEANT4 [35] simulation of the CMS detector response.These factors are combined to obtain an overall correction factor, ε trk = ε eff /(1 − ε mis ), which is used to account for detector effects on the total number of reconstructed tracks.The ε trk factor is calibrated not only in terms of p T and η, but also as a function of the detector occupancy.The occupancy is estimated by the total number of clusters registered in the silicon pixel tracker N pixel , where a weak linear decline of ε trk by up to 7% over an increase of N pixel by 30% is ob-served.In the analysis, each track is assigned a weight of 1/ε trk (η, p T , N pixel ) to account for track reconstruction effects.

Measurement method
The main experimental observable of this analysis is the mean transverse momentum ⟨p T ⟩ of charged particles in an event as a function of N ch , where ⟨p T ⟩ and N ch are measured within the same η and p T ranges (otherwise, rapidity-dependent entropy fluctuations would lead to a reduced signal [17]).Charged particle p T spectra for p T > 0.3 GeV are measured for events in 50 GeV intervals of E HF T, sum from 3400 GeV to 5200 GeV, with tracking efficiency and misreconstruction effects corrected.To avoid any bias in estimating ⟨p T ⟩ and N ch , it is necessary to extrapolate the measured p T spectra to the full p T range.The resulting ⟨p T ⟩ values (mean of the p T spectra) from all E HF T, sum intervals are then plotted against the corresponding N ch values (integral of the p T spectra) to form the final observable.The E HF T, sum variable essentially serves as a centrality estimator to vary the initial medium entropy density and temperature.In particular, as the E HF T, sum values are obtained in a forward η range that does not overlap with the range used to measure the corresponding ⟨p T ⟩ and N ch values, potential biases are avoided.For example, hard processes originating early in the collision tend to fragment into large numbers of high-p T particles, yet these particles may not reflect an increase in the entropy and temperature of the QGP medium.
The extrapolation of the p T spectra to the full p T range is performed by fitting a Hagedorn function [36] to the measured p T spectra over the range of 0.4 < p T < 4.5 GeV in each E HF T, sum interval.This method is found to provide an excellent description of the data [37] and models (TRAJECTUM and HYDJET).The chosen p T range for the fitting is varied to the evaluate corresponding uncertainties.The fitted functions are then used to extrapolate the missing portions of the p T spectra in the low-p T region.
As the extraction of the speed of sound mainly depends on the relative variation of ⟨p T ⟩ with respect to N ch (see Eq. ( 1)), normalized quantities, ⟨p T ⟩ norm = ⟨p T ⟩/⟨p T ⟩ 0 and N norm ch = N ch /N 0 ch , are used as the primary observables, where the ⟨p T ⟩ 0 and N 0 ch represent the mean transverse momentum and charged-particle multiplicity in a reference event class.Here, the centrality range chosen for the reference event class only needs to be close to that used for the speed of sound determination, and 5% most central events (as determined by E HF T, sum and denoted "0-5%") is used.By normalizing both ⟨p T ⟩ and N ch by their values in the reference event class, most of the systematic uncertainties can be minimized.The ⟨p T ⟩ 0 and N 0 ch values obtained are found to be in good agreement with the ALICE results in the 0-5% centrality range [37,38].Figure 2 shows the event fraction distribution as a function of the normalized multiplicity.
To extract the speed of sound, the expression that describes ⟨p T ⟩ norm as a function of N norm ch is taken from Ref. [17], as where, and σ can be constrained by fitting the measured multiplicity distribution using the procedure described in Ref. [39].The multiplicity distribution at fixed values of b is modeled using a Gaussian function.Integrating over b gives a minimum bias multiplicity distribution which can be fitted to data.As shown in Fig. 2, this fit provides a good description of the data.The results of this fit can be used to estimate the Gaussian mean and width at b = 0, yielding N knee ch = 1.1 and σ = 0.0272 with negligible uncertainties.Using the extracted N knee ch and σ values, a fit to the measured ⟨p T ⟩ norm as a function of N norm ch is performed using Eq. ( 2), thereby extracting the speed of sound.In practice, we limit the fit to the very high-multiplicity region of N norm ch > 1.14, as will be discussed in detail later.
The dominant sources of systematic uncertainties for the measured ⟨p T ⟩ norm and N norm ch values originate from the tracking correction and the extrapolation to the full p T range.As mentioned earlier, using normalized quantities minimizes the majority of the systematic uncertainties.Systematic uncertainties are directly evaluated for the normalized quantities, as well as for ⟨p T ⟩ 0 and c 2 s .The tracking correction uncertainty is evaluated by varying the default track selections to a set of looser or tighter values.The maximum deviation with respect to the default results is taken as a systematic uncertainty, which is found to be ±0.01GeV in ⟨p T ⟩ 0 and ±0.002 in the fitted c 2 s value.The p T extrapolation uncertainty is estimated by varying the range of measured spectra fitted by the Hagedorn function to a lower limit of 0.3 or 0.5 GeV and an upper limit of 4 or 5 GeV.The resulting systematic uncertainty is found to be at most ±0.023GeV for ⟨p T ⟩ 0 and ±0.012 for the c 2 s value.Systematic uncertainties for c 2 s associated with the choice of the lower fit limit in N norm ch are estimated by varying the limit from 1.13 to 1.17, resulting in an uncertainty of ±0.010 in c 2 s .Total uncertainties are obtained by adding the various sources in quadrature.Systematic uncertainties for ⟨p T ⟩ norm are extracted point-bypoint as a function of N norm ch .

Results
The observed multiplicity dependence of the average transverse momentum, both normalized by their values in the 0-5% centrality class, is presented in Fig. 3, within the kinematic range of |η| < 0.5 and extrapolated to the full p T range in central PbPb events.Hydrodynamic simulations from the TRAJECTUM [19,40,41] and Gardim et al. [17] models are also shown for comparison.Both models use an EoS from lattice QCD calculations [42].The TRAJECTUM model is a computational framework to simulate the full evolution of heavy ion collisions, which includes the modeling of initial stages, a viscous hydrodynamic phase with transport coefficients, and a hadronic gas phase.Parameters of the TRAJECTUM model are constrained by a global Bayesian analysis of a variety of experimental observables [19], where the band shown corresponds to uncertainties within the allowed range of TRAJECTUM configuration parameters.The model of Gardim et al. [17], besides the hydrodynamic phase, also considers the preequilibrium dynamics and hadronic interactions after thermal freeze-out.No uncertainties are evaluated for this model as only a single set of model parameters is used.
The ⟨p T ⟩ norm value first shows a very weak declining trend toward a local minimum around N norm ch ∼ 1.05.At higher multiplicities, corresponding to ultra-central PbPb events, a steep rise is observed, which is consistent with the expected increase in temperature with entropy density, as schematically illustrated in Fig. 1.The observed trend, including the minimum around N norm ch ∼ 1.05, is qualitatively consistent with the prediction by the TRAJECTUM model.A slightly steeper rise at high multiplicities is observed for the TRAJECTUM simulation when compared with the data.This suggests that the speed of sound used in the model may be slightly larger than is found in the QGP.However, this difference is not significant within experimental and theoretical uncertainties.The model by Gardim et al. also predicts a rise of ⟨p T ⟩ norm at very high multiplicities, with a slope similar to that observed in the data.However, it shows a flat trend at lower multiplicities instead of the local minimum structure around N norm ch ∼ 1.05 as seen in the data and the TRAJECTUM model.The origin of the observed local minimum is not currently understood.
To directly extract the speed of sound, the multiplicity dependence of the ⟨p T ⟩ norm data in Fig. 3 is fitted by Eq. (2).Because the observed local minimum is not captured by the simplified model in Eq. ( 2), the fit is performed only in the high-multiplicity range with N norm ch > 1.14.The final result of the squared speed of sound is found to be c 2 s = 0.241 ± 0.002 (stat) ± 0.016 (syst) in natural units.The same fit is also performed to the prediction from the TRAJECTUM model, resulting in c 2 s = 0.283 ± 0.045, where the model uncertainty is again determined within the allowed parameter space constrained by a global Bayesian analysis [19].
To constrain the EoS, a simultaneous determination of c 2 s and its corresponding temperature 0.8 0.85 0.9 0.95  ).Bars and the red band correspond to statistical and systematic uncertainties, respectively.Hydrodynamic simulations from the TRAJECTUM model [19] and the model by Gardim et al. [17] are also shown for comparison.The dashed line is a fit to the data using Eq. ( 2) in the range of N ch /N 0 ch > 1.14.
is necessary.Based on the hydrodynamic simulations discussed in Refs.[16,17], the effective temperature (T eff ) of the QGP phase is found to be given approximately by ⟨p T ⟩/3, with T eff = ⟨p T ⟩/3.07 quoted in Ref. [16] based on a soft equation of state.While the scaling factor relating T eff to ⟨p T ⟩ can depend on specific model assumptions, the theoretical uncertainty in this value is believed to be small compared to the quoted experimental uncertainties, thereby having no impact on the main conclusions drawn in this paper.In essence, T eff represents the initial temperature that a uniform fluid at rest would have if it possessed the same amount of energy and entropy as the QGP fluid does when it reaches its freeze-out state, the point at which the quarks become bound into hadrons.Due to longitudinal expansion and cooling, the T eff value is generally lower than the initial temperature of the QGP fluid.Nevertheless, it still characterizes a temperature in the QGP phase, to which the extracted c 2 s value based on the final-state ⟨p T ⟩ and N ch corresponds.Possible effects of shear and bulk viscosity are investigated in Ref. [16] and found to not impact this framework, as the shear viscosity increases ⟨p T ⟩ by about the same amount that the bulk viscosity decreases it.The ⟨p T ⟩ 0 value is measured to be 658 ± 25 (syst) MeV, leading to a T eff value for the ultra-central PbPb data of 219 ± 8 (syst) MeV (it varies by at most 2% toward the very end of N ch distribution within the 0-5% centrality range).The statistical uncertainty is orders of magnitude smaller than the quoted systematic uncertainties.s and T eff , while statistical uncertainties are smaller than the marker size.Values extracted from the TRAJECTUM simulation [19] following the same fitting procedure as the data and from the earlier work [16] are presented as the other colored boxes.The curve shows the prediction of c 2 s as a function of T from lattice QCD calculations [6].The dashed line at the value of 1/3 corresponds to the upper limit for noninteracting, massless gas ("ideal gas") systems [42].s value extracted in Ref. [16], and lattice QCD predictions of the c 2 s value as a function of T [6].The new CMS data allow for an unprecedented level of precision in the experimental determination of the speed of sound in an extended volume of QGP matter.The results exhibit excellent agreement with the lattice QCD prediction, with comparable uncertainties.Thus, our findings provide compelling and direct evidence for the formation of a deconfined QCD phase at LHC energies.

Conclusion
In summary, this study presents a measurement with a new hydrodynamic probe in ultrarelativistic nuclear collisions that results in the most precise determination to date of the speed of sound in an extended volume of quark-gluon plasma matter.By determining the dependence of the average transverse momentum on the total multiplicity for charged particles in nearly head-on lead-lead collisions at a center-of-mass energy per nucleon pair of 5.02 TeV, a squared speed of sound of 0.241 ± 0.002 (stat) ± 0.016 (syst) in natural units is determined.The effective medium temperature, estimated using the mean transverse momentum, is 219 ± 8 (syst) MeV.The excellent agreement of lattice quantum chromodynamics predictions with the experimental results provides strong evidence for the existence of a deconfined phase of matter at extremely high temperatures.Bucharest, Romania 93 Now at an institute or an international laboratory covered by a cooperation agreement with CERN 94 Also at Texas A&M University at Qatar, Doha, Qatar 95 Also at Kyungpook National University, Daegu, Korea 96 Also at another institute or international laboratory covered by a cooperation agreement with CERN 97 Also at Universiteit Antwerpen, Antwerpen, Belgium 98 Also at Yerevan Physics Institute, Yerevan, Armenia 99 Also at Northeastern University, Boston, Massachusetts, USA 100 Also at Imperial College, London, United Kingdom 101 Also at Institute of Nuclear Physics of the Uzbekistan Academy of Sciences, Tashkent, Uzbekistan

Figure 1 :
Figure 1: Conceptual representation of temperature vs. entropy density from mid-central to ultra-central heavy ion collisions.

Figure 2 :
Figure 2: The event fraction distribution as a function of the charged-particle multiplicity, N ch , within the kinematic range of |η| < 0.5 and extrapolated to the full p T range, in PbPb collisions at √ s NN = 5.02 TeV.The N ch value is normalized by its value in the 0-5% centrality class (N 0 ch ).The curve represents a fit to the data using the S. Das et al. model[39].tiplicity distribution at b = 0, normalized by N 0 ch .In Fig.2, the N knee ch value corresponds to the vicinity of the location beyond which the knee-shaped distribution starts rapidly falling.For the region of N norm ch

Figure 3 :
Figure 3: The average transverse momentum of charged particles, ⟨p T ⟩, as a function of the charged-particle multiplicity, N ch , within the kinematic range of |η| < 0.5 and extrapolated to the full p T range in PbPb collisions at √ s NN = 5.02 TeV.Both ⟨p T ⟩ and N ch are normalized by their values in the 0-5% centrality class (⟨p T ⟩ 0 and N 0 ch).Bars and the red band correspond to statistical and systematic uncertainties, respectively.Hydrodynamic simulations from the TRAJECTUM model[19] and the model by Gardim et al.[17] are also shown for comparison.The dashed line is a fit to the data using Eq.(2) in the range of N ch /N 0 ch > 1.14.

Figure 4 :
Figure 4: The speed of sound, c 2 s , as a function of the effective temperature, T eff , with the CMS data point obtained from ultra-central PbPb collision data at √ s NN = 5.02 TeV.The size of the red box indicates systematic uncertainties of c 2s and T eff , while statistical uncertainties are smaller than the marker size.Values extracted from the TRAJECTUM simulation[19] following the same fitting procedure as the data and from the earlier work[16] are presented as the other colored boxes.The curve shows the prediction of c 2 s as a function of T from lattice QCD calculations[6].The dashed line at the value of 1/3 corresponds to the upper limit for noninteracting, massless gas ("ideal gas") systems[42].

Figure 4
Figure 4 depicts c 2 s as a function of T eff , with the CMS data point obtained from ultra-central PbPb collision data at √ s NN = 5.02 TeV.The results are compared to the TRAJECTUM model, the c 2s value extracted in Ref.[16], and lattice QCD predictions of the c 2 s value as a function of T[6].The new CMS data allow for an unprecedented level of precision in the experimental determination of the speed of sound in an extended volume of QGP matter.The results exhibit excellent agreement with the lattice QCD prediction, with comparable uncertainties.Thus, our findings provide compelling and direct evidence for the formation of a deconfined QCD phase at LHC energies.
Individuals have received support from the Marie-Curie program and the European Research Council and Horizon 2020 Grant, contract Nos.675440, 724704, 752730, 758316, 765710, 824093, and COST Action CA16108 (European Union); the Leventis Foundation; the Alfred P. Sloan Foundation; the Alexander von Humboldt Foundation; the Science Committee, project no.22rl-037 (Armenia); the Belgian Federal Science Policy Office; the Fonds pour la Formation à la Recherche dans l'Industrie et dans l'Agriculture (FRIA-Belgium); the Agentschap voor Innovatie door Wetenschap en Technologie (IWT-Belgium); the F.R.S.-FNRS and FWO (Belgium) under the "Excellence of Science -EOS" -be.h project n. 30820817; the Beijing Municipal Science & Technology Commission, No. Z191100007219010 and Fundamental Research Funds for the Central Universities (China); the Ministry of Education, Youth and Sports (MEYS) of the Czech Republic; the Shota Rustaveli National Science Foundation, grant FR-22-985 (Georgia); the Deutsche Forschungsgemeinschaft (DFG), under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306, and under project number 400140256 -GRK2497; the Hellenic Foundation for Research and Innovation (HFRI), Project Number 2288 (Greece); the Hungarian Academy of Sciences, the New National Excellence Program -ÚNKP, the NKFIH research grants K 124845, K 124850, K 128713, K 128786, K 129058, K 131991, K 133046, K 138136, K 143460, K 143477, 2020-2.2.1-ED-2021-00181, and TKP2021-NKTA-64 (Hungary); the Council of Science and Industrial Research, India; ICSC -National Research Center for High Performance Computing, Big Data and Quantum Computing, funded by the EU NexGeneration program (Italy); the Latvian Council of Science; the Ministry of Educa-tion and Science, project no.2022/WK/14, and the National Science Center, contracts Opus 2021/41/B/ST2/01369 and 2021/43/B/ST2/01552 (Poland); the Fundac ¸ão para a Ciência e a Tecnologia, grant CEECIND/01334/2018 (Portugal); the National Priorities Research Program by Qatar National Research Fund; MCIN/AEI/10.13039/501100011033,ERDF "a way of making Europe", and the Programa Estatal de Fomento de la Investigaci ón Científica y Técnica de Excelencia María de Maeztu, grant MDM-2017-0765 and Programa Severo Ochoa del Principado de Asturias (Spain); the Chulalongkorn Academic into Its 2nd Century Project Advancement Project, and the National Science, Research and Innovation Fund via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, grant B37G660013 (Thailand); the Kavli Foundation; the Nvidia Corporation; the SuperMicro Corporation; the Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).