Effects of saturation and fluctuating hotspots for flow observables in ultrarelativistic heavy-ion collisions

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Effects of saturation and fluctuating hotspots for flow observables in ultrarelativistic heavy-ion collisions © Authors. Published by the American Physical Society Published version Hirvonen, Henry; Kuha, Mikko; Auvinen, Jussi; Eskola, Kari J.; Kanakubo, Yuuka; Niemi, Harri Hirvonen, H., Kuha, M., Auvinen, J., Eskola, K. J., Kanakubo, Y., & Niemi, H. (2024). Effects of saturation and fluctuating hotspots for flow observables in ultrarelativistic heavy-ion collisions. Physical Review C, 110(3), Article 034911. https://doi.org/10.1103/physrevc.110.034911 2024 PHYSICAL REVIEW C 110, 034911 (2024) Effects of saturation and fluctuating hotspots for flow observables in ultrarelativistic heavy-ion collisions Henry Hirvonen, Mikko Kuha , Jussi Auvinen ,KariJ.Eskola , Yuuka Kanakubo , and Harri Niemi University of Jyväskylä, Department of Physics, P.O. Box 35, FI-40014 University of Jyväskylä, Finland and Helsinki Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland (Received 3 July 2024; accepted 29 August 2024; published 20 September 2024) We investigate the effects of saturation dynamics on midrapidity flow observables by adding fluctuating hotspots into the novel Monte Carlo EKRT (MC-EKRT) event generator for high-energy nuclear collisions. We demonstrate that the intensity of the saturation effects significantly affects the ratio between the flow coefficients v3and v2at the LHC. Adding a hotspot substructure to the nucleons enhances the saturation effects and improves the agreement with the measured data. We show that the collision-energy dependence of the flow coefficients obtained using the MC-EKRT initial states with hotspots is improved in comparison with the earlier event-by-event EKRT model. In addition, we present the results for the charged hadron multiplicity distribution in Pb+Pb collisions at the LHC, and show that the minijet multiplicity originating fluctuations of the saturation scale included in MC-EKRT, as well as the presence of hotspots, are necessary for describing the measured large-multiplicity tail in the distribution. DOI: 10.1103/PhysRevC.110.034911 I. INTRODUCTION The highest-energy nucleus-nucleus collisions, ultrarelativistic heavy-ion collisions, which are currently performed at the CERN Large Hadron Collider (LHC) and at the Brookhaven National Laboratory (BNL) Relativistic Heavy Ion Collider (RHIC), aim at determining the properties of the nearly net-baryon-free hot quark-gluon plasma (QGP). One also strives for a detailed understanding of the stronginteraction dynamics that is responsible for the creation and further evolution of the QGP in these collisions. See, e.g., Ref. [1]forareview. According to lattice simulations of quantum chromodynamics (QCD), the theory of the strong interaction, the strongly interacting matter takes the form of the QGP at high temperatures of T150–160 MeV [2,3] at a vanishing baryochemical potential. Quarks and gluons can be produced in ultrarelativistic heavy-ion collisions from the kinetic energy of the colliding nuclei so copiously that the effective temperature (energy over particle ratio) of the system clearly exceeds 160 MeV. In these conditions, the normal formation of the color-confined, color-singlet bound states, hadrons, is momentarily inhibited, and a nearly thermalized QGP, where the degrees of freedom are colored gluons, quarks, and antiquarks, can be formed. The subsequent spacetime evolution stages of such a QCD matter—the expansion and cooling of Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. the QGP, the crossover transition to a hadron gas, followed by the expansion and cooling of the hadron gas—as well as the simultaneous appearance of the QGP and hadron-gas phases in different density regions of the expanding system, are describable in terms of relativistic dissipative fluid dynamics [4–27]. While QCD is a cornerstone of the standard model of particle physics, relativistic fluid dynamics has become a standard tool in the analysis of heavy-ion observables. The determination of the QCD matter properties, such as its equation of state and transport properties such as the shear and bulk viscosities, from the measured LHC and RHIC observables is a highly challenging task. Clearly, a precise determination requires a simultaneous analysis of as many heavy-ion observables as possible, from as many collision systems and collision energies as possible—a global analysis of heavy-ion observables [11,14,17,28]. A proper statistical analysis, Bayesian inference [26,29–40] is necessary for setting well-defined uncertainties to the extracted matter properties. Interestingly, neural networks are currently making it possible to include also statistics-expensive observables, such as complicated rare flow correlators, into the global analysis [41,42] (see also Ref. [43]). The mentioned global analyses of heavy-ion observables are based on a fluid-dynamical description, which takes initial densities and flow velocities of the produced QCD matter as initial conditions. One either parametrizes these initial conditions [26,29–33,35–40,44] or tries to compute them from a QCD dynamical model for the initial production of gluons and quarks [14,17,28,38]. In both cases there is some number of fit parameters that characterize the initial states, and these will obviously be correlated with the actual QCD-matter properties extracted from the data via Bayesian inference. It is therefore important to model the QCD-matter initial states based on 2469-9985/2024/110(3)/034911(11) 034911-1 Published by the American Physical Society HENRY HIRVONEN et al. PHYSICAL REVIEW C 110, 034911 (2024) QCD dynamics as far as is possible, in order to understand the dominant particle production mechanism, to reduce the uncertainties in the extraction of the initial states, and to have predictive power for moving from one system to another. The EKRT (Eskola-Kajantie-Ruuskanen-Tuominen) model [17,45–47], which treats the nuclear collisions as collisions of parton clouds, and supplements a perturbative QCD (pQCD) calculation for the production of few-GeV partons (minijets) [48,49] with a collinear factorization-inspired QCD saturation mechanism [17,47] for regulating the small-pT minijet production (pTbeing transverse momentum), is an example of such a QCD-based initial state modeling with predictive power. The event-by-event (EbyE) version of the model, EbyE-EKRT [17], has been quite successful in explaining a large collection of heavy-ion bulk observables at the LHC and RHIC [17,28,34,50–52]. The latest progress here is the novel Monte Carlo (MC)-EKRT event generator, introduced recently in Ref. [53], and employed in the present paper. The new features in MC-EKRT [53] relative to EbyE- EKRT [17] are that now the produced partonic system contains local fluctuations of the minijet multiplicity, which in turn induce dynamical fluctuations to the saturation controlling the initial parton production. Also per-nucleon conservation of energy and valence-quark numbers are accounted for. MC-EKRT also introduces a new type of spatially dependent nuclear parton distribution functions (snPDFs) that are specific to the nucleon configuration in each event and can cope with the largest density fluctuations of the nucleon densities. Thanks to these new features, MC-EKRT gives initial conditions for full (3+1)-dimensional [(3+1)D] EbyE fluid dynamics, and thus enables the studies of rapidity-dependent observables, such as rapidity distributions of yields and flow coefficients of charged hadrons in Pb+Pb collisions at the LHC and at the highest-energy Au+Au collisions at RHIC, see Ref. [53]. In this paper, we employ the new MC-EKRT framework for computing event-by-event initial conditions for (2+1)D dissipative shear- and bulk-viscous second-order transient fluid dynamics in the midrapidity unit of 5.02 and 2.76 TeV Pb+Pb collisions at the LHC. In particular, we study the sensitivity of the flow coefficients vnto the model details, such as the nucleonic width and substructure, the Gaussian smearing in coupling the individual minijets to continuous fluid dynamics, as well as the order in which we do the minijet filtering based on saturation and conservation of energy. In addition, we show how the added minijet multiplicity fluctuations are the piece formerly missing from EbyE-EKRT in explaining the behavior of the charged multiplicity distributions in the most central collisions. The recently developed neural networks for predicting flow observables directly from the initial energy density event by event [41,42], are also utilized. As the main result of this paper, we show that a detailed simultaneous description of the vn’s requires saturation to be the driving QCD mechanism for initial parton production. In particular, this result calls for further nucleonic substructure, hotspots, to be introduced in MC-EKRT. We also implement these in MC-EKRT and discuss their interesting interplay with saturation, in describing the v2/v3ratio as well as in explaining the measured charged multiplicity distributions. II. MC-EKRT INITIAL STATE FOR FLUID DYNAMICS A. Minijet sampling The MC-EKRT event generator of Ref. [53] produces partonic initial states, i.e., saturated systems of gluons and quarks with pTp0∼1 GeV, that can be fed as initial conditions to (3+1)D event-by-event fluid-dynamical simulations. The generation of such MC-EKRT initial states proceeds via the following steps (for details, see Ref. [53]). First, the nucleon configurations of the colliding (here spherically symmetric) nuclei Aand Bare generated by sampling the standard two-parameter Woods-Saxon distribution, and by requiring an exclusion radius of 0.4 fm. A squared impact parameter b2 AB for the A+Bcollision, defining the distance between the centers of masses of the colliding nuclei, is sampled from a uniform distribution. In the absence of hotspots (i.e., without subnucleonic density fluctuations), the A+Bcollision is triggered using MC Glauberlike blackdisc nucleons with a trigger cross section identical to the inelastic nucleon-nucleon cross section σNN inel , which is obtained from the measured total and elastic nucleon-nucleon cross sections as a function of the nucleon-nucleon center-of- momentum system (CMS) energy √sNN [54,55]. Once the A+Bcollision is triggered, MC-EKRT does not consider nucleonic subcollisions at all but pictures the entire nuclear collision as a collision of two extensive parton clouds. For distributing the parton subclouds spatially around each nucleon, MC-EKRT assumes a Gaussian thickness function, TN(¯s)=1 2πσ2 N exp −|¯s|2 2σ2 N,(1) with a width parameter σN=σN(√sNN) that is obtained from exclusive photoproduction of J/ in photon-proton collisions at HERA [56,57]. Then, multiple dijet production, i.e., the number of independent dijets with jet transverse momentum pT⩾p0=1 GeV, that is assigned to originate from each ab pair, is sampled from a Poissonian probability distribution withamean ¯ Nab jets =TNN(¯ bab)σab jet (p0,√sNN,{¯sa},{¯sb}),(2) where TNN(¯ bab) is the nucleonic overlap function and ¯ bab is the impact parameter between the nucleons aand b, while σab jet is the integrated pQCD (mini)jet cross section, which MCEKRT computes using the novel snPDFs for aand b, and all possible leading-order (LO) partonic 2 →2 subprocesses. A CMS-energy-dependent multiplicative Kfactor is introduced to σab jet as a free fit parameter, to account for the missing higherorder contributions. The (mini)jet cross section depends on the transverse momentum cutoff parameter p0,ontheCMS energy √sNN, as well as on the transverse locations ¯saand ¯sb of aand bin the nucleon configurations of Aand B, indicated here with {¯sa}and {¯sb}. As explained in detail Ref. [53], the novel snPDFs are now nucleon-configuration specific and account for the nuclear 034911-2 EFFECTS OF SATURATION AND FLUCTUATING … PHYSICAL REVIEW C 110, 034911 (2024) modifications of each nucleon’s PDFs caused by all other nucleons in the nucleus. In other words, the MC-EKRT snPDFs are nucleon specific and nucleon-configuration specific. Also noteworthy is that these novel snPDFs can fully cope with the event-by-event density fluctuations, which was not the case with the formerly developed spatial nPDFs, such as those in Ref. [58]. The MC-EKRT snPDFs are normalized (averaging over all nucleons in each nucleus and over a large number of nuclei) to the spatially averaged nuclear PDF modifications of the EPS09LO set [59], and CT14LO [60] are employed for the free proton PDFs. Finally, the transverse location for each produced dijet is sampled from the product of the two overlap functions TN, whose transverse integral gives the usual overlap function TNN. The kinematic variables and the flavor chemistry of the produced partons, along with identifying the valence quark-consuming processes, is sampled from the differential jet sub-cross-sections, as explained in Ref. [53]. B. Minijet filtering The next, and decisive, step in MC-EKRT is the filtering of the excessive candidate dijets, based on the EKRT saturation [17,45–47] and conservation of energy and valence quark numbers. As explained in Refs. [17,47,53] saturation here is expected to occur when all the higher-order (n>2) →2 parton processes start to dominate over the 2 →2 ones. For maintaining collinear factorization at the highest values of jet transverse momenta, the filterings are performed in the order of decreasing factorization scale, which here is the jet pT. Then, the highest-pTpartons can remain in the system while the lower-pTones may get filtered away. For the saturation filtering, MC-EKRT assigns a transverse radius 1/(κsat pT) for each dijet candidate, where κsat is a packing factor, a free parameter to be fitted from the data. The transverse position of each candidate dijet is kept track of, and a candidate dijet gets filtered away if it overlaps with any of the previously accepted dijets. As shown in Ref. [53], after the saturation filtering the pTdistribution of surviving partons is not anymore sensitive to the original cutoff parameter p0but now saturation is the dynamical and local regulation mechanism for these distributions. This is the major difference to the traditional minijet eikonal models (and models alike), which are employed in event generators describing multiparton interactions, such as HIJING [61]. Similarly, MC-EKRT keeps track of all the longitudinal momentum fractions and valence quarks drawn out from their mother nucleons by the candidate dijets. If the candidate dijet would make its mother nucleon exceed its energy or valencequark budget, again checking the dijet candidates in the order of decreasing pT, then that dijet candidate gets filtered away. In the EKRT framework, in the spirit of suggesting saturation as the dominant QCD mechanism that regulates and controls initial parton production in highest-energy nuclear collisions, the default is to do the saturation filtering first, and only then the energy and valence-quark number conservation filterings. There is, however, an option in the code, which we utilize, and consequences we study in this paper, to have all the filterings done simultaneously. C. Nucleon substructure and hotspot trigger The fluctuating substructure to the nucleons of the MCEKRT framework is implemented as follows. While there is clear evidence that the nucleon substructure is necessary for describing the measured incoherent J/ψ photoproduction [62], the situation is less clear in heavy-ion collisions. The global analyses performed in Refs. [37,44] provide a slight preference towards the inclusion of the nucleon substructure, but the evidence is not conclusive. However, these analyses use the TRENTo [63] initial-state model, in which the effect of the substructure can partly be compensated with other initialstate parameters. In the MC-EKRT model, the addition of the nucleon substructure enhances the saturation effects since it confines the minijet production into more localized transverse regions. This leads to a change in the initial geometry, which might have an impact on the flow observables. The nucleon substructure is implemented by introducing Gaussian hotspots to the nucleon thickness function: TN(¯s)=1 Nh Nh  i=1 1 2πσ2 h exp −¯s−¯sh i 2 2σ2 h,(3) where Nhis the number of hotspots, and σhis the width of the hotspot. In this paper, Nh=3 is always used when the nucleon substructure is enabled. The hotspot locations ¯sh iare sampled from a two-dimensional Gaussian distribution with a width σs. The total nucleon width σNis then related to the hotspot widths via σ2 N=σ2 s+σ2 h. Therefore, only two of the three widths are independent. As in Refs. [53,57,64], the energy dependence of the total nucleon width is parametrized as σN=√bwith b/GeV−2=b0+4α Pln W W0,(4) where W=√sNN, and b0,α P, and W0are fit parameters. In the present paper, our default choice of parameters, based on the H1 measurements [64], are b0=4.63, α P=0.164, and W0=90 GeV. This corresponds to σN=0.517 fm for 2.76 TeV, and σN=0.532 fm for 5.023 TeV collision energies. In principle, the nucleon substructure needs to be accounted for when performing the triggering of the nuclear collision event [65,66] since otherwise there might be events where the collision is accepted even though there is no hadronic interaction. As mentioned before, without any substructure, the triggering is done by assuming hard-sphere scattering between two nucleons. The event is accepted if the distance dNN min between any nucleons a∈Aand b∈Bsatisfy dNN min <σNN inel π,(5) where σNN inel is the inelastic nucleon-nucleon cross section. The same kind of geometrical criterion can be extended to account for the locations of the hotspots. That is, the triggering with the nucleon substructure is done based on the minimum 034911-3 HENRY HIRVONEN et al. PHYSICAL REVIEW C 110, 034911 (2024) distance between two colliding hotspots dHS min, i.e., dHS min <σHS π,(6) where σHS is an effective hotspot-hotspot cross section fitted to reproduce the same nucleus-nucleus cross section as obtained with condition (5). Therefore, the value of σHS will depend on the hotspot sampling width σsand the collision system. Even though in principle hotspot triggering could have a notable impact, we have noticed that in most cases all the measured observables remain nearly unchanged in the 0–80 % centrality range. The largest effects are most visible in the most peripheral charged particle multiplicity region, where usually no measured data are given. In the 60–80 % centralities, the differences in charged particle multiplicities are only a few percent at most. However, since in MC-EKRT we sample dijets from the same nucleon configuration until at least one is produced in a collision, the addition of hotspot triggering there speeds up the generation of the initial states. D. Initialization of fluid dynamics The initial condition of fluid dynamics is the energymomentum tensor Tμν at some initial proper time τ0. However, MC-EKRT produces a list of massless partons with known momentum rapidities yi, transverse momenta pTi, and transverse coordinates x⊥,0i. Thus, the partons need to be propagated to the τ0surface and converted to the components of the energy-momentum tensor. Here we assume that all the partons are produced at the longitudinal location zi=0at time t=0, and that they propagate as free particles to the proper time τ0=0.2 fm. Therefore, spacetime and momentum rapidities are equivalent, i.e., ηs,i=yi. The spacetime coordinates of the parton iare then [τ0,x⊥i(τ0),η s,i] where x⊥i(τ0)=x⊥,0i+τ0pTi/pTi. The components of the energy-momentum tensor in the τ−ηscoordinates are obtained as in Ref. [53], Tαβ (xα)= id2pTdypαpβ pτ 1 τcosh(y−ηs) ×δ(2)(x⊥−x⊥i)δ(ηs−ηs,i)δ(2)(pT−pTi) ×δ(y−ηs),(7) where the four-momentum pα=(pτ,pT,pη) at a spacetime location xα=(τ,x⊥,η s) is given by pα=⎛ ⎝ pTcosh(y−ηs) pT τ−1pTsinh(y−ηs)⎞ ⎠.(8) Depositing all energy and momentum of a parton into a single cell on a hydro grid as suggested by the delta functions appearing in Eq. (7) would lead to extreme fluctuations in energy and momentum densities. To obtain smooth density distributions, smearing is required. Here we are performing (2+1)D hydrodynamic simulations, where a natural choice is to let all partons that are produced in the midrapidity window ycontribute to the fluid dynamical initial state. That is, in Eq. (7) we replace δ(ηs−ηs,i)→θ(y/2−|ηs,i|)/y, where θis the Heaviside theta function. Here we use y= 1.0, but we have tested that the final results are practically insensitive to the choice of yas long as 0.5⩽y⩽2.0. The smearing in the transverse (x,y) plane is performed by replacement δ(2)(x⊥−x⊥i)→g⊥(x⊥;x⊥i), where g⊥(x⊥;x⊥i)=C⊥ 2πσ2 ⊥ exp −(x⊥−x⊥i)2 2σ2 ⊥(9) is a Gaussian distribution with transverse smearing width σ⊥, which is normalized as d2x⊥g⊥(x⊥;x⊥i)=1.(10) The computation cost is reduced by imposing a ±3σ⊥cutoff on the smearing range, and the coefficient C⊥takes care of the unit normalization. As in Ref. [53], we only consider the local rest frame energy density ewhen initializing the fluid dynamical system, i.e., we neglect the initial transverse velocity and the initial components of the shear-stress tensor. Therefore, the initialization is determined by Tττ(τ0,x⊥,y)=1 τ0y i pTig⊥(x⊥;x⊥i)θ(y/2−|yi|), (11) which in this case coincides with e. The remaining components are then obtained, using the equation of state, as Tij = P(e)δij. Finally, we emphasize that even if we utilize only the midrapidity minijets in computing the above initial conditions, the underlying MC-EKRT event generation is fully 3D. Thus, the midrapidity initial conditions are influenced also by the finite-rapidity effects in saturation and in energy conservation. III. FLUID SIMULATION FRAMEWORK The simulations performed in this paper focus on midrapidity observables and therefore we assume that the longitudinal expansion of the system is boost invariant. The same framework as in Ref. [28] is used, i.e., we evolve the initially formed strongly interacting matter using dissipative fluid dynamics, and compute the final particle spectra at the dynamical decoupling surface. Additionally, the neural networks trained in Ref. [41] are utilized for significantly decreasing the computation time of the simulations. In this section, we give a brief recapitulation of each aspect of the framework. A. Fluid dynamics Fluid dynamics is based on the local conservation laws for energy, momentum, and conserved charges. Here we neglect the conserved charges, in which case the conservation law for the energy-momentum tensor, ∂μTμν =0, controls the dynamics. The energy-momentum tensor can be decomposed with respect to four-velocity uμas Tμν =euμuν−Pμν +πμν ,(12) 034911-4 EFFECTS OF SATURATION AND FLUCTUATING … PHYSICAL REVIEW C 110, 034911 (2024) where μν =gμν −uμuνis a projection operator, P= −1 3μνTμν is the total isotropic pressure, e=Tμν uμuνis the local rest frame energy density, and πμν =Tμνis the shear-stress tensor. The angular brackets denote the symmetric, traceless part of the tensor that is orthogonal to the fluid four-velocity. Here the fluid velocity is defined in the Landau frame, i.e., Tμ νuν=euμ. The bulk viscous pressure is defined as the deviation of the isotropic pressure Pfrom the equilibrium pressure P0, i.e., =P−P0. The equilibrium pressure is given by the equation of state (EoS) of the QCD matter at zero baryon density, P0=P0(e). In this work, we use the s95p-v1 parametrization [67] for the EoS, which includes the partial chemical decoupling at Tchem =155 MeV. The partial chemical decoupling is implemented by adding temperaturedependent chemical potentials for each hadron in the hadronic part of the EoS [68–70]. The conservation laws together with the EoS are enough to solve the evolution in equilibrium, but additional constraints are needed when dissipative effects are present. The dissipative parts of the energy-momentum tensor are the shear-stress tensor and the bulk viscous pressure. In the formalism by Israel and Stewart [71], the equations of motion for dissipative parts take a form τ d dτ+=−ζθ −δθ +λπ πμν σμν ,(13) τπ d dτπμν+πμν =2ησμν +2τππμ αωνα −δπππμν θ−τπππμ ασνα +ϕ7πμ απνα+λπσ μν,(14) where θ=∇μuμis the expansion rate, σμν =∇μuνis the strain-rate tensor, and ωμν =1 2(∇μuν−∇νuμ) is the vorticity tensor. The first-order transport coefficients ηand ζare called shear and bulk viscosity, respectively. In a 14-moment approximation to the massless gas [72–75], the first-order transport coefficients are related to the shear and bulk relaxation times as τπ=5η e+P0 ,τ =151 3−c2 s2 (e+P0)−1 ζ, (15) and the remaining second-order transport coefficients are δ =2 3τ,λ π =8 51 3−c2 sτ,δ ππ =4 3τπ, τππ =10 7τπ,ϕ 7=9 70P0 ,λ π =6 5τπ,(16) where csis the speed of sound. The specific shear viscosity η/sand specific bulk viscosity ζ/sare from the η/s=dyn parametrization introduced in Ref. [28]. B. Decoupling and particle spectra The fluid dynamic evolution is continued until reaching the kinetic decoupling surface. Here the decoupling surface is determined by the dynamical decoupling conditions Kn =τπθ=CKn (17) γτ π R=CR,(18) where Kn is the Knudsen number, γis the Lorentz gamma factor, and the coefficients CKn and CRare proportionality constants of O(1), which are fitted to the measured data. Here, values CKn =0.8 and CR=0.15 are used according to Ref. [28]. The size of the system Ris defined as R=A π,(19) where Ais the area in the transverse (x,y) plane where Kn < CKn. Additionally, the decoupling is forced to happen in the hadronic phase of the QCD matter, i.e., when T<150 MeV. Given these conditions the decoupling surface is determined using the Cornelius algorithm [76]. At the decoupling surface with the directed surface element dμ, the Lorentz-invariant particle spectrum for particle type iis computed according to the Cooper-Frye integral, Ed3Ni d3k= dμkμfi(x,k),(20) where Eand kμare particles energy and four-momentum, respectively. The distribution function for particle species i is decomposed into in- and out-of-equilibrium parts as fi= f0i+δfi, where the equilibrium part is given by f0i(x,k)=exp kμ iuμ−μi T±1−1 ,(21) where +(−) sign is for fermions (bosons), and μiis the chemical potential. Here, the viscous corrections to the equilibrium distribution are of the form [6,77–79] δfi=−f0i˜ f0i Cbulk Tm2 3Ek−1 3−c2 sEk +f0i˜ f0i 2T2(e+P0)πμνkμkν,(22) with ˜ f0i=1±f0i(+for bosons and −for fermions) and the coefficient 1 Cbulk = i gim2 i 3Td3k (2π)3k0f0i˜ f0im2 i 3Ek−1 3−c2 sEk, (23) where giis the degeneracy factor. After computing the spectra from Eq. (20), the two- and three-body decays of unstable particles are computed as in Ref. [80]. C. Neural networks To reduce the computational cost of the simulations, deep convolutional neural networks trained in Ref. [41] are utilized here for predicting final-state event-by-event observables at midrapidity. Each neural network takes the discretized initial energy density profile in the transverse-coordinate (x,y) plane as an input, and outputs one pT-integrated observable. Separate neural networks are used to predict flow coefficients vn, 034911-5 HENRY HIRVONEN et al. PHYSICAL REVIEW C 110, 034911 (2024) FIG. 1. Neural network validation test for the flow coefficients v2{2},v3{2},andv4{2}in 5.023 TeV Pb+Pb collisions. The networks were trained with the EbyE-EKRT data from Ref. [28] as described in Ref. [41]. The hydro results and the neural network validation results were obtained from 10k MC-EKRT initial states, which included hotspots and multiplicity fluctuations that were not present in the training data. The measured data are from the ALICE Collaboration [81]. charged particle multiplicities dNch/dη, and mean transverse momenta [pT]. Predicting flow observables with neural networks is many orders of magnitude faster than performing full hydrodynamic simulations. For example, predicting results for 107events takes only around 20 h with Nvidia Tesla V100 GPU. As the training data for the neural networks is from Ref. [28], and the predictions made using these networks emulate the same dynamics as the training data, the viscosities η/s, and ζ/s, and other parameters affecting the fluid-dynamical evolution are the same ones as in Ref. [28]. In Ref. [41], it was demonstrated that the neural networks work accurately when using the EbyE version of the EKRT model. However, it is nontrivial that the accuracy of the neural networks, which are trained by the EbyE-EKRT data from Ref. [28], and not from MC-EKRT, would extend to the MCEKRT initial states with hotspots, where the initial geometry can be significantly different. Therefore, the neural networks were validated by generating 10 000 MC-EKRT initial states and comparing the neural network predictions against (2+1)D fluid dynamical simulations for the 5.023 TeV Pb+Pb collision system. The validation tests for the flow coefficients v2,v3, and v4are shown in Fig. 1. The initial-state parameters used in the validation test were κsat =2.5, K=2.2, σ⊥=0.4 fm, and σh=0.2 fm. The obtained excellent agreement between the fluid dynamical simulations and neural network predictions illustrates the versatility of the neural networks with different initial conditions. Additionally, we have verified that the accuracy of the employed neural networks remains good for other training observables as well. IV. RESULTS In this section, we present the results of fluid-dynamical simulations with MC-EKRT initial states for midrapidity bulk observables, and compare the results against the earlier EbyE EKRT work [28]. All the fluid dynamical results are generated using our neural networks, and they contain 50 000 collision events, except the multiplicity distribution results, which are obtained from 150 000 events. As discussed in Sec. III C, the neural network results correspond to the fluid dynamical simulations with the matter properties and decoupling parameters from Ref. [28]. Therefore, any differences between the presented results are due to differences in the initial states. When examining the effects of the initial-state through final-state observables, it is important to remember that some observables are highly sensitive to the properties of the matter. For instance, the magnitude of flow coefficients is significantly influenced by the shear viscosity to entropy density ratio η/s. In contrast, the ratios of flow coefficients are less sensitive to such details, particularly the ratio between v3and v2, which can provide valuable insights into the geometry and structure of the initial state [82]. The effect of the Gaussian smearing width σ⊥is demonstrated in Fig. 2, where the ratios of the flow coefficients v2,v3, and v4in 5.023 TeV Pb+Pb collision system are shown as a function of centrality for different smearing widths. The MC-EKRT initial-state parameters are set to κsat =1.4, and K=2.5. Nucleon substructure is not included in these plots. As can be seen in the left panel, the magnitude of flow is sensitive to the Gaussian smearing width σ⊥. However, σ⊥ has only little impact on the ratios between the flow coefficients, as shown by the middle and right panels. Therefore, the parameter σ⊥is influencing the flow coefficients in a similar manner as the shear viscosity. Here, and in what follows, we adjust σ⊥to obtain the measured v2in midcentral collisions for all different MC-EKRT results. However, we want to emphasize that this is only done to illustrate the capabilities and uncertainties of MC-EKRT. To get the best overall fit to all different observables, a global analysis is needed, but this is beyond the purpose of this study. An intriguing aspect of the MC-EKRT model is the interplay between the saturation and conservation-law filters. The impact of different filters on the flow coefficients in 5.023 TeV Pb+Pb collisions is illustrated in Fig. 3. In all these scenarios, a value of K=2.5 is used, while the saturation 034911-6 EFFECTS OF SATURATION AND FLUCTUATING … PHYSICAL REVIEW C 110, 034911 (2024) FIG. 2. The effect of the Gaussian smearing width σ⊥on the two-particle flow coefficient v2{2}(left panel) and the ratios v3{2}/v2{2} (middle panel), and v4{2}/v2{2}(right panel) in 5.023 TeV Pb+Pb collisions. No nucleon substructure is included here. The experimental data for the ratios are computed based on the ALICE measurements for the two-particle flow coefficients [81]. parameter κsat is adjusted to achieve roughly identical charged particle multiplicities in central collisions. This corresponds to κsat =1.3 for the saturation-only case, and κsat =1.4for the other cases. The nucleon width is set according to the default parametrization from Eq. (4), i.e., σN=0.53 fm, and no nucleon substructure is introduced. For the saturation-first case σ⊥=0.3 fm, for the case with all filters at the same time σ⊥=0.4 fm, and for the saturation-only case σ⊥=0.3fm. The most notable feature in Fig. 3is the significant impact of saturation on the ratio between v3and v2. The case with only saturation reproduces the measured v2and v3most accurately, while the simultaneous application of all the filters leads to a clear underestimation of v3. When saturation is applied before other filters, the results approach those of the saturation-only scenario, as anticipated. The discrepancies in the v3/v2ratio arise from the geometrical differences in saturation and momentum conservation. Saturation does not allow geometrical overlap in the transverse plane. This leads to a more evenly distributed energy density profile. Energy conservation, on the other hand, gives no direct geometrical constraints. The stronger the saturation the more the eccentricity ε2is suppressed compared to the eccentricity ε3.The reduced eccentricity ε2can be compensated by decreasing the smearing width σ⊥so that the elliptic flow v2remains nearly unchanged, while ε3increases. This is reflected in the shown flow coefficients. It is also noteworthy that the v3/v2ratio is very similar between the MC-EKRT model with only saturation and the EbyE-EKRT model, which does not explicitly include momentum conservation. Since strong saturation appears to be necessary for matching the measured v3/v2ratio, we will now focus exclusively on the scenarios where saturation is applied first, followed by the conservation filters. This approach is also theoretically justified because, in principle, saturation should inherently account for conservation laws. However, achieving this would require implementing saturation through momentum-conserving multiparton distributions to all orders, which is not practically feasible. FIG. 3. The flow coefficients vn{2}as a function of centrality for 5.023 TeV Pb+Pb collisions. The simulation results with different MC-EKRT filter settings are compared against the ALICE measurements [81], and the EbyE-EKRT results from Ref. [28]. No nucleon substructure was included here. 034911-7 HENRY HIRVONEN et al. PHYSICAL REVIEW C 110, 034911 (2024) FIG. 4. The flow coefficients vn{2}as a function of centrality for 5.023 TeV Pb+Pb collisions. The MC-EKRT results with and without nucleon substructure are compared against the ALICE measurements [81], and the EbyE-EKRT results from Ref. [28]. Since saturation is sensitive to the nuclear overlap TATB (nuclear thickness function TAis the sum of TNs), the hotspots introduce interesting dynamics. With the hotspots, TAcan reach ≈10 times higher values than with the average nucleon geometry. Therefore, one would expect the saturation strength and the v3/v2ratio to increase when hotspots are included. The effect of hotspots on the flow coefficients is illustrated in Fig. 4, which compares two different hotspot parametrizations. The first parametrization uses the default nucleon width from parametrization Eq. (4), together with hotspots with width σh=0.2 fm. In this case, the MC-EKRT parameters are set to κsat =2.5, K=2.2, and σ⊥=0.4 fm. For the second parametrization, the nucleon width is obtained from Eq. (4), but this time a significantly stronger energy dependence with α=0.6 is used. This corresponds to a nucleon width σN= 0.75 fm for 5.023 TeV collision energy. This nucleon width is in line with the many global analyses, where values in the range ∼0.6–1.0 fm are preferred [33,37,44,83]. With a wider nucleon, a narrower hotspot with σh=0.15 fm is used together with parameters κsat =2.5, K=2.4, and σ⊥=0.25 fm. The saturation-first case from Fig. 3is here left as a reference curve. As expected, the addition of hotspots appears to increase the v3/v2ratio. The best overall fit to the measurements is obtained with the narrow hotspots, i.e., σh=0.15 fm, corresponding thus to the strongest saturation. In this case, the centrality dependence of v2, and v3matches nearly perfectly to the ALICE measurements [81], while maintaining a good agreement for v4. These findings suggest that the interplay between hotspots and saturation is crucial for the simultaneous description of the flow coefficients and especially of the v3/v2 ratio. In Fig. 5, the flow coefficients are shown for 2.76 TeV Pb+Pb collisions. The different curves correspond to the same cases as in Fig. 4,buttheKfactor is adjusted to obtain a reasonable agreement with the measured charged particle multiplicity. The obtained values are K=2.5fortheσh=0.2 fm case, while the σh=0.15 fm and the no-hotspots cases both use K=2.7. The agreement between the data and the results is quite similar to the 5.023 TeV collision energy results. At both energies, the narrow-hotspot case with σh=0.15 fm can describe the measured flow coefficients well, while the centrality dependence of v2is slightly off for the σh=0.2fm case. From Figs. 4and 5it can be seen that MC-EKRT with the nucleon substructure captures the energy dependence of the flow coefficients significantly better than the EbyE-EKRT model. In Fig. 6, the charged particle multiplicity as a function of centrality is shown for the same initial-state parametrizations in 2.76 TeV and 5.023 TeV Pb+Pb collisions. The agreement between the results and the ALICE measurements [84,85] is good in all cases, even though there are some minor discrepancies in the centrality behavior. The initial state without hotspots seems to produce slightly too weak a centrality dependence, while, with the hotspots, the centrality dependence is a bit too steep. However, these are small differences, and further improvements could be obtained by fine tuning the matter properties and initial state parameters. The MC-EKRT approach adds minijet-multiplicity- originating saturation-scale fluctuations to the EKRT initial state. These fluctuations, together with hotspot fluctuations, should in principle increase the hadron multiplicity fluctuations in the most central collisions. This effect is studied in Fig. 7, where the charged hadron multiplicity distributions from MC-EKRT with and without hotspots are compared against the EbyE-EKRT results, which do not contain multiplicity-originating fluctuations of the saturation scale or hotspots. To make the results comparable with the V0 amplitude measured by ALICE [86], they are normalized to have approximately the same mean as the V0 amplitude. AsshownalsoinRef.[17], the EbyE-EKRT results almost completely miss the high-multiplicity tail in the distribution. The addition of the further saturation scale fluctuations indeed enhances the high-multiplicity tail in the distribution, and therefore improves the agreement with the measurements as one would expect. The addition of the hotspots is important also for this observable, as it increases the fluctuations and high-multiplicity tail further, leading to a very good agreement with the ALICE data. 034911-8