Dynamic modeling for heavy-ion collisions

Quantifying the phase structure of hot and dense QCD matter is one of the primary goals in relativistic nuclear physics. Lattice QCD calculations Bhattacharya:2014ara ; Bazavov:2011nk provide conclusive information that hot nuclear matter at zero net baryon density transits from the hadron gas to the Quark-Gluon Plasma (QGP) phase via a smooth crossover at the pseudo-critical temperature $T_{pc}=156.5\pm 1.5$ MeV HotQCD:2018pds . Moreover, the same type of phase structure remains for small-to-moderate values of the net baryon chemical potential $\mu_{B}$ HotQCD:2018pds . It was conjectured that the quark-hadron transition turns from a crossover to a first-order phase transition at some finite baryon chemical potential, suggesting the existence of a critical point Asakawa:1989bq ; Bzdak:2019pkr . However, the first principles lattice QCD calculations can only provide limited guidance on the existence and location of the QCD critical point in the $T-\mu_{B}$ phase diagram because of the fermion sign problem at finite net baryon density deForcrand:2009zkb . Meanwhile, exploration of the dense quark matter properties at large net baryon density is of particular importance since the gravitational waves detection, emerging from black hole/neutron star mergers, now give stringent constraints on the properties of the compact stars, including the nuclear matter equation of state (EoS) LIGOScientific:2017vwq ; LIGOScientific:2018cki ; Baym:2017whm ; Dexheimer:2020zzs .

To establish quantitative connections between the QCD phase structure and measurements from relativistic heavy-ion collisions over an extensive collision energy range, we need to model the entire dynamical evolution of the heavy-ion collisions. Relativistic hydrodynamics, incorporated with nuclear matter EoS, viscosity, and initial state fluctuations, has been a precision tool to understand the macroscopic dynamics of the strongly coupled QGP. The fluid dynamic description is connected with hadronic transport approaches to microscopically describe the system’s evolution in the dilute phase. Such a hydrodynamics + hadronic transport hybrid theoretical framework has successfully described and even predicted various flow correlation measurements with remarkable precision Shen:2014vra ; Bernhard:2016tnd ; Schenke:2020mbo ; JETSCAPE:2020mzn . Phenomenological studies with precise flow measurements of the hadronic final state from the Beam Energy Scan (BES) program can elucidate the collective aspects of the baryon-rich QGP and constrain the QGP EoS and transport properties, such as its viscosity and charge diffusion coefficients Denicol:2018wdp ; McLaughlin:2021dph . Strange quarks produced in relativistic nuclear collisions are an interesting probe for studying the evolution of the collisions. Since the collision system is strangeness neutral, the $s\bar{s}$ pair production mechanism is sensitive to the properties of strongly-interacting matter, especially at the early time Pratt:2018ebf . The partonic and hadronic production channels are very different and may signal the onset of deconfinement and the quark-gluon plasma Blume:2017icv . Furthermore, unstable resonance states of strange baryons can provide detailed information about hadronic reactions in the late-stage of heavy-ion collisions Oliinychenko:2021enj . This proceeding will highlight the recent developments of the hybrid dynamical frameworks, emphasizing the strangeness-related observables.

Relativistic heavy-ion collisions with collision energies $\sqrt{s_{\mathrm{NN}}}\gtrsim 20$ GeV can probe the nuclear matter properties with a net baryon chemical potential up to $\sim 300$ MeV, where the QGP and hadronic phases are connected with a smooth crossover Andronic:2017pug ; Adamczyk:2017iwn . At $\sqrt{s_{\mathrm{NN}}}\sim\mathcal{O}(10)$ GeV, heavy-ion collisions violate longitudinal boost invariance, and the overlapping time for the two nuclei to pass through each other is significant compared to the total collision lifetime Karpenko:2015xea ; Shen:2017bsr ; Shen:2020jwv . Therefore, it is crucial to develop initial conditions with non-trivial longitudinal dynamics. The complex 3D collision dynamics can be approximated by parametric energy depositions based on the collision geometry Hirano:2005xf ; Bozek:2017qir ; Shen:2020jwv ; Sakai:2020pjw . Non-trivial dynamics could be included by modeling the energy loss in individual nucleon-nucleon collisions based on either string deceleration Shen:2017bsr ; Bialas:2016epd or full transport simulations Pang:2012he ; Karpenko:2015xea ; Du:2018mpf . Further theory developments exist to understand the early-stage baryon stopping from the Color Glass Condensate-based approaches in the fragmentation region Li:2018ini ; McLerran:2018avb and a holographic approach at intermediate couplings Attems:2018gou . In the meantime, dynamical initialization schemes have been developed to interweave initial-state and hydrodynamics on a local basis to model the extended interaction region between the two colliding nuclei at low collision energies Shen:2017ruz ; Shen:2017bsr ; Du:2018mpf ; Akamatsu:2018olk . These schemes have also been applied to study small systems and jet-medium interactions Okai:2017ofp ; Kanakubo:2019ogh .

Figure 1 shows a series of snapshots of one central Au+Au collision at $\sqrt{s_{\mathrm{NN}}}=39$ GeV flowing through the QCD phase diagram. The numerical simulation is based on 3D MC-Glauber with string deceleration model Shen:2017bsr . The first row presents the early dynamical initialization stage $\tau\leq 1.5$ fm/$c$, during which the two colliding nuclei pass through each other. The collisions among nucleons continuously deposit energy-momentum currents to form the strongly coupled QGP. Because most of the baryon charge sources are at the end of strings, their deposition times are later than those of the energy-momentum currents. Hence, the hydrodynamic medium first heats up at $\mu_{B}\sim 0$. Around $\tau\sim 1.2$ fm/$c$, the baryon charges start to be doped into the system and drive the system to flow to a large $\mu_{B}$ region in the phase diagram. For central collisions at $\sqrt{s_{\mathrm{NN}}}=39$ GeV, the duration of energy-momentum and net baryon charge depositions last for about 1.5-2 fm/$c$. During this phase, the fireball reaches $T_{\mathrm{max}}\sim 250$ MeV and $\mu_{B,\mathrm{max}}\sim 200$ MeV at the mid-rapidity. For fireball in the forward space-time rapidity region $1<\eta_{s}<1.8$, the peak temperature of the medium drops to about 200 MeV while the medium reaches out to a larger $\mu_{B}$ region, $\mu_{B}\sim 300$ MeV. For $\tau\geq 2$ fm/$c$, the fireball expands hydrodynamically in 3D and lives up to $10.5$ fm/$c$ before all the fluid cells convert to hadrons in the dilute hadronic phase.

Because most hybrid dynamical frameworks assume grand-canonical ensemble (GCE) at the particlization, the hadron chemistry is mainly controlled by the distributions of temperature and chemical potentials on the constant energy density particlization hyper-surface. Furthermore, the late-stage hadronic transport model can further modify the relative particle abundances via inelastic scatterings and baryon-anti-baryon annihilation. Therefore, the yield ratios of identified particles provide strong constraints on the particlization condition in hybrid simulations Monnai:2019hkn ; Alba:2020jir . Assuming GCE, one can map the net electric charge $(Q)$, strangeness $(S)$, and baryon $(B)$ chemical potentials to identified particle ratios as follows, $\pi^{-}/\pi^{+}\propto\exp(-2\mu_{Q})$, $K^{-}/K^{+}\propto\exp[-2(\mu_{Q}+\mu_{S})]$, and $\bar{p}/p\propto\exp[-2(\mu_{Q}+\mu_{B})]$. The ratios $K^{-}/\pi^{-}\propto\exp(-\mu_{S})$ and $K^{+}/\pi^{+}\propto\exp(+\mu_{S})$ provide direct information about the strangeness chemical potential. The yield ratios of multi-strange baryons to their antiparticles $\bar{\Lambda}/\Lambda\propto\exp[-2(\mu_{B}-\mu_{S})]$, $\bar{\Xi}^{+}/\Xi^{-}\propto\exp[-2(\mu_{B}-2\mu_{S}-\mu_{Q})]$, and $\bar{\Omega}/\Omega\propto\exp[-2(\mu_{B}-3\mu_{S})]$ contain different weights on $\mu_{S}$.

Figure 2 shows the comparison between the preliminary calculations from the 3D string deceleration model Shen:2017bsr and the identified particle yield ratios in central Au+Au collisions at the RHIC BES energies. The model calculations with a constant switching energy density $e_{\mathrm{sw}}=0.26$ GeV/fm³ reproduce the collision energy dependence of the data from 7.7 GeV to 200 GeV. The key theoretical ingredients that lead to this good description are initial-state baryon stopping plus the constraints $n_{S}=0$ and $n_{Q}=0.4n_{B}$ for Au+Au collision Monnai:2019hkn ; Monnai:2021kgu . The strangeness neutrality condition $n_{S}=0$ results in a linear correlation between the strangeness chemical potential and $\mu_{B}$, $\mu_{S}\sim\mu_{B}/3$ Monnai:2021kgu . The collision energy dependence of $\mu_{S}$ can be constrained with the $K^{+}/\pi^{+}$ and $K^{-}/\pi^{-}$ ratios. Below $\sqrt{s_{\mathrm{NN}}}=5$ GeV, the $K^{+}/\pi^{+}$ ratio becomes smaller because of the canonical suppression of the strange quark production Andronic:2017pug ; STAR:2021hyx . The small $\mu_{Q}$ from the constraint $n_{Q}=0.4n_{B}$ quantitatively reproduces the $\pi^{-}/\pi^{+}$ ratios. The model calculation underestimates the ratios among multi-strangeness baryons for $\sqrt{s_{\mathrm{NN}}}<40$ GeV, which hints that these strangeness baryons could have a higher chemical freeze-out energy density than the non-strange hadrons Bellwied:2018tkc illustrated in Fig. 2d.

Over the past decade, the RHIC Beam Energy Scan program has excited a wave of theory developments towards a (3+1)D paradigm of dynamical modeling of relativistic nuclear collisions. The advancements in 3D initial-state models, equation of state at finite net baryon density, dynamical initialization schemes enable us to quantify initial baryon stopping and study the collectivity of the QGP in a baryon-rich environment. To confront precision flow measurements from the upcoming RHIC BES program phase II, our community has adopted the Bayesian Inference method as a standard approach to systematically constrain the QGP’s thermodynamic and transport properties and initial-state fluctuation spectrum Bernhard:2019bmu ; JETSCAPE:2020mzn ; Nijs:2020roc . Meanwhile, the increasing theoretical and numerical complexities in 3D hybrid models could pose significant challenges in future high-dimensional Bayesian analysis.

The current RHIC fixed-target experiment and the future FAIR/NICA experimental programs push the exploration of the QCD phase diagram beyond $\mu_{B}=400$ MeV, where guidance from lattice QCD on EoS is very limited. Searching for experimental signals of critical point and first-order phase transition requires the theoretical frameworks to be calibrated by particle production and flow measurements to provide reliable baseline expectations without critical phenomena Zhao:2020irc ; Oliinychenko:2020znl ; Pratt:2020ekp ; Vovchenko:2021kxx . Further studying the effects of out-of-equilibrium evolution of critical fluctuations Rajagopal:2019xwg ; Du:2020bxp ; Dore:2020jye ; Nahrgang:2020yxm ; Du:2021zqz on net proton high-order cumulants and light nuclei productions Sun:2017xrx as functions of collision energy and rapidity intervals would shed light on the QCD phase structure at large baryon density.

This work is support in part by the National Science Foundation (NSF) under grant number PHY-2012922 and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-SC001346 and within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration. The calculations in this work used resources provided by the Open Science Grid Pordes:2007zzb ; Sfiligoi:2009cct , which is supported by the NSF award 2030508.