Charged Particle Spectra and Yields

Relativistic heavy ion collision experiments are designed for the search and study of the QGP. In head-on relativistic heavy ion collisions, two nuclei can be represented as two thin disks approaching each other at high speed because of the Lorentz contraction effect in the moving direction. During the initial stage of collisions, the energy density is higher than the critical energy density obtained from lattice QCD calculations, whereby quarks and gluons are deconfined from nucleons and form QGP. The large cross section of interaction may lead to the thermalization of the QGP. At this stage, high transverse momentum jets and heavy-flavor pair are produced because of the large momentum transfer. Subsequently, the QGP expands and cools down, entering into the mixed-phase expansion. The chemical freeze-out point is formed after the inelastic interaction stops, which means that thereafter particle yields and ratios do not change. After chemical freeze-out, the elastic interaction between hadrons and resonance decays [31] change the p_T distribution of particles. The particles finally freeze out of the system after the elastic interaction stops, at the so-called kinetic freeze-out point. Studying the bulk properties of the system, such as the spectra, the yields (dN/dy), particle ratios, and freeze-out properties, can provide insight into the particle production mechanisms and evolution of QCD matter.

In experimental observations, first invariant yields of various particles are presented as a function of transverse momentum p_T. Figure 6 shows the invariant yields of pions (π^±), kaons (K^±), protons (p), anti-protons ($\overline{\text{p}}$), phi-mesons (ϕ), lambda baryons (Λ), anti-lambda baryons ($\overline{\text{Λ}}$), cascades (Ξ^-), anti-cascades (${\overline{\text{Ξ}}}^{+}$), Omegas (Ω^-), anti-Omegas (${\overline{\text{Ω}}}^{+}$), deuterium (d), and anti-deuterium ($\overline{\text{d}}$). The results are shown for Au+Au collisions at $\sqrt{s_{\text{NN}}}=19.6\text{ GeV}$ in four collision centralities: 0–10%, 20–40%, 40–60%, and 60–80% [32-35]. The invariant yields show a decrease with increasing p_T, moving from central to peripheral collisions. The curves represent the blast wave fits to the spectra [36].

Fig. 6Full-screen View

(Color online) Invariant yields of identified particles measured as a function of transverse momentum in various collision centralities in Au+Au collisions at $\sqrt{s_{\text{NN}}}=19.6\text{ GeV}$ [32-35]. Filled points represent particles and open points are antiparticles. The curves are fits to the data points based on the thermal blast wave distribution [36]

The yields, dN/dy, are obtained by integrating these measured spectra and interpolating through fitting functions where the measurements are not available. Figure 7 shows the energy dependence of particle yields for π^±, K^±, p, $\overline{\text{p}}$, ϕ, Λ, $\overline{\text{Λ}}$, Ξ^-, ${\overline{\text{Ξ}}}^{+}$, Ω^-, ${\overline{\text{Ω}}}^{+}$, d, and $\overline{\text{d}}$. Results from STAR BES-I [32-35] are compared with previously published STAR results at higher energies and those of other world experiments. The reader is referred to the topical review in [6] on data collection. The yields of anti-baryons increase rapidly with increasing collision energy, demonstrating the increasing contribution of pair production. However, the yields of baryons and K⁺ show a non-trivial energy dependence in the BES energy range, indicating the interplay of baryon stopping/association and pair production.

Fig. 7Full-screen View

(Color online) Energy dependence of particle yields dN/dy in central collisions [6]. Solid points represent particles, and open points are the antiparticles. The data are normalized by the number of participants

The hadron yields retain the footprint of the hot and dense hadronic matter during the evolution of the collision system, presumably because the system undergoes a crossover at phase transition [32]. As observed experimentally, the relative abundances of hadrons follow the thermal distribution at common temperatures and baryonic chemical potentials, such that the thermal fit can fix the temperature and baryonic chemical potential [32]. The temperatures of chemical freeze-out (T_ch for central Au+Au collisions at different collision energies are shown in Fig. 8. With increasing energy, the T_ch increases and becomes constant at ～ 160 MeV after $\sqrt{s_{\text{NN}}}=11.5\text{ GeV}$. The parameters extracted from net-proton higher moments [38] are also presented in the figure. They are consistent with the results from the hadron yield fit. The extracted parameters from BES data are also consistent with the results from the lattice QCD calculations and thermal fit to the global hadron yield data [37]. In covering the RHIC BES program, STAR fixed target program, and future experimental facilities shown in the figure, a more precise description of the QCD phase diagram is presented.

Fig. 8Full-screen View

(Color online) A summary of the chemical freeze-out temperature T_ch(μ_B) distribution [37]. Data points are from the 0–5% central Au+Au collisions at STAR BES [32, 38]

One of the foundations of the BES program is the promise of a sweeping variation of the chemical potential across the QCD phase diagram through changes in the beam energy of heavy ion collisions, whereby the chemical potential is extracted empirically from the final-state particle distributions. An important subject in itself is how the baryons are shifted from target and projectile rapidity to midrapidity. A puzzling feature of ultra-relativistic nucleus-nucleus collisions is the experimental observation of substantial baryon asymmetry in the central rapidity (midrapidity) region both at the energies of RHIC [39, 40, 32] and large hadron collider (LHC) ($\sqrt{s_{\text{NN}}}=900\text{ GeV}$) [41, 42]. Such a phenomenon is impressive, as the baryon number is strictly conserved; therefore, net-baryon number cannot be created in the system and must come from the colliding targets and projectile. In the conventional picture, the valence quarks carry the baryon quantum number in the nucleus. At sufficiently high energies these valence quarks are expected to pass through, ending up far from midrapidity in the fragmentation regions [43]. The RHIC BES program covers a wide range of baryon stopping, over an order of magnitude of net proton yields at midrapidity [40, 32].

Figure 9 presents the net-proton yields at midrapidity in Au+Au collisions at $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$ to 200 GeV [44]. For all centralities in heavy ion collisions, the midrapidity net-baryon density follows an exponential distribution with the variable $\delta y=y_{\text{b}}-y_{\text{cm}}$, where y_b is the beam rapidity, and y_cm is the center-of-mass rapidity. This variable $\delta y$ is referred to as the “rapidity loss” which for the midrapidity protons produced in a collider experiment is equal to beam rapidity: $\delta y=y_{\text{b}}$ as $y_{\text{cm}}=0$. A single collision energy therefore gives rise to a single point in Fig. 9. The data points at each centrality can be fitted with an exponential function A×exp(-α_Bδy). The baryon stopping is often characterized by average rapidity loss [45], which is characterized by a complicated beam energy dependence and is usually skewed by the large proton yields close to beam rapidity, leading to the conclusion [45] that the “rapidity loss” of projectile baryons at RHIC breaks the linear scaling observed at lower energies. Another way of characterizing the baryon stopping is to use the $\overline{\text{p}}/\text{p}$ ratio [46, 41, 42]. Both pair production and baryon stopping contribute to this ratio. Most of the dynamic models of heavy ion collisions parametrize the baryon stopping to reproduce the experimental data although at a fundamental level, a lack of understanding on how baryons are stopped prevails. A recent modeling of heavy ion collisions indicates that the inclusion of the baryon junction is essential for describing net-proton density at RHIC [47]. Clearly some of the earlier implementations of baryon junctions [48, 49], which attempted to match the earlier experimental data with certain parameter tunes, do not reproduce the experimental results presented in Fig. 9.

Fig. 9Full-screen View

(Color online) Exponential dependence of midrapidity (y≈0) baryon density per participant pair in heavy ion collisions with y_b equal to the rapidity difference between beam and detector midrapidity (δy) [44]. An exponential fit function of A×exp(-α_Bδy) is also included. The figure is taken from [44]

Strangeness Production

Strange hadrons serve as excellent probes for the physics of QCD phase boundary and search for the onset of deconfinement. Strangeness enhancement in heavy ion hadron collisions has long been suggested as a signature of the quark-gluon plasma [50-52], motivating its measurement in many experiments at different accelerator facilities. In general, the yields of strange hadrons in nuclear collisions are close to those expected from statistical models [53-55]. The precise measurement of these yields in phase-I of RHIC BES experiments has led to a better understanding of the strange quark production mechanisms in nuclear collisions and a more reliable extraction of the chemical freeze-out parameters [32], as shown in Fig. 8. In the higher beam energies, formation of a thermalized system is expected and strangeness is abundantly produced. However, at lower beam energies, strangeness production is less, requiring special attention and local treatment of the canonical ensemble. This part is further discussed in relation to Fig. 11, with the ϕ(1020) meson with zero net strangeness number (S=0) offering a unique opportunity to scrutinize the thermodynamic properties of strange quarks in the hot and dense QCD environment [56].

Fig. 11Full-screen View

(Color online) (a) ϕ/K^- and (b) $\varphi /{\text{Ξ}}^{-}$ ratios in central heavy ion collisions as a function of collision energy, $\sqrt{s_{\text{NN}}}$, compared with various thermal and transport model calculations. Figure taken from [78]

The precise measurement of strange hadron production at different p_T ranges and centralities in heavy ion collisions are also crucial for a better understanding of the production mechanism and medium properties created in the system. At high p_T, the nuclear modification factor R_CP of various particles at top RHIC energy is observed to be much less than unity [57-59], indicating a significant energy loss of the scattered partons in the dense nuclear matter, known as “jet quenching” [60]. At intermediate p_T, the baryon to meson ratios, p/π and $\text{Λ}/{\text{K}}_s^0$ are found to be larger than unity and much higher than those observed in the peripheral A+A and elementary collisions. This baryon to meson ratio enhancement can be explained by the recombination/coalescence models which require constituent quarks in the partonic medium to coalesce into hadrons, or soft and hard partons to recombine into hadrons [61-63]. Thus, measurements of R_CP and baryon to meson ratios of strange hadrons are one of the corner-stone pieces of evidence for the formation of the strongly interacting QGP medium. The precise measurement of these variables in heavy ion collisions at lower beam energies can potentially reveal the medium properties at finite μ_B and help locate the collision energy for the onset of deconfinement.

Apart from light hadrons, Fig. 7 also shows the energy dependence of strange particle yields at midrapidity for K^±, ϕ, $\text{Λ}\left(\overline{\text{Λ}}\right)$, ${\text{Ξ}}^{-}\left({\overline{\text{Ξ}}}^{+}\right)$ and ${\text{Ω}}^{-}\left({\overline{\text{Ω}}}^{+}\right)$ from central heavy ion collisions. Results from STAR BES-I [34, 33] are compared with previously published STAR results at higher energies and other corresponding world data including experiments at the Alternating Gradient Synchrotron (AGS) and Conseil Européen pour la Recherche Nucléaire (CERN) [58, 64-73]. The yields $\text{d}N/\text{d}y/\langle N_{\text{part}}/2\rangle$ of the anti-hyperons ($\overline{\text{Λ}}\mathrm{,}{\overline{\text{Ξ}}}^{+}\mathrm{,}{\overline{\text{Ω}}}^{+}$) and ϕ meson increase rapidly with increasing energy, with a non-trivial energy dependence on Λ, ${\text{Ξ}}^{-}$ and Ω^- yields. The ${\text{Ξ}}^{-}$ and Ω^- yield first increases with energy from 7.7 GeV to 19.6 GeV, then remains almost constant up to energies around 39 GeV, then rising again toward higher energies. The Λ yield first decreases from 7.7 GeV to 39 GeV, then increases toward higher energies. The Λ behavior is similar to the trend displayed by proton in these measured energy regions [32], reflecting a significant increase in baryon density at lower collision energies. The observed Λ behavior can be the result of the interplay of the pair production of Λ-$\overline{\text{Λ}}$ and associated Λ production along with a kaon, with the former increasing strongly with an increase in collision energy and the latter increasing strongly with an increase in net-baryon density.

Figure 10 shows the energy dependence of $\text{Λ}\left(\overline{\text{Λ}}\right)$ and ${\text{Ξ}}^{-}\left({\overline{\text{Ξ}}}^{+}\right)$ midrapidity yield ratio to that of pions in central Au+Au collisions from RHIC STAR BES, and the existing data from various experiments [33, 64-67, 73-75] as well as the calculations from hadronic transport models (UrQMD 1.3, hadron-string dynamics (HSD)) and statistical hadron gas model (SHM) [76, 77, 55]. The STAR BES data are in good agreement with the trend displayed by the existing experimental data. The hadronic models (UrQMD 1.3 and HSD) seem to reproduce the Λ/π data, indicating that hadronic rescatterings might play an important role in hyperon production in heavy ion collisions at this energy range. However the default UrQMD (v1.3) fails to reproduce the $\text{Ξ}$/π ratio due to a smaller $\text{Ξ}$ yield in the model [77]. On the other hand, the SHM model predictions agree well with data across the entire energy range from AGS to top RHIC. The SHM model used here is based on a grand-canonical ensemble and assumes chemical equilibrium. The energy dependence of the parameters T_ch and μ_B in the model are obtained by a smooth parametrization of the original fitting parameters to the midrapidity particle ratios from heavy ion experiments at super proton synchrotron (SPS) and RHIC. Both Λ/π and ${\text{Ξ}}^{-}/\pi$ ratios show a maximum at ～8 GeV, which seems to be consistent with the maximum net-baryon density at freeze-out at this collision energy.

Fig. 10Full-screen View

(Color online) Energy dependence of Λ, $\overline{\text{Λ}}$, Ξ^-, and ${\overline{\text{Ξ}}}^{+}$ midrapidity yield ratio to that of pions (1.5 (${\pi }^{+}+{\pi }^{-}$)) in central Au+Au collisions from STAR BES, compared with the existing data from various other experiments [33]

Thermodynamic properties of strange quarks play an important role in understanding the QCD matter equation of state (EOS) in high-density regions. In statistical thermal models, grand canonical ensemble (GCE) and canonical ensemble (CE) statistical descriptions are applied differently to conserve the strangeness number in computing the final state particle yields. It has been argued that at lower energies, strangeness number should be conserved locally on an event-by-event basis as described by CE, which leads to a reduction in the yields of hadrons with non-zero strangeness number (“Canonical Suppression”) but not for the ϕ(1020) meson with zero net strangeness number (S=0) [79]. Fig. 11 shows the measurements of ϕ/K^- and $\varphi /{\text{Ξ}}^{-}$ ratio in the central heavy ion collisions as a function of collision energy [78, 33, 80, 81] compared with various thermal and transport model calculations [76, 77, 82]. As shown in the plot, both GCE and CE models describe the measured ratios $\sqrt{s_{\text{NN}}}$ at greater than 7.7 GeV, whereas clearly GCE fails when the collision energies approach the production threshold (2.89 GeV for ϕ and 3.25 GeV for Ξ^-). The measurements favor CE calculations with a small strangeness correlation length (r_c), necessitating more detailed investigation for precise and differential data.

In addition to the thermal model, transport model calculations from modified UrQMD with high mass strange resonances can reasonably reproduce the data in Fig. 11, implying that the feed-down effect is relevant [82, 83]. In heavy ion collisions, the near/sub-threshold production of multi-strange hadrons can be achieved from the multiple collisions of nucleons, produced particles, and short-lived resonances. However, particle production below the free nucleon-nucleon (NN) threshold is expected to be sensitive to the stiffness of the nuclear EoS at high density [84].

Figure 12 panel (a) shows the nuclear modification factor, R_CP, of ${\text{K}}_S^0$, in Au+Au collisions at STAR BES from 7.7 GeV to 39 GeV [33]. For p_T ≈ 4 GeV/c, the ${\text{K}}_S^0$ R_CP is below unity at $\sqrt{s_{\text{NN}}}=39\text{ GeV}$, which is similar to the observation at top RHIC energy although the lowest R_CP value is larger here. The ${\text{K}}_S^0$ R_CP at p_T > 2 GeV/c keeps increasing with decreasing collision energies, indicating that the partonic energy loss effect becomes less important. Eventually, the value of ${\text{K}}_S^0$ R_CP tends to increase at large p_T at $\sqrt{s_{\text{NN}}}=11.5\text{ GeV}$ and 7.7 GeV although the maximum accessible p_T is smaller at these two energies. This suggests that the cold nuclear matter effect (Cronin effect) starts to take over at these energies and enhances all the hadron yields at intermediate p_T (to ≈ 3.5 GeV/c). Similar to the observation for identified charged hadrons, the energy evolution of strange hadron R_CP reflects the decreasing partonic effect with decreasing beam energies [60].

Fig. 12Full-screen View

(Color online) (a) ${\text{K}}_S^0$ nuclear modification factor, R_CP, at midrapidity in Au+Au collisions at STAR BES from 7.7 GeV to 39 GeV [33]. (b) $\overline{\text{Λ}}/{\text{K}}_S^0$ ratio as a function of p_T in central Au+Au collisions at STAR BES [33]. (c) baryon-to-meson ratio, Ω/ϕ, as a function of p_T in central Au+Au collisions from STAR BES [33]

Figure 12 panel (b) shows the $\overline{\text{Λ}}/{\text{K}}_S^0$ ratios as a function of p_T in central Au+Au collisions at STAR BES from $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$ to 39 GeV [33]. The $\overline{\text{Λ}}$ is chosen because it is a newly produced baryon in the baryon-rich medium created by the lower BES energies. An enhancement of baryon-to-meson ratios is observed at intermediate p_T in central A+A collisions compared to peripheral A+A or p+p collisions at the same energy for energies $\sqrt{s_{\text{NN}}}\ge 19.6\text{ GeV}$. The maximum value of $\overline{\text{Λ}}/{\text{K}}_S^0$ reaches the maximum value of unity at p_T ≈ 2.5 GeV/c for most central collisions, whereas in peripheral collisions, the maximum value is significantly lower, only about 0.3–0.5 which is not shown in the plot. The enhancement of baryon-to-meson ratio in central collisions in these energies is interpreted as being a consequence of hadron formation through parton recombination and parton collectivity. Therefore, the baryon-to-meson ratios are expected to be sensitive to the parton dynamics of the collision system. Unfortunately, for $\sqrt{s_{\text{NN}}}\le 11.5\text{ GeV}$, the statistics on different centralities and maximum p_T are limited. Hence, whether the baryon-to-meson enhancement still persists at these low energies remains unclear with the current data.

Figure 12 panel (c) shows the baryon-to-meson ratio, Ω/ϕ, as a function of p_T in central Au+Au collisions from $\sqrt{s_{\text{NN}}}=11.5\text{ GeV}$ to 39 GeV and $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$ 0–60% centrality [33]. For energies $\sqrt{s_{\text{NN}}}\ge 19.6\text{ GeV}$, the measured data follow each other closely and also the previous measurement from 200 GeV, which is consistent with coalescence and recombination dynamics over a broad p_T range of 1–4 GeV/c [34]. The ratios at $\sqrt{s_{\text{NN}}}\le 11.5\text{ GeV}$ seem to deviate from the trend observed at higher beam energies. In particular, the ratios at 11.5 GeV appear to take a downward turn around p_T of 2 GeV/c, whereas those at higher beam energies such as 39 GeV, peak at p_T of 3 GeV/c or above. Since the Ω and ϕ particles have small hadronic rescattering cross sections, the change in these Ω/ϕ ratios may indicate a significant change in the hadron formation dynamics and/or the strange quark p_T distribution at the lower energies.

Collectivity

Collective observables, including radial and anisotropic flow, are powerful tools for extracting parameters of the EOS and understanding the properties of the medium created by high-energy nuclear collisions [4, 85, 86]. In this section, the energy dependence of $v_1\mathrm{,}{v}_2$, their scaling, and EOS parameters are discussed.

The elliptic flow scaled by the number of constituent quarks (NCQ), $v_2/n_q$, for the copiously produced hadrons π^± (squares), K^± (crosses), p and $\overline{\text{p}}$ (circles), is shown as a function of the scaled transverse kinetic energy (m_T-m₀)/m_q in Fig.13. The data are from 10–40% mid-central Au+Au collisions at RHIC. Data points from collisions at 27 GeV and 54.4 GeV are shown as open and closed symbols, respectively. The colored dashed lines, also displayed in the figure, represent the scaling fit to the data from pions, kaons, and protons in Au+Au collisions at 7.7 GeV, 14.5 GeV, 27 GeV, 54.4 GeV, and 200 GeV for both positively and negatively charged particles [90, 91]. Although the overall quark number scaling is evident, the best scaling is reached by the RHIC top energy $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ collisions [92]. As collision energy decreases, the scaling deteriorates. Particles and antiparticles are no longer consistent with the single-particle NCQ scaling due to the mixture of transported and produced quarks [91]. More detailed discussions on the effects of transported quarks on collectivity can be found in [93, 94]. An important evidence for QGP formation in high energy collisions at RHIC, is the observed NCQ scaling that originates from partonic collectivity [95, 92, 96]. Interestingly, in the analysis of the elliptic flow of light nuclei in low- and intermediate-energy nuclear reactions, a similar scaling law exists, i.e., the elliptic flow of light nuclei is scaled according to the number of constituent nucleons [97]. Inspired by [97], the STAR experiment of relativistic heavy ion collisions [98] also confirms the nucleon number scaling of the elliptic flow of light nuclei, that is, it fulfills the theoretical prediction of [97]. The similarity between NCQ scaling of elliptic flows and nucleon-number scaling of light nuclei lies in the merger mechanism of hadron formation or nucleosynthesis, whereas the difference lies in whether the merger is at the quark or nucleon level. Notably, the LHC-ALICE collaboration reported measurements of higher-order anisotropic flows [99], providing for the first time experimental measurements of triangular flow v₃. The NCQ scaling of higher order collective flows is theoretically discussed in [100] and confirmed in experimental measurements [101], which can also be regarded as a further probe of QGP. Experimental efforts on the measurement of v₃ in BES energies have also been conducted [102]. Based on the data [102], detailed tests on NCQ scaling seem promising.

Fig. 13Full-screen View

(Color online) Elliptic flow (v₂) scaled by number of constituent quarks (nq); ($v_2/n_q$) is shown as a function of scaled transverse kinetic energy $\left(m_T-m_0\right)/n_q$ for pions, kaons, and protons from Au+Au collisions in 10–40% centrality at $\sqrt{s_{\text{NN}}}=3\text{ GeV}$, 27 GeV, and 54.4 GeV for positively charged particles (left panel) and negatively charged particles (right panel). Colored dashed lines represent the scaling fit to data from Au+Au collisions at 7.7 GeV, 14.5 GeV, 27 GeV, 54.4 GeV, and 200 GeV from the STAR experiment at RHIC [87-89]. Statistical and systematic uncertainties are shown as bars and gray bands, respectively

At low energies of $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ Au+Au collisions, a totally different scaling behavior emerges, as shown in Fig. 13. Unlike that observed in high-energy collisions, all v₂ values are negative, which is a characteristic of nuclear shadowing in such non-central collisions. There is no sign of NCQ scaling at this low energy [103]. These results clearly indicate different properties for the matter produced. With the baryonic mean field, hadronic transport model calculations from jet AA microscopic transport model (JAM) [104] and UrQMD [76, 77] reproduce the observed negative values of v₂ for protons as well as Λs. In other words, in Au+Au collisions at $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$, partonic interactions no longer dominate, and baryonic scatterings take over, indicating that predominantly hadronic matter is created in such low-energy collisions.

Now, focusing on the p_T-integrated results, we jointly examine v₁ and v₂. The collision energy dependence of directed and elliptic flow is summarized in Fig. 14, where panel (a) shows the slope of the p_T-integrated directed flow at midrapidity, $\text{d}{v}_1/\text{d}y|y=0$, for π, K, p, Λ, and multi-strange hadrons ϕ and Ξ^- from Au+Au collisions in the 10–40% centrality interval. Here, K and π represent the combined results of K^± and ${\text{K}}_S^0$, and π^±, respectively. Panel (b) displays the p_T-integrated v₂ at midrapidity for π, K, p, and Λ as open squares, filled triangles, filled circles, and open circles, respectively. Due to partonic collectivity in Au+Au collisions at high energy [106], all observed v₁ slopes and v₂ at midrapidity are found to be negative and positive, respectively, which is opposite to what is observed at low energy. This can be seen in the 3.0 GeV Au+Au collision results shown in Fig. 14. The early strong partonic expansion leads to positive v₂ with NCQ scaling in high-energy collisions, whereas at 3.0 GeV, both weaker pressure gradients and shadowing of the spectators result in negative v₂ values without scaling. The results from calculations using the hadronic transport models JAM and UrQMD are also shown as colored bands in the figure. By including the baryonic mean field, both the JAM and UrQMD models reproduce the trends for both $\text{d}{v}_1/\text{d}y{|}_{y=0}$ and v₂ for baryons, including protons and Λ. The consistency of the transport models (JAM and UrQMD) with the baryonic mean field for all measured baryons implies that the dominant degrees of freedom at a collision energy of 3.0 GeV are from the interacting baryons. The signatures for the transition from partonic dominance to hadronic and then to baryonic dominance regions have been discussed in [4, 105, 32, 93] for the ratios of K⁺/π⁺ and net-particle v₁ slopes, respectively. The data from 3.0 GeV Au+Au collisions clearly reveal that baryonic interactions dictate the collision dynamics.

Fig. 14Full-screen View

(Color online) Collision energy dependence of the directed flow slope dv₁/dy|_y=0 for p, Λ, charged πs, and kaons (including K^± and ${\text{K}}_S^0$), ϕ, and Ξ^- (top panel) [105, 93]. The bottom panel shows the elliptic flow v₂ for protons and πs from heavy ion collisions [91, 92]. Statistical and systematic uncertainties are shown as bars and gray bands, respectively. The JAM and UrQMD results are displayed as colored bands: golden, red, and blue bands represent the JAM mean-field, UrQMD mean-field, and cascade mode, respectively

The results of collectivity, EOS, and phase structure are closely connected. By comparing measurements with calculations, the parameters of the EOS for each collision can be readily extracted [108, 109]. As an example, Fig. 15 shows the ratio of shear viscosity to entropy as a function of scaled temperature [110]. In the left panel, chemical freeze-out temperature from each energy [32] is used and normalized to that from Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$. In the high energy limit, $\sqrt{s_{\text{NN}}}=39-200\text{ GeV}$, the ratio, namely the quantum limit, reaches unity, implying that the medium created in such collisions is dominated by partonic interactions with a minimum value of 4π η/s. At lower collision energies, on the other hand, hadronic interactions become dominant, and the medium shows a rapid increase in the viscosity-to-entropy ratio. The right panel shows the temperature evolution of the shear viscosity-to-entropy ratio as a function of the scaled temperature T/T_C. Here, T_C represents the critical temperature in the calculation [108, 109]. The entire curve is extracted from the experimental results of R_AA and v₂ from Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$. The observed V-shaped feature is quite similar to what is described in [111] for a system dominated by electromagnetic interactions. The phase transition is universal and independent of the degrees of freedom of the medium under study. This unique feature is a clear evidence of the crossover transition in strong interactions. For a comprehensive discussion on the shear viscosity and phase transition in nucleon and quark levels, the reader is referred to a recent review [112].

Fig. 15Full-screen View

(Color online) The effective values of shear viscosity-to-entropy ratio, 4π η/s, shown as a function of the scaled temperature. The horizontal dashed line indicates the quantum lower limit. Left panel: 4π η/s extracted from the energy dependence of the measured v₂ [107] and v₃ [102], shown as the scaled chemical freeze-out temperature T_ch/T_ch(200 GeV). Right panel: temperature evolution of 4π η/s, extracted from Bayesian analyses [108, 109]

Chirality

Quark interactions with topological gluon configurations can induce chirality imbalance and local parity violation in QCD [113-115]. In relativistic heavy ion collisions, this can lead to observable electric charge separation along the direction of the strong magnetic field produced primarily by spectator protons [116-118]. This is called the chiral magnetic effect (CME). A CME-induced charge separation, if observed, would confirm a fundamental property of QCD. Measurements of electric charge separations can provide a means for studying non-trivial QCD topological structures and are therefore of paramount importance. Extensive theoretical and experimental efforts have been devoted to the search for CME [118-120].

The commonly used observable to measure charge separation is the three-point correlator difference [124], $\Delta {\gamma }_{\text{OS}}-{\gamma }_{\text{SS}}$. Here $\gamma =\langle \cos \left(\alpha +\beta -2{\psi }_2\right)\rangle$, α, and β are the azimuthal angles of the two charged particles and ψ₂ is that of the second-order harmonic plane; γ_OS denotes the γ of opposite electric charge sign (OS), and γ_SS denotes that of the same-sign pairs (SS). The first γ measurements were carried out by the STAR collaboration in Au+Au collisions at top RHIC energy in 2009 [121], which resulted in significant Δγ observations. Further measurements were made at lower RHIC energies by STAR [122] and at higher LHC energies by ALICE [123]. Fig. 16 shows the γ_OS and γ_SS correlators as a function of the collision centrality in Au+Au collisions at $\sqrt{s_{\text{NN}}}=7.7-200\text{ GeV}$ at RHIC and in Pb+Pb collisions at 2.76 TeV at LHC. At high collision energies, charge-dependent signals are observed, and γ_OS is larger than γ_SS. The difference between γ_OS and γ_SS, i.e. Δγ, decreases with increasing centrality, which is consistent with the expectation of a magnetic field strength decrease with increasing centrality. At the low collision energy of $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$, the difference between γ_OS and γ_SS disappears, which is consistent with the disappearance of the CME in the presumably hadronic dominant interactions at this energy. Thus, these results are qualitatively consistent with the CME expectation.

Fig. 16Full-screen View

(Color online) γ correlators as functions of centrality in Au+Au collisions at $\sqrt{s_{\text{NN}}}=7.7-200\text{ GeV}$ from STAR [121, 122] and in Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ from ALICE [123]

One of the difficulties in interpreting the positive Δγ is whether the CME is the major charge-dependent background contribution to the observable [125-127], such as those from resonance decays and jets. The Δγ variable is ambiguous between differentiating an OS pair from the CME back-to-back perpendicular to ψ₂ and an OS pair from a resonance decay along ψ₂. More resonances are produced along the ψ₂ than perpendicular to it, with the relative difference quantified by the elliptical anisotropy parameter v₂ of the resonances. (Jet correlations also exhibit azimuthal anisotropy because of jet quenching effects in heavy ion collisions [128].) The CME background arises from the coupling of this elliptical anisotropy and genuine particle correlations from resonance decays and jets, among others. Calculations using the blast wave parameterizations of the measured particle production data can essentially reproduce the entirety of the measured γ correlations [127].

The CME and the v₂-related background are driven by different physics: the CME is sensitive to the magnetic field which is mostly perpendicular to the spectator plane, whereas the v₂-related background is connected to the participant plane. In non-central heavy ion collisions, the participant plane is generally aligned with the reaction plane, and the Δγ measurement is thus entangled by the two contributions of possible CME and v₂-induced background. In small-system p+A or d+A collisions, the participant plane is determined purely by geometry fluctuations that are uncorrelated with the magnetic field direction [129]. As a result, any CME signal would average to zero in small-system collisions. Background sources, on the other hand, contribute similarly to small-system collisions and heavy ion collisions. In Fig. 17 (left), the first Δγ measurements in small system p+A collisions from CMS [129] are displayed. Within a margin of uncertainty, the results of p+Pb and Pb+Pb collisions exhibit the same magnitude and trend as a function of multiplicity. Fig. 17 (right) shows the Δγ measurements in small system p/d+A collisions from STAR [130]. The trends of the magnitudes are similar, decreasing with increasing multiplicity. These results indicate that strong correlations are present in small systems, contributing to the γ correlators. The nature of these correlations may be attributable to genuine three-particle correlations, which would explain the peripheral heavy ion data but would be insufficient for mid-central heavy ion data as they are strongly diluted by event multiplicity. The nature of some of the correlations may be from flow as there are indications of collective flow in these small systems [131, 132], especially at LHC energies [133]. Nevertheless, the small system results suggest the complex nature of the backgrounds which must be rigorously removed before addressing the important physics of the chiral magnetic effect.

Fig. 17Full-screen View

(Color online) The Δγ correlators as functions of multiplicity in p+Pb and Pb+Pb collisions from the LHC-compact moon solenoid (CMS) [129] (left) and in p/d+Au and Au+Au collisions from RHIC-STAR [130] (right)

Since the major background is induced by v₂, examining the Δγ observable with varying v₂ while holding the expected CME signal constant, is of interest. The event shape engineering (ESE) method is performed based on the magnitude of the flow vector to possibly access the geometry of the initial participant. The ESE selection of events is not expected to affect the magnetic field when restricted to a specified narrow centrality. The different dependencies of the CME signal and background on v₂ (q₂) can possibly be used to disentangle the CME signal from the background. Using the ESE method, the ALICE experiment showed that the CME fraction in the measured Δγ is consistent with zero [134].

Examining the Δγ observable with varying magnetic field while keeping the v₂ relatively constant is also of interest. To gauge the differences in the magnetic field relative to v₂, isobaric and U+U collisions have been proposed [138]. Isobaric collisions were proposed to study two systems with similar v₂ but different magnetic field strengths [138], such as ${}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}$, which have the same mass number but different charge (proton) number. One would thus expect that v₂ would be very similar at midrapidity in ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions, but the magnetic field, proportional to the nuclei electric charge, could vary by 10%. The variation of the magnetic field strength between ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions provides an ideal way to disentangle the signal of the chiral magnetic effect from v₂ related background, as v₂ related backgrounds are expected to be very similar between these two systems.

Figure 18 shows the ratio of Δγ/v₂ in Ru+Ru over Zr+Zr collisions from the isobar analysis [135-137], as an observable. The CME-sensitive observable ratios lie below unity leading to the conclusion that no predefined CME signatures, such as a larger-than-unity Ru+Ru over Zr+Zr ratio of Δγ/v₂, are observed in this blind analysis. This is rather counterintuitive at a first glance but can be explained using nuclear structure considerations. In fact, the ⁹⁶Zr nucleus was predicted to be larger than ⁹⁶Ru because of its thicker neutron skin, resulting in a slightly smaller energy density and fewer particles being produced in Zr+Zr than in Ru+Ru collisions [139-141]. The larger ⁹⁶Zr nucleus also provides smaller eccentricity at a given centrality and thus smaller v₂ [139, 140]. Although the non-identical v₂ is properly considered in the blind analysis observable Δγ/v₂, non-identical event multiplicities are not. After properly factoring in the multiplicity, the isobar ratios of NΔγ/v₂ from various analyses shown in Fig. 18 indicate a positive signal of a few standard deviations [135, 142]. However, non-flow contamination exists in the Δγ/v₂ ratio variable [143]. One such contamination is the aforementioned genuine three-particle correlations because Δγ is measured by the three-particle correlator in STAR TPC. Another contamination is attributed to the fact that two-particle v₂ cumulant measurements are contaminated by non-flow correlations and such v₂ values are used to compute the Δγ from the three-particle correlator measurement. Rigorous studies of non-flow contamination have been conducted in post-blind analyses, and improved background baselines are derived [136, 137]. Figure 18 shows the measured isobar ratios of Δγ/v₂ from the blind analysis along with the estimated background baselines from the post-blind analysis. The results show that the isobar ratios are consistent with the baselines, indicating that no statistically significant CME signals have been observed in the isobar data.

Fig. 18Full-screen View

(Color online) Compilation of results from isobar analysis. Only contrasting results from the two isobar systems are shown. The results are expressed in terms of Ru+Ru over Zr+Zr collision measurement ratios. Solid dark symbols indicate CME-sensitive measures, whereas open light symbols indicate their counterpart insensitive to CME. The vertical lines indicate statistical uncertainties whereas boxes indicate systematic uncertainties. The colors in the background are intended to separate different types of measures. The fact that CME-sensitive observable ratios lie below unity leads to the conclusion that no predefined CME signatures are observed in this blind analysis [135]. The estimated background baselines from non-flow contamination for the four cumulant measurements of the isobar Δγ/v₂ ratios are indicated by horizontal bars (central values) and shaded areas (total uncertainties, the quadratic sum of the statistical and systematic uncertainties on the background baseline estimates) [136, 137]

The STAR isobar data, without any clear evidence for a possible CME-related signal difference possibly arising from the charge difference (44 in Ru versus 40 in Zr), have provided important lessons in the experimental search for CME. First, the difference in nuclear shape and/or neutron skin between isobaric nuclei can induce percent-level background variations, which cannot be easily estimated using theoretical calculations or controlled with experimental constraints. Thus, searches for small differences in the CME signal arising from the magnetic field variation in isobar collisions is extremely challenging. Second, the strength of the magnetic field plays a critical role in the CME signal; therefore, larger nuclei would be preferable in the search for a possible CME-induced signal in Δγ correlations. Third, a better understanding of the background sources in the Δγ correlator is required to suppress this background from elliptic flow and non-flow correlations.

The major background source in the CME observable Δγ is induced by elliptic flow (v₂). The original event shape engineering approach [144, 134, 145] uses particles from separate rapidity or pseudorapidity regions to define event classes. This approach can select event shapes sensitive to the eccentricity of the initial overlapping participants and the corresponding geometrical fluctuations. However, for particles of interest used for measuring the CME-sensitive observable Δγ in a different rapidity region, the event-by-event v₂ background has contributions from both eccentricity and particle emission pattern fluctuations. Petersen and Muller [146] pointed out that emission pattern fluctuations dominate the event-by-event v₂ fluctuations. Recently, Xu et al. proposed a novel event shape selection (ESS) approach to suppress the background in the CME Δγ measurement [147]. They found that to suppress the apparent flow-induced background in Δγ, the combined event-by-event information from eccentricity and emission pattern fluctuations from particles of interest should be used to select azimuthally round shape events for correlator measurements. With this ESS approach, the suppression of the flow-related background becomes possible. Using a multi-phase transport (AMPT) [148] and anomalous viscous fluid dynamics (AVFD) [149] model simulations, Xu et al. [147] showed that the most effective ESS approach is to use particle pairs to construct the event shape variable, thereby forming event shape classes for the CME sensitive correlator to calculate the zero elliptic flow at the limit for particles of interest. This is consistent with the expectation that the background in Δγ has significant contributions from particle pair emissions coupled with elliptic flow.

The RHIC BES-II also provides a unique venue for the CME search, covering the center of mass energies from 7.7 GeV to 27 GeV. At these beam energies, the STAR event plane detector (EPD), added during the BES-II program, can register spectator protons from the colliding beams. This capability allows an accurate estimation of the reaction plane, enhancing the sensitivity to the magnetic field direction and suppressing non-flow contributions to the background. For Au+Au collisions at the top RHIC energy, spectator neutrons may be detected by the zero-degree calorimeter (ZDC) although the corresponding event plane resolution is not as good as that in the BES-II data. Theoretical calculations expect the initial magnetic field to be smaller in Au+Au collisions from BES-II than that from the top RHIC energy. However, the dynamics of the QGP formation and time evolution of the magnetic field in QGP as a function of collision energy have not been fully understood. Recent STAR measurements of the deflection of charged particles by the magnetic field in heavy ion collisions indicate significant imprints of magnetic-field effects at these BES-II energies [150]. The STAR collaboration reported preliminary results on the CME search from the RHIC BES-II data at the 2023 Quark Matter Conference, demonstrating a promising approach for focusing on Au+Au collisions using an innovative experimental technique for background suppression [151].

As aforementioned, the Δγ measurement in heavy ion collisions is entangled by two contributions, one from the CME and the other from the v₂ induced background. These are sensitive to differences in planes, allowing the measurement of Δγ. The background is related to v₂, as determined by the participant geometry, and therefore is the largest with respect to the participant plane (ψ_PP). The CME-driven charge separation is along the magnetic field direction (ψB), unlike ψ_PP. The ψB and ψ_PP are generally correlated with the impact parameter direction, ψ_RP, and therefore correlated among themselves. While the magnetic field is mainly produced by spectator protons, their positions fluctuate, thus ψB is not always perpendicular to ψ_RP. The position fluctuations of participant nucleons and spectator protons are independent, whereby ψ_PP and ψB fluctuate independently about ψ_RP. Notably, a new approach has been proposed to measure Δγ with respect to ψ_SP and ψ_PP to disentangle the CME signal from the v₂ background [152, 153]. This is exploited by STAR in measuring Δγ with respect to the first-order harmonic plane from the ZDC and second-order harmonic plane from the TPC. Because the former aligns better with the spectator plane and the latter aligns better with the participant plane, these measurements contain different amounts of sensitive flow backgrounds and magnetic field-sensitive CME signal in the harmonic plane, enabling the extraction of a possible CME.

STAR reported such measurements in Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [154], as shown in Fig. 19. The charge separation, with the flow background removed, is consistent with zero in peripheral collisions. In mid-central collisions, on the other hand, an intriguing indication of finite CME signals is observed on the order of 1–3σ standard deviations.

Fig. 19Full-screen View

(Color online) The flow-background removed CME signal fraction $\langle f_{\text{CME}}\rangle$ in 50–80% (open markers) and 20–50% (solid markers) centrality Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [154]. Results are shown for the full-event (FE) analysis method with two $p_{\perp }$ ranges and sub-event (SE) analysis method with two Δη gaps. Error bars show statistical uncertainties; the caps indicate systematic uncertainties

In RHIC 2023-2025, STAR is expected to collect approximately 20 B events, which is about a factor of 10 more compared to the data used for Fig. 19. More precise results are expected in the near future. New analyses utilizing event-shape-engineering with particle pair anisotropy and invariant mass are ongoing, and results are expected soon.

Criticality

In high energy nuclear collisions where the baryon density is vanishingly small, the transition from QGP to hadronic matter is a smooth crossover [155]. At finite density and lower temperatures, the transition is speculated to be first-order, with an associated phase boundary. The point that connects the smooth crossover and first-order phase boundary is the QCD critical point [156]. Since 2010, RHIC has conducted two rounds of beam energy scan campaigns primarily aimed at investigating the QCD critical point. The BES programs cover an energy range from $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ to 3.0 GeV, corresponding to a baryonic chemical potential of 20 MeV≤μ_B≤750 MeV. As of the summer of 2022, both BES-I and BES-II have been completed.

High-order cumulants of protons and net-protons (event-by-event number: net-p = $\text{p}-\overline{\text{p}}$) distributions are used in the search for the QCD critical point [157] because of their high sensitivity to correlation length. The experimental results shown as a function of the collision energy are depicted in Fig. 20. Overall, the ratios of C₄/C₂ for net-protons from collider mode ($\sqrt{s_{\text{NN}}}\ge 7.7\text{ GeV}$) [158, 38] and protons from the fixed-target mode decrease as collision energy decreases as dictated by the baryon number conservation. Both the hadronic resonance gas model and hadronic transport model UrQMD [76, 77] calculations reproduce this trend. As a function of collision energy, a rise and then fall of the net-proton C₄/C₂ (or $\kappa {\sigma }^2$) is expected to indicate the critical behavior near the critical point in the QCD phase diagram. While results of C₄/C₂ ratios from BES-I have shown dip like energy dependence around 20 GeV, the statistics at lower collision energies are too poor to draw any conclusions regarding this prediction. Note that at low energies, or equivalently, in the high baryon density region, both the HADES ($\sqrt{s_{\text{NN}}}=2.4\text{ GeV}$) and STAR ($\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ [159, 160]) high moment proton results are below the Poisson baseline, and the non-critical hadronic transport model calculations reproduce the data at the high baryon density region. This implies that this energy regime is dominated by hadronic interactions. To look for an oscillation pattern in the energy dependence of the ratio of C₄/C₂, the second phase of the beam energy scan (BES-II) was conducted at RHIC. The BES-II data analyses are under way.

Fig. 20Full-screen View

(Color online) Collision energy dependence of the ratios of cumulants, C₄/C₂, for proton (squares) and net-proton (red circles) for top 0–5% Au+Au collisions at RHIC [38]. The points for protons are shifted horizontally for clarity. The new result for proton from $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ collisions is shown as a filled square. High acceptance dielectron spectrometer (HADES) data of $\sqrt{s_{\text{NN}}}=2.4\text{ GeV}$ 0–10% collisions is also shown. Results from HRG and transport model UrQMD [76, 77] are shown

In addition, it is predicted that density fluctuations near the QCD critical point can be probed via an intermittency analysis in relativistic heavy ion collisions [162, 163]. Figure 21 shows the energy dependence of the scaling exponent (ν) for identified charged hadrons in Au+Au collisions for two different collision centralities (0–5% and 10–40%) [161]. In the most central collisions, ν exhibits a non-monotonic behavior as a function of collision energy, reaching a minimum around $\sqrt{s_{{}_{\text{NN}}}}=20-30\text{ GeV}$. In contrast, for 10–40% central collisions, ν remains approximately constant with increasing $\sqrt{s_{{}_{\text{NN}}}}$. The observed non-monotonic energy dependence of ν in the most central collisions can indicate the density fluctuations induced by the QCD critical point. However, at $\sqrt{s_{{}_{\text{NN}}}}\le 11.5\text{ GeV}$, large systematic and statistical uncertainties exist for ν. Higher statistics data are required from the BES-II program to confirm this energy dependence. The measured value of ν is significantly smaller than the theoretical predictions of ν= 1.30 from the Ginzburg-Landau (GL) theory and ν= 1.0 from the 2D Ising model. These theoretical values are derived from calculations over the entire phase space without constraints on acceptance, whereas the experimental measurements are limited to the available transverse momentum space; ν is anticipated to increase if measured over the entire phase space, particularly when including higher p_T regions. Therefore, theoretical calculations that consider a reduced transverse momentum phase space and equivalent experimental acceptance, are required to understand the measured scaling exponent. The transport-based UrQMD model is unable to calculate ν because of the absence of power-law scaling of $\Delta F_q\left(M\right)\propto \Delta F_2{\left(M\right)}^{{\beta }_q}$. Consequently, models that exhibit such power-law scaling are required to produce a non-critical baseline for comparison with experimental data.

Fig. 21Full-screen View

(Color online) Energy dependence of the scaling exponent (ν) for identified charged hadrons (h^±) in Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}=7.7-200\text{ GeV}$ [161]. Red circles and blue squares represent ν in the most central collisions (0–5%) and mid-central collisions (10–40%), respectively. The statistical and systematic errors are shown in bars and brackets, respectively

Recently, a study was conducted on the information entropy [164] of the net-proton multiplicity distribution using the hybrid model of ideal fluid dynamics plus UrQMD [165]. The ratios of the net-proton information entropies of the UrQMD result with EoS:CH (chiral+hadronic gas EoS with first-order transition and critical endpoint) were compared with those of UrQMD with EoS:BM (bag model EoS with strong first-order phase transition between QGP and hadronic phase), STAR experimental data, and UrQMD results without hydrodynamic EOS. The comparisons show that the STAR experimental data extracted from [38] display an enhancement in collision at 20 GeV with respect to the baseline entropy without hydrodynamics, which is consistent with the minimum $\kappa {\sigma }^2$ value reported in [158]. On the other hand, the UrQMD simulations with the EoS:BM and EoS:CH EOS also show slightly pronounced enhancements, at energies higher than approzimately 30 GeV, consistent with recent observations of $N_{\text{t}}{N}_{\text{p}}/N_{\text{d}}^2$ [166] and the analysis of the intermittency scaling exponent [161], with a peak or dip around $\sqrt{s_{\text{NN}}}=20-30\text{ GeV}$, which could indicate the presence of a critical end point (CEP) in these collisions. Thus, the information entropy can also be regarded as an indication of an alternative observable for studying the QGP phase transition.

Global Polarization of QCD Matter

In non-central relativistic heavy ion collisions, huge orbital angular momenta (OAM) and vorticity fields are produced in QGP [167]. These can lead to hadron polarization and spin alignment along the direction of the system OAM through spin-orbit couplings [168-170] or spin-vorticity couplings [171, 172], a phenomenon known as global polarization. Such polarization phenomena in relativistic heavy ion collisions possess some unique features which are different from those of conventional observations. For example, the measurement is not mediated by a magnetic field, as in the well-known Barnett effect of magnetization due to rotation [173]. The global spin polarization of particles is directly observed in relativistic heavy ion collisions, which is not possible in ordinary matter. Also, QGP at very high energies is almost neutral by charge conjugation. However, if it were precisely neutral, the observation of polarization by magnetization would be impossible because particles and antiparticles would have opposite magnetic moments. In fact, Λ and $\overline{\text{Λ}}$ in relativistic heavy ion collisions at high energy have almost the same mean polarization, suggesting that polarization is a strong interaction driven phenomenon. If the electromagnetic field were responsible for this effect, the signs of its mean spin vector components would be opposite. Hence, for non-relativistic matter, the fact that it is impossible to resolve polarization by rotation and magnetization, is at the very heart of the Barnett effect [173] and its reverse Einstein-de Hass effect [174]. In relativistic matter, because of the existence of antiparticles, the rotation and magnetization effects can be distinguished, and QGP is the first relativistic system through which this distinction has been observed [175].

The global polarization of hyperons can be determined from the angular distribution of hyperon decay products in hyperon’s rest frame with respect to the system OAM: $\frac{\text{d}N}{\text{d}\cos {\theta }^{*}}\propto 1+{\alpha }_{\text{H}}{P}_{\text{H}}\cos {\theta }^{*}\mathrm{,}$ (1) where ${\alpha }_{\text{H}}$ is the hyperon decay parameter; P_H is the hyperon polarization; and θ^* is the angle between the polarization vector and direction of the daughter nucleon momentum in the hyperon rest frame. Since the system OAM is perpendicular to the reaction plane, the global polarization can be measured via the distribution of the azimuthal angle of the hyperon decay nucleon in the hyperon rest frame with respect to the reaction plane. The reaction plane is defined by the direction of the incoming nuclei (beam direction) and impact parameter vector ($\widehat{b}$) [176]. The reader is referred to [177, 178] for the analysis details and focus on the results presented here.

Figure 22 shows the first measurement of P_H at $\sqrt{s_{\text{NN}}}=62.4\text{ GeV}$ and 200 GeV in the STAR experiment, which are consistent with zero [179]. The later STAR measurements at $\sqrt{s_{\text{NN}}}=3\text{ GeV}$, 7.7 GeV–39 GeV [175, 181], with higher statistics at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [177], indicate statistically significant global polarization P_H>0, whereas high-statistics ALICE measurements at $\sqrt{s_{\text{NN}}}=2.76\text{ TeV}$ and 5.02 TeV demonstrate a P_H consistent with zero at the LHC energies [180]. The P_H is observed to increase with collision centrality, in agreement with the increasing system OAM from central to peripheral collisions. Fig. 22 also shows a measurement of Λ polarization shifting to lower energies of Au+Au collisions at $\sqrt{s_{\text{NN}}}=2.4\text{ GeV}$ and Ag+Ag collisions at $\sqrt{s_{\text{NN}}}=2.55\text{ GeV}$ in the HADES experiment [182]. An increasing trend of P_H with decreasing $\sqrt{s_{\text{NN}}}$ is observed. The collision energy dependence of experimental data can be reasonably described by theoretical calculations, as displayed in the figure, including hydrodynamic [184, 187], transport [186], and chiral kinetic simulations [185]. Some models also predict that P_H would vanish at $\sqrt{s_{\text{NN}}}=2m_{\text{N}}$, and thus P_H may peak around 3.0 GeV [187, 188]. Investigating the dependencies of P_H on the hyperon transverse momentum p_T and rapidity y is also interesting, as different models display opposite trends for the high rapidity region [189, 185]. The available measurements mostly cover midrapidity, and the observed P_H versus p_T and y is constant within a margin of uncertainty. Future measurements in the large rapidity region would be of special interest, particularly after the STAR forward detector upgrade. There are also discussions on the collision system dependence of P_H, for example in smaller colliding systems [190]. Recently, the STAR experiment measured the Λ global polarization in isobar Ru+Ru and Zr+Zr collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ and observed Λ polarization along the beam direction relative to the second and third harmonic event planes originating from the local vorticity [191].

Fig. 22Full-screen View

(Color online) Global Λ and $\overline{\text{Λ}}$ polarization as a function of $\sqrt{s_{\text{NN}}}$ in mid-central heavy ion collisions [175, 177, 179-183]. For clarity, data points of the same collision energy from updated measurements are slightly shifted along the x-axis. Calculations using the hybrid (UrQMD+vHLLE) [184], chiral-kinetic transport (Chiral kinetic) [185], and multi-phase transport models [186] are comparable to the higher $\sqrt{s_{\text{NN}}}$ data, whereas the hydrodynamics 3-fluid model with different equations of state predicts a sharply rising P_Λ at lower $\sqrt{s_{\text{NN}}}$ [187]

All particles and antiparticles of the same spin should have the same global polarization assuming that OAM is the only driving source of polarization. A difference can arise from the effects of the initial magnetic field, given the fact that particles and their antiparticles have opposite magnetic moments. In addition, different particles can be produced at different times or regions as the system freezes out, or through meson-baryon interactions. The measurements of Λ and $\overline{\text{Λ}}$ polarization in $\sqrt{s_{\text{NN}}}=7.7-39\text{ GeV}$ show no difference within current uncertainties. Therefore, to establish the global nature of the polarization, it is very important to measure the polarization of different particles, and if possible, particles of different spins.

Global polarization leaves its imprint on vector mesons such as ϕ(1020) and K^*0(892). Unlike Λ and $\overline{\text{Λ}}$ hyperons that can undergo weak decay with parity violation, the polarization of vector mesons cannot be directly measured since they mainly decay through the strong interaction, in which parity is conserved. Nevertheless, the spin state of a spin-1 vector meson can be described by a 3× 3 spin density matrix with unit trace [192]. The diagonal elements of this matrix, namely, ${\rho }_{11}\mathrm{,}{\rho }_{00}$ and ${\rho }_{-\text{1}-\text{1}}$, are probabilities of the spin component and take on the values of 1, 0, and -1, respectively, along a quantization axis, which is selected for the projection of OAM to have well-determined values. When the three spin states have equal probability of being occupied, all three elements are 1/3, and there is no spin alignment. If ${\rho }_{00}\ne 1/3$, the spin of the vector meson is aligned with the spin quantization direction. For a vector meson decaying into two spin-0 daughters, the angular distribution of one of its decay products in the vector meson rest frame can be written as $\frac{\text{d}N}{\text{dcos}{\theta }^{*}}\propto \left(1-{\rho }_{00}\right)+\left(3{\rho }_{00}-1\right){\text{cos}}^2{\theta }^{*}\mathrm{,}$ (2) where θ^* is defined as in Eq. (1), that is, by the polar angle between the quantization axis and momentum direction of one of the decay products. For our study of global spin alignment, the quantization axis was chosen to be the direction of the system OAM, which is perpendicular to the reaction plane. By fitting the angular distribution of decay products using Eq.(2), one can infer the ρ₀₀ value.

The search for global spin alignment of ϕ(1020) and K^*0(892) mesons for Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ started in parallel with the search for Λ polarization. Due to limited statistics at that time, no significant results were reported [193].

Figure 23 presents the ϕ(1020) meson spin alignments in Au+Au collisions at beam energies between $\sqrt{s_{\text{NN}}}=11.5\text{GeV}$ and 200 GeV [197]. The STAR measurements presented in Fig. 23 are for centralities between 20% and 60% where a maximum OAM of the collision system is expected. The quantization axis is the normal to the 2nd-order event plane [176] (a proxy for the reaction plane), determined using STAR charged particle information. The ϕ-meson results are presented for transverse momentum 1.2 GeV/c < p_T < 5.4 GeV/c, and ρ₀₀ which are significantly above 1/3 for collision energies of 62 GeV and below for this species, indicating finite global spin alignment. The ρ₀₀ for ϕ mesons, averaged over beam energies of 62 GeV and below, is 0.3512±0.0017 (stat.) ±0.0017 (syst.). Taking the total uncertainty as the quadrature sum of statistical and systematical uncertainties, our results indicate that the ϕ-meson ρ₀₀ is above 1/3 with a significance of 7.4 σ [197].

Fig. 23Full-screen View

(Color online) Global spin alignment of ϕ and K^*0 vector mesons in heavy ion collisions. The measured matrix element ρ₀₀ as a function of beam energy for the ϕ and K^*0 vector mesons within the indicated windows of centrality, transverse momentum (p_T), and rapidity (y). The open symbols indicate the ALICE results [194] for Pb+Pb collisions at 2.76 TeV and p_T values of 2.0 GeV/c and 1.4 GeV/c for the ϕ and K^*0 mesons, respectively. The blue solid curve is a fit to the data in the range of $\sqrt{s_{\text{NN}}}=19.6\text{ GeV}$ to 200 GeV, based on a theoretical calculation using the ϕ-meson field [195]. Parameter sensitivity of ρ₀₀ to the ϕ-meson field is shown in [196]. The blue dashed line is an extension of the solid curve for the fitted parameter $G_s^{\left(y\right)}$. The black dashed line represents ${\rho }_{00}=1/3$

Figure 23 also presents the beam-energy dependence of ρ₀₀ for K^*0 within 1.0 GeV/c <p _T<5.0 GeV/c. We observe that the ρ₀₀ for K^*0 is largely consistent with 1/3, in marked contrast to the case for ϕ. The ρ₀₀ for K^*0, averaged over beam energies of 54.4 GeV and below, is 0.3356 ± 0.0034 (stat.) ± 0.0043 (syst.), and the deviation from 1/3 has a ～0.42σ significance [197]. Measurements from the ALICE collaboration for Pb+Pb collisions at $\sqrt{s_{\text{NN}}}$ = 2.76 TeV [194], taken from the data points [194] closest to the mean p_T for the range of 1.0 GeV/c <p _T<5.0 GeV/c, are also shown for comparison in Fig. 23. They are consistent with 1/3 considering the large statistical uncertainties.

According to the quark coalescence for hadron production in heavy ion collisions, the Λ polarization depends linearly on the quark polarization, whereas the vector meson polarization displays a quadratic dependence [168, 169]. One would therefore expect the polarization for ϕ to be smaller than the one measured for Λ. However, the measured ρ₀₀ of ϕ is orders of magnitude larger than the expected from the same vorticity that caused the measured Λ and $\overline{\text{Λ}}$ polarizations in the same collisions. Contributions from electromagnetic fields and other possible conventional mechanisms are also orders of magnitude smaller than expected from the data [195, 198-200]. However, a new mechanism of vector meson spin alignment due to a strong force field is capable of describing both the ρ₀₀ of ϕ and K^*0 [196], pointing out that the difference of Λ polarization and vector meson spin alignment can be understood as follows. The Λ polarization provides information on the mean values of quark polarization, whereas the ρ₀₀ provides information on the correlation of quark and antiquark polarization inside the vector meson [201]. Thus measurements of vector meson spin alignment provide a novel way for probing the quark spin correlations. This information may also be accessible via the measurements of hyperon-hyperon and hyperon-antihyperon spin correlations [201]. These different scenarios open an exciting discovery potential for the spin polarization measurements. For example, the strong force correlation may be expected to provide new information about the short distance structure of QGP and the nature of QCD phase diagram [202, 203].

Light Cluster Formation

Light nuclei and hypernuclei are loosely bound objects of nucleons and hyperons with binding energies of several MeV. Their formation in heavy ion collisions provides important information on the properties of nuclear matter at high densities and temperatures, such as the nucleon-nucleon/hyperon interactions. The equation of state may offer insights into the inner structure of compact stars.

The production of light nuclei in relativistic nucleus-nucleus collisions has been studied since the early 1960s [204], and their production mechanisms are still being debated [205-207]. The thermal/statistical and nucleon coalescence models are two widely recognized and effective methods for explaining the production of light nuclei in high-energy heavy ion collisions. In the thermal model, the formation of light nuclei is similar to that of hadrons, with the yields calculated based on particle masses and the thermodynamic properties near the chemical freeze-out of the collision system [205, 206]. The coalescence model assumes that light nuclei emerge through a combination of nucleons coming close to each other at the time of kinetic freeze-out [208-210].

Based on the coalescence model, the compound yield ratio $N_{\text{t}}\times N_{\text{p}}/N_{\text{d}}^2$ of tritons (N_t), deuterons (N_d), and protons (N_p) is predicted to be sensitive to the neutron density fluctuations, making it a promising observable to search for in the signature of CEP and/or a first-order phase transition in heavy ion collisions [212-215]. The expected signature of CEP is its non-monotonic variation as a function of collision energy.

Figure 24 shows the charged-particle multiplicity $\text{d}{N}_{\text{ch}}/\text{d}\eta$ (|η| < 0.5) dependence of the yield ratio $N_{\text{t}}\times N_{\text{p}}/N_{\text{d}}^2$ in Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}=7.7-200\text{ GeV}$, combining all centrality bins [166]. The yield ratio exhibits a charge-particle multiplicity scaling, regardless of collision energy and centrality. The shaded bands are the corresponding results from the calculations of hadronic transport AMPT and MUSIC+UrQMD hybrid models, in which neither critical point nor first-order phase transition is included. These two models are employed to generate the nucleon phase space at kinetic freeze-out, with light nuclei forming via nucleon coalescence. The overall trend of the experimental data is well described by the model calculations. The light blue dashed line is the result calculated from the thermal model at chemical freeze-out ($T_{\text{ch}}=157\text{ MeV}$ at 200 GeV) for central Au+Au collisions, which overestimates the experimental data by more than a factor of two at $\text{d}{N}_{\text{ch}}/\text{d}\eta \sim 600$ which could be due to the effects of hadronic rescattering during hadronic expansion. The black dot-dashed line is a fit to the data based on the coalescence model. The lower panel of Fig. 24 shows that most of the measurements are within a significance of 2 σ from the coalescence baseline, except for the enhancements observed for the yield ratios from the 0-10% most central Au+Au collisions at 19.6 GeV and 27 GeV with a significance of 2.3 σ and 3.4 σ, respectively. Notably, in the net-proton higher moments and charged particle intermittency measurements, non-monotonic behaviors are observed around a collision energy of $\sqrt{s_{{}_{\text{NN}}}}=20\text{ GeV}$. Further studies on dynamical modeling of heavy ion collisions based on a realistic equation of state are required to confirm that the enhancements are due to large baryon density fluctuations near the critical point. These systematic measurements of triton yields and yield ratios over a broad energy range provide important insights into the production dynamics of light nuclei, enhancing our understanding of the QCD phase diagram.

Fig. 24Full-screen View

(Color online) The yield ratio $N_{\text{t}}\times N_{\text{p}}/N_{\text{d}}^2$ as a function of charged-particle multiplicity $\text{d}{N}_{\text{ch}}/\text{d}\eta$ (|η|<, 0.5) in Au+Au collisions at $\sqrt{s_{\text{NN}}}=7.7-200\text{ GeV}$ for various collision centralities [166]. The black dot-dashed line denotes the coalescence-inspired fit. The significance of the deviation relative to the fit is shown in the lower panel. The results calculated from the thermal model are shown as the blue long-dashed line. Calculations from AMPT and MUSIC+UrQMD hybrid models are shown as shaded bands [211]

Similar to the number of constituent quark scaling of hadron flow, the light nuclei flow is expected to exhibit an approximate scaling with the mass number A scaling under the coalescence assumption [97] $v_n^A\left(p_{\text{T}}\mathrm{,}y\right)/A\approx v_n^p\left(p_{\text{T}}/A\mathrm{,}y\right)\mathrm{.}$ (3) However, unlike quarks, whose flow cannot be directly measured, both proton and light nuclei flow can be directly measured in experiments to validate the coalescence model. Figure 25 shows the light nucleus v₁ slopes dv₁/dy|_y=0, which are utilized to characterize the strength of v₁, scaled by the atomic mass number, as a function of collision energy from the $\sqrt{s_{\text{NN}}}=3-40\text{ GeV}$ STAR experiment [216-218]. Overall, the $d{v}_1/dy{|}_{y=0}$ decreases monotonically with increasing collision energy for both protons and light nuclei. At $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$, the dv₁/dy|_y=0 follow an approximate scaling with the atomic mass number A. The transport model calculation with a baryon mean field and a coalescence afterburner qualitatively reproduces the measurements for both protons and light nuclei, as indicated by the short lines near the data points. The results indicate that the light nuclei are likely formed via the coalescence of nucleons in $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ Au+Au collisions, with the baryonic interactions dominating the collision dynamics [216]. At $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$, the A scaling still holds for v₁ of the deuteron. However, in moving to higher energies, the dv₁/dy|_y=0 values for protons become negative, while the corresponding value for deuterons is still positive albeit with larger uncertainties [217]. This discrepancy in the scaling behavior of light nuclei dv₁/dy|_y=0 at energies below 7.7 GeV and above 11.5 GeV may indicate a different production mechanism or system evolution, as the QGP is expected to form at higher energies, and interactions occur at the partonic level [92, 91].

Fig. 25Full-screen View

(Color online) Light nucleus v₁ slopes dv₁/dy|_y=0 scaled by the atomic mass number as a function of collision energy in 10–40% mid-central Au+Au collisions [216-218]. For clarity, the data points are shifted horizontally. Results of the JAM model in the mean-field mode plus coalescence calculations are shown as color bars

Hypernuclei are nuclei containing at least one hyperon. As such, they are excellent experimental probes to study the hyperon-nucleon (Y–N) interaction [219, 220], an important ingredient in the EOS of dense nuclear matter [221, 222]. Similar to light nuclei production in heavy ion collisions, statistical thermal hadronization [205] and coalescence models [221] have been proposed to describe hypernuclei formation. Although thermal model calculations primarily depend only on the freeze-out temperature and the baryon-chemical potential, the Y–N interaction plays an important role in the coalescence approach through its influence on the dynamics of hyperon transportation in the nuclear medium, as well as its connection to the coalescence criterion for hypernuclei formation from hyperons and nucleons [221].

Figure 26 shows the ${}_{\text{Λ}}^3\text{H}$ and ${}_{\text{Λ}}^4\text{H}$ midrapidity yields for central Au+Au collisions of $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ in comparison with the measurements at LHC, with the insets showing the $\text{d}N/\text{d}y\times B\mathrm{.}R\mathrm{.}$ as a function of B.R.. We observe that the ${}_{\text{Λ}}^3\text{H}$ yield in Au+Au collisions at $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ is significantly enhanced compared to the yield at LHC, likely driven by the increase in baryon density at low energies. Calculations from the thermal model [205], which adopt the canonical ensemble for strangeness, mandatory at low beam energies, are compared with the experimental data. Interestingly, while the ${}_{\text{Λ}}^3\text{H}$ yields at 3.0 GeV and 2.76 TeV are well described by the model, the ${}_{\text{Λ}}^4\text{H}$ yield is underestimated by approximately a factor of four. Coalescence calculations using the Dubna cascade model (DCM), which is an intra-nuclear cascade model for describing the dynamical stage of the reaction [221], are consistent with the ${}_{\text{Λ}}^3\text{H}$ yield while underestimating the ${}_{\text{Λ}}^4\text{H}$ yield, whereas coalescence using JAM calculations is consistent with both. Note that in the DCM model, the same coalescence parameters are assumed for two hypernuclei, whereas in the JAM model, parameters are tuned separately for ${}_{\text{Λ}}^3\text{H}$ and ${}_{\text{Λ}}^4\text{H}$ to fit the data. The calculated hypernuclei yields are expected to depend on the choice of the coalescence parameters [221]. The recent calculations from parton-hadron-quantum-molecular-dynamics (PHQMD) [224], a microscopic transport model which utilizes a dynamical description of hypernuclei formation, are consistent with the measured yields within a margin of uncertainty. Compared to the JAM model which adopts a baryonic mean field approach, baryonic interactions in PHQMD are modeled by density dependent 2-body baryonic potentials. Meanwhile, the UrQMD-hydro hybrid model overestimates the yields at 3.0 GeV by an order of magnitude [223]. STAR measurements possess the capability to distinguish between different production models and provide new baselines for the strangeness canonical volume in thermal models and coalescence parameters in transport plus coalescence models. Such constraints can be utilized to improve the model estimation of the production of exotic strange matter in the high baryon density region.

Fig. 26Full-screen View

(Color online) Beam energy dependent (a) ${}_{\text{Λ}}^3\text{H}$ and (b) ${}_{\text{Λ}}^4\text{H}$ yields at |y|<0.5 in central heavy ion collisions compared to theoretical model calculations. The data points assume a branching ratio (B.R.) of 25(50)% for ${}_{\text{Λ}}^3\text{H}$ (${}_{\text{Λ}}^4\text{H}$) → ³He (⁴He) + π^-. The insets show their yields at |y|<0.5 times the B.R. as a function of B.R. [223]

The STAR experiment reported the first observation of the v₁ of hypernuclei ${}_{\text{Λ}}^3\text{H}$ and ${}_{\text{Λ}}^4\text{H}$ in 3.0 GeV Au+Au collisions [225], as shown in Fig. 27. The mass dependence of dv₁/dy|_y=0 for Λ and hypernuclei is similar to that of light nuclei, increasing linearly with the particle mass, following a baryon mass number scaling. Although the dv₁/dy|_y=0 values for hypernuclei are systematically lower compared to those for nuclei of equivalent mass numbers, this discrepancy may be attributed to the fact that the dv₁/dy|_y=0 for Λ is lower than that for protons. The calculations using the transport model plus an afterburner qualitatively reproduce the data within uncertainties, suggesting that the hypernuclei are produced via coalescence of hyperon and light nuclei core in such heavy ion collisions. If hypernuclei are formed through the coalescence process, both their v₁ and yield can be affected by the interactions involving hyperons and nucleons (Y-N), which are essential for understanding the inner structure of compact stellar objects. The linear fits to the extracted dv₁/dy|_y=0 in Fig. 27 show comparable slopes considering uncertainties for both light nuclei and hypernuclei, but their central values are slightly different. This difference may originate from the differences in nucleon-nucleon and Y-N interactions. Thus, more precise measurements with increased statistics, especially at high baryon density, will be crucial in elucidating the hypernuclei production mechanisms and hyperon-nucleon interactions in the future.

Fig. 27Full-screen View

(Color online) Mass dependence of light nuclei and hypernuclei v₁ slopes dv₁/dy|_y=0 from the $\sqrt{s_{\text{NN}}}=3\text{ GeV}$ 5%-40% centrality Au+Au collisions [225].The dashed lines are the results of a linear fit to the measured light nuclei and hypernuclei dv₁/dy|_y=0, respectively. The calculations of transport models plus coalescence afterburner are shown as gold and red bars for the JAM model, and blue bars for the UrQMD model

Heavy Flavor Hadron Production

Heavy flavor hadrons are hadrons with at least one constituent heavy flavor quark. They are penetrating probes of QGP. Heavy flavor quarks are predominantly produced through initial hard scattering processes in heavy ion collisions owing to their large masses. These initial hard processes occur before the formation of QGP. Consequently, heavy flavor quarks experience the entire evolution of QGP created in heavy ion collisions. Heavy flavor quarks interact with the deconfined quarks, mainly light flavor quarks, and gluons when they transit QGP and approach thermalization. Their thermal relaxation time is expected to be comparable to or longer than the lifetime of the QGP created in heavy ion collisions. Heavy flavor quarks may acquire collectivity from the collectively expanding hot medium. The collectivity of heavy flavor quarks is sensitive to the hot medium transport properties, especially to the parameter known as the heavy flavor diffusion coefficient 𝒟s [226].

Significant elliptic flow (v₂) was observed for charmed meson D⁰ in Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ by the STAR collaboration [96]. The D⁰s are fully reconstructed via the two body decay of charged pions and kaons with a branching ratio of (3.95 ± 0.03)%. The random combinatorial background of pions and kaons originating from primary vertices are significantly suppressed by precise measurements of the distance of closest approach (DCA) between tracks and primary vertex owing to the relatively large $c\tau$ of D⁰ mesons (≈123 μm). The precise measurements of DCA were provided by the heavy flavor tracker installed in STAR during 2014 and 2016. The v₂ results for D⁰ mesons in 10–40% central Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ are compared with those of light flavor hadrons (K_S, Λ, ${\text{Ξ}}^{-}$) [227]. At p_T<2 GeV/c, the v₂ for D⁰ mesons is found to be smaller than that of light flavor hadrons, exhibiting the mass-ordering behavior expected from hydrodynamics. At p_T>2 GeV/c, D⁰ v₂ is consistent with that of light flavor mesons such as K_S. The comparison of $v_2/n_q$ as a function of $\left(m_{\text{T}}-m_{\text{0}}\right)/n_q$, where nq is the number of constituent quarks in a hadron, shows that among these hadrons D⁰ elliptic flow follows the universal trend as in light hadrons. These comparisons indicate that charm quarks acquire significant flow through interaction with the strongly coupled QGP created in 10–40% Au+Au collisions at the top energy RHIC. Recent phenomenological models constrained by the D⁰ v₂ measurement as well as measurements of heavy flavor quarks v₂ using single electrons from heavy flavor hadron decays (HFE) suggest that the dimensionless charm quark spatial diffusion coefficient $2\pi T{\mathcal{D}}_s$ falls within the range of 2-5 in the vicinity of the critical temperature ([228] and references therein). This value is consistent with theoretical calculations from the quenched lattice QCD within large uncertainties. The dependence of the heavy flavor quarks diffusion coefficient Ds on heavy flavor quark momentum, as well as on temperature and baryon chemical potential of QGP is yet to be determined. The measurements of heavy flavor quarks collectivity in Au+Au collisions at energies below the RHIC top energy enabled by the RHIC BES program can shed new light on the temperature and baryon chemical potential dependence of the QGP transport parameter 𝒟s.

The elliptic flow of heavy flavor hadrons from RHIC BES program was measured in Au+Au collisions at $\sqrt{s_{\text{NN}}}=27\text{ GeV}$, 39 GeV, 54.4 GeV, and 62.4 GeV [229, 230]. The 39 GeV and 62.4 GeV data were collected in 2010 during the first phase of the RHIC BES program. The number of events used for the analyses were 87 and 38 million, respectively. The 27 GeV and 54.4 GeV data were collected in 2017 and 2018 between the first and second phases of the RHIC BES program. The number of events passing the event-level criterion were 240 and 570 million, respectively. With an order of magnitude difference of the number of events, combined with the significant energy dependence of heavy flavor hadron production cross section, the precision of the measurements at different energies varies considerably. The results from the 54.4 GeV collisions have the best precision. While the full reconstruction of heavy flavor hadrons in these data is not possible because of a lack of silicon vertex devices, the electrons from heavy flavor hadron decays are used as proxy for the heavy flavor hadrons. Electrons are identified using the inverse velocity calculated from the path length measured by the STAR TPC, with the TOF measured by the vertex position detector providing the start time measurement and the TOF detector providing stop time measurement. The electron candidates are further selected by the ionization energy loss in the TPC gas. The number of electrons are corrected for purity. Photonic electrons (PE) are the dominant source of background for heavy flavor decay electrons. These photonic electrons are produced via Dalitz decay of light mesons such as π⁰ and η and photon conversion in the detector material. The yield of non-photonic electrons (NPE) is calculated as follows: $N^{\text{NPE}}=N^{\text{INC}}-N^{\text{PE}},$ (4) where N^INC and N^PE represent the yield of inclusive and photonic electrons, respectively. The PE candidates are selected via the invariant mass distribution of inclusive electron candidates and partner electrons from the same event. The yield of photonic electrons can be expressed as: $N^{\text{PE}}=\left(N^{\text{UL}}-N^{\text{LS}}\right)/\epsilon ,$ (5) where N^UL and N^LS are the raw yield of unlike- and like-sign pairs, respectively, and ε is the partner electron acquiring efficiency. The v₂ of inclusive electron and photonic electrons is extracted using the event-plane η-sub method. The v₂ of NPE is calculated by: $N^{\text{NPE}}{v}_2^{\text{NPE}}=N^{\text{INC}}{v}_2^{\text{INC}}-N^{\text{PE}}{v}_2^{\text{PE}}-{\displaystyle \sum f_{\text{h}}\times N^{\text{INC}}}{v}_2^h,$ (6) where $v_2^{\text{INC}}$, $v_2^{\text{PE}}$, and $v_2^h$ are the v₂ of inclusive electrons, photonic electrons, and hadrons contaminated by inclusive electron candidates, respectively; fh is the hadron contamination fraction.

In addition to photonic electrons, other major background sources for heavy flavor decay electrons are from the decay of kaons ($K_{e3}$) and vector mesons (ρ, ω, and ϕ), which are eliminated using the following equation: $v_2^{\text{HFE}}=v_2^{\text{NPE}}\left(1+f_{{\text{K}}_{e3}}+f_{\text{VM}}\right)+f_{K_{e3}}\times v_2^{K_{e3}}+f_{\text{VM}}\times v_2^{\text{VM}}\mathrm{,}$ (7) where $f_{{\text{K}}_{e3}}$ and f_VM are the estimated yield ratios of electron decays of kaons and vector mesons, respectively, with respect to HFE. The residual non-flow contribution is estimated according to the HFE-hadron correlation in p+p collisions and hadron multiplicity in Au+Au collisions.

Figure 28 shows the elliptic flow coefficient v₂ of heavy flavor decay electrons as a function of p_T at midrapidity (|y|<0.8) in Au+Au collisions at $\sqrt{s_{\text{NN}}}=54.4\text{ GeV}$. The error bars and boxes depict statistical and systematic uncertainties, respectively. The hatched areas indicate the estimated non-flow contributions. Significant v₂ is observed at 0.5 GeV/c < p_T < 2 GeV/c; the average v₂ in the p_T range of 1.2 GeV/c to 2.0 GeV/c is 0.094±0.008(stat.)±0.014(syst.), whereas the estimated upper limit of non-flow is only 0.02. The red curve represents the projected charm quark decay v₂ assuming that open charm hadron v₂ follows NCQ scaling along with other light hadrons. Because charm quark is the dominant contributor of HFE in this pT range, the significant v₂ and the consistency between the data and the red curve indicate that charm quarks interact with the hot medium and may reach local thermal equilibrium in Au+Au collisions even though the center-of-mass energy is nearly a factor of four lower than the RHIC top energy.

Fig. 28Full-screen View

(Color online) v₂ as a function of p_T for heavy flavor decay electrons at midrapidity in Au+Au collisions at $\sqrt{s_{\text{NN}}}$ = 54.4 GeV [229], compared with TAMU [231] and PHSD [232, 233] model calculations

The two bands shown in Fig. 28 are calculated using two phenomenological models, TAMU [231] and PHSD [232, 233]. Both models assume that heavy flavor quarks interact with the QGP medium elastically. This assumption is generally accepted in the low p_T region. The elastic scattering is implemented in different ways in the two models. In the TAMU model, the microscopic elastic interaction between heavy flavor quarks and quarks/gluons in the hot, dense medium are evaluated using non-perturbative T-matrix calculations. The heavy flavor quark transport coefficient calculated is then fed into macroscopic Langevin simulation of heavy quark diffusion through the background medium and modeled by ideal 2+1D hydrodynamics. In the PHSD model, heavy flavor quarks interact with the off-shell massive partons in the QGP medium. The masses and widths of the partons in the QGP medium and the scattering probability are provided by the dynamical quasi-particle model. In both models, the heavy flavor quarks hadronize through both coalescence and fragmentation. In the PHSD model, the hadronized heavy flavor hadrons subsequently interact with other hadrons in the hadronic phase. Although the calculations from both TAMU and PHSD models are systematically lower than the measurements, the deviation is only 1-2 σ at p_T>0.5 GeV/c when considering the estimated upper limit of non-flow contribution. Furthermore, neither model considers the contribution from charm baryons, whose yield is measured to be evenly enhanced in heavy ion collisions relative to that of mesons [234]. This contribution will slightly increase HFE v₂ at p_T>1 GeV/c.

Heavy quarkonium is a bound state of heavy flavor quark and its antiquark. The pairs of heavy flavor quark and its antiquark are produced predominantly by the initial scattering in heavy ion collisions and are tightly bound together, which makes them less sensitive to interactions with other particles. However, it is believed that the color potential of the bound states is subject to modification when QGP is formed, resulting in the dissociation of heavy quarkonium [242-245]. The suppression of quarkonium yield in heavy ion collisions arising from a modification of the potential is considered as the ‘smoking-gun’ signature of deconfinement in QGP. This suppression is sensitive to the temperature profile of QGP because the modification of the potential between a heavy quark and its antiquark in QGP is sensitive to the temperature of the medium.

The suppression of J/ψ in heavy ion collisions was extensively studied in experiments at CERN SPS [246]. The production yield of quarkonium in heavy ion collisions was found to be affected by cold nuclear matter (CNM) effects. Suppression of the J/ψ yield beyond the expected CNM effects was observed in central Pb+Pb collisions at 17.3 GeV based on the results from proton and nucleus collisions and was considered as evidence of deconfinement in QGP [235].

However, the first quarkonium measurement in heavy ion collisions at RHIC was very puzzling. The J/ψ suppression, quantified by the nuclear modification factors, and its centrality dependence measured in Au+Au collisions at $\sqrt{s_{\text{NN}}}=200\text{ GeV}$ [247] was found to be consistent with that observed in Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=17.3\text{ GeV}$. However, J/ψ suppression is expected to be stronger at higher collision energies where the initial temperature is believed to be systematically higher. Furthermore, the suppression was observed to be significantly stronger at forward rapidity than at midrapidity, where the energy density is believed to be higher compared with that of forward rapidity. An additional production mechanism, (re)combination of charm quark and anti-charm quark in the QGP medium, was proposed to solve this problem [248, 249]. Unlike the static and dynamic screening of the potential between heavy quark and its antiquark, the (re)combination mechanism can enhance quarkonium yield in heavy ion collisions. Because the yield of quarkonium from the (re)combination mechanism is approximately proportional to the square of total heavy quark cross section, which exhibits a significantly increasing trend as the energy of collisions increases, the contribution from (re)combination mechanism should display a clear center-of-mass energy dependence. The RHIC BES program provides a unique opportunity for varying the initial temperature and number of charm quark pairs in an event, which should shed new light on the production mechanism of quarkonium in heavy ion collisions.

The data for the J/ψ production study were collected during the RHIC BES in 2010 by the STAR experiment at $\sqrt{s_{\text{NN}}}=39\text{ GeV}$ and 62.4 GeV [236]. The total number of minimum bias triggered events are 182 million and 94 million, respectively. The J/ψs are reconstructed considering their decays into electron-positron pairs. The electron daughters are identified by combining the information from the STAR TPC and TOF. The random combinatorial background is reconstructed using the mixed-event technique. The invariant mass spectrum for unlike-sign pairs from mixed-events are normalized to that for like-sign pairs and subtracted from that for unlike-sign pairs. The combinatorial background-subtracted invariant mass spectrum is fit to the J/ψ template using Monte Carlo simulations plus a linear function for the residual background to extract the J/ψ yields. The nuclear modification factors (R_AA) are calculated as follows: $R_{\text{AA}}=\frac{1}{T_{\text{AA}}}\frac{{\text{d}}^2{N}_{\text{AA}}/\text{d}{p}_T\text{d}y}{{\text{d}}^2{\sigma }_{\text{pp}}/\text{d}{p}_{\text{T}}\text{d}y}\mathrm{,}$ (8) where ${\text{d}}^2{N}_{\text{AA}}/\text{d}{p}_{\text{T}}\text{d}y$ is the efficiency and acceptance corrected J/ψ p_T spectrum measured in A+A collisions; T_AA is the nuclear overlap function from Glauber Monte Carlo simulations; and ${\text{d}}^2{\sigma }_{\text{pp}}/\text{d}{p}_{\text{T}}\text{d}y$ is the J/ψ production cross section in p+p collisions at the same energy as that of A+A collisions. The cross sections of J/ψ in p+p collisions at $\sqrt{s}=39\text{ GeV}$ and 62.4 GeV are derived by interpolating J/ψ data worldwide because there are no measurements available for p+p collisions at $\sqrt{s}=39\text{ GeV}$ and 62.4 GeV, and previous measurements near these two energies from the Intersecting Storage Ring collider experiments show discrepancies [250].

Figure 29 shows the center-of-mass energy dependence of J/ψ nuclear modification factors measured at midrapidity in central heavy ion collisions from SPS, RHIC, and LHC. The data from RHIC were measured in Au+Au collisions, whereas the data from SPS and LHC were measured in Pb+Pb collisions [235-239]. The error bars and boxes representing statistical and systematic uncertainties indicate that the J/ψ R_AA remains constant from $\sqrt{s_{\text{NN}}}=17.3\text{ GeV}$ to 200 GeV and substantially increases from RHIC top energy to LHC.

Fig. 29Full-screen View

(Color online) J/ψ nuclear modification factors (R_AA) as a function of center-of-mass energy in central heavy ion collisions. The solid circles represent the measurements at SPS, RHIC, and LHC [235-239], and the curves in the left (right) panel reflect the calculations from transport model I [240] (II [241])

The curves in the left and right panels of Fig. 29 depict the calculations from two transport models [240, 241]. The dot-dashed lines represent the contribution from primordial J/ψs which are affected by the static/dynamic color-screening of the potential in QGP medium and CNM effects, whereas the dashed lines represent the contribution from (re)combination. The solid lines represent the sum of the two components. Although both transport models can describe the data, except for the transport model I at SPS energy, the decomposed contributions from primordial and (re)combination are quite different, further indicating the limitations on the understanding of J/ψ production mechanisms in heavy ion collisions, which should be urgently addressed before extracting the properties of QGP. STAR has collected considerably more data samples at different energies during the second phase of the RHIC BES program. Preliminary results show that J/ψ suppression can be precisely measured in Au+Au collisions at 54.4 GeV, and the suppression measurements can be extended to energies down to 14.6 GeV, an energy below the SPS top energy. These new data will shed new light on the production mechanism of J/ψ in heavy ion collisions.

Di-lepton Production

Photons and dileptons (e⁺e^- or μ⁺μ^-) emerge at various stages throughout the space-time evolution of the nuclear medium formed in ultra-relativistic heavy ion collisions. As penetrating electromagnetic probes, dileptons remain unaffected by strong interactions, preserving undistorted information on their sources. These sources are manifested differently in various lepton-pair invariant mass (Mll) regions, which are typically categorized into three classes. In the low-mass region (LMR), below the ϕ mass (Mll< 1.1 GeV/c²), contributions predominantly arise from decays of light mesons (π⁰, η, ρ⁰, ω, ϕ). Investigation of ρ⁰ spectra modifications allows probing the in-medium hadronic properties, which are particularly sensitive to mechanisms of chiral symmetry restoration in QCD matter [251]. The expected modifications in dilepton yields within the LMR provide insights into the medium’s lifetime and transition from hadronic to partonic degrees of freedom [252]. In the intermediate-mass region (IMR), which lies between the ϕ and J/ψ masses ($M_{ll}\simeq 1.2-3\text{GeV}/c^2$), the invariant mass spectrum appears as a continuum arising from both heavy flavor decays and QGP thermal radiation. This provides an opportunity to directly measure the average temperature of the QGP [253] by extracting the inverse slope of the mass spectra, which remain unaffected by the blue shift of the expanding system. In the high-mass region (HMR), defined as Mll≥ 3 GeV/c², primary sources contributing to the dielectron spectrum are heavy flavor/quarkonium decays and the Drell-Yan process. Detailed discussions on the associated physics for HMR have been provided in the preceding section and are not reiterated here.

To achieve precise measurements of the aforementioned dileptons, detectors with large, uniform acceptance and excellent lepton identification capabilities are crucial. The integration of the TOF detector has paved the way for dilepton measurements at STAR. Specifically, by combining timing measurements from the TOF detector with momentum and ionization energy loss ($\langle \text{d}E/\text{d}x\rangle$) measurements from the TPC, robust identification of electrons over a wide p_T range is achieved. This identification is characterized by high efficiency and purity, facilitating comprehensive dielectron analysis. The identified electron and positron candidates are paired by opposite and same sign charges, called unlike-sign and like-sign pairs, respectively. The like-sign pairs are used to statistically model both the combinatorial and correlated backgrounds. Moreover, the subtraction of decay products from light mesons (known as the “cocktail”) is achieved through simulations. These dielectron spectra are of particular interest, as they are anticipated to carry radiation signatures from various stages of heavy ion collisions prior to freeze-out.

The acceptance-corrected excess dielectron mass spectra, following the careful removal of background contributions, have been thoroughly measured by the STAR collaboration across various collision energies [254, 255], as illustrated in Fig. 30. Accompanying these measurements are model calculations [253, 252, 256] depicting the total thermal radiation (solid lines), which consider contributions from both in-medium hadronic processes (dashed lines) and the QGP phase (dotted lines). Remarkably, the model predictions provide a coherent framework for interpreting the measured dielectron spectra across a wide energy range and invariant mass regions. In the low-mass region, the predominant hadronic radiation is primarily attributed to the in-medium ρ broadening, stemming from its interactions with the hadronic medium, particularly baryons. Notably, this model also yields a consistent description of the invariant mass spectrum of dimuon pairs measured by the NA60 experiment at the SPS [257]. The observed in-medium ρ broadening serves as a compelling indicator of the partial restoration of chiral symmetry within the hot QCD medium [258]. However, in the intermediate-mass region, contribution from QGP radiation is anticipated although current measurements still lack precision in this regime. Consequently, the search for and exploration of QGP thermal radiation remain pivotal future endeavors in dilepton experiments at both RHIC and LHC.

Fig. 30Full-screen View

(Color online) The dielectron invariant mass spectra in Au+Au collisions at $\sqrt{s_{\text{NN}}}=27\text{ GeV}$, 39 GeV, 62.4 GeV, and 200 GeV measured by the STAR collaboration [254, 255]. The spectra are shown after subtracting the hadronic background contributions (cocktail), and compared with theoretical model calculations [253, 252, 256]. The theoretical predictions represent the total thermal radiation (blue solid lines), including contributions from both in-medium hadronic processes (red dashed lines) and QGP thermal radiation (red dotted lines). The figure is sourced from [254]

To quantitatively compare the excess in the LMR, the integrated excess yield of dielectrons in the mass region 0.3 GeV/c² < Mll < 0.7 GeV/c² is normalized to the charged pion yield dN/dy to cancel out the volume effect. Fig. 31 shows the collision energy dependence of the integrated dielectron yield, as measured by HADES [259], NA60 [257], and STAR [260, 255, 254] collaborations. The figure also includes theoretical model calculations depicting the dielectron yields (dashed blue lines) and the fireball lifetime (solid red lines) [253]. Impressively, the model provides a commendable description of the energy dependence, illustrating a modest increase from the SPS to the top RHIC energy. This observed increase effectively tracks the fireball lifetime over a broad spectrum of collision energies. Notably, the STAR measurements presented here pertain to BES phase I. However, the subsequent analyses of BES-II data extend these measurements from 19.6 GeV down to 7.7 GeV, providing fresh insights into the properties of the hot medium within the high baryon density regime.

Fig. 31Full-screen View

(Color online) The collision energy dependence of the integrated dielectron yield, normalized by the charged pion yield dN/dy, in the mass region 0.3 GeV/c² < Mll < 0.7 GeV/c², as measured by HADES [259], NA60 [257], and STAR collaborations [260, 255, 254]. This comparative analysis is supplemented by model calculations of dielectron yields (dashed blue lines) and fireball lifetime (solid red lines) [253]. The figure is sourced from [5]

Traditionally, LMR dileptons have been used to explore the in-medium broadening of the ρ meson and its association with chiral symmetry restoration. The impact on the invariant mass distribution of dileptons is often overlooked, as it is considered a trivial thermal factor incorporated into models for data comparison [252, 253]. Recent observations by STAR [263] indicate that the broadening of the ρ meson are so extensive that the LMR dileptons can be used to determine the temperature of the thermal source responsible for LMR radiation. To extract this temperature, a fitting function that combines the in-medium resonance structure with the continuum thermal distribution is applied to the measured mass spectrum. In a vacuum, the mass line shape of the ρ decaying into dileptons is represented by a relativistic Breit-Wigner function, f^BW(M). Within a hot QCD medium, this line shape is modified, multiplied by the Boltzmann factor, $e^{-M/T}$, to account for phase space effects [263]. Furthermore, if ρ is completely dissolved in the medium, its mass spectral structure spreads out and approaches a smooth distribution similar to the dielectron continuum from QGP thermal radiation, described by $M^{3/2}{e}^{-M/T}$ [253]. Fig. 32 illustrates temperatures derived from BES-I dielectron data as a function of the baryon chemical potential μ_B. The chemical freeze-out temperature T_ch and μ_B are determined by applying statistical thermal models to the yields of hadron production. T_ch from several statistical thermal models and the QCD critical temperature T_C from lattice QCD [54, 32, 262] are shown in Fig. 32 as open circles and a shaded band, respectively. Similarly, temperatures extracted from previously published low-mass thermal dielectron spectra [261, 259] are presented. Notably, temperatures extracted from BES-I and SPS LMR closely align with T_ch from statistical thermal models and T_C from lattice QCD. This alignment suggests that dielectron emission at LMR is mainly influenced by ρ broadening during the phase transition (or mixed phase), and the chemical freeze-out temperature at RHIC BES energies lies at the phase transition boundary. Recent analyses by the STAR collaboration have extracted temperatures from IMR thermal dileptons for Au+Au collisions at $\sqrt{s_{\text{NN}}}=27\text{ GeV}$ and 54 GeV, as reported in [263]. The extracted temperatures are significantly higher than those from statistical thermal models and lattice QCD calculations, indicating that IMR dileptons primarily originate from the earlier partonic stage of the collisions. Further details can be found in [263].

Fig. 32Full-screen View

(Color online) Temperatures versus baryon chemical potential, with temperatures extracted from the in-medium ρ mass spectra of LMR. The earlier QGP stage region from NA60 [261] (diamonds) and LMR from HADES data [259] (inverted triangle) are also depicted. Chemical freeze-out temperatures extracted from statistical thermal models (SH, GCE, SCE) [54, 32] are represented as open circles. The QCD critical temperature T_C at finite μ_B, predicted by LQCD calculations [262], is shown as a yellow band. The figure is adapted from STAR [263]

In relativistic heavy ion collisions, dileptons emerge not only from hadronic processes but also through the interaction of the intense electromagnetic fields accompanying the colliding ions, known as the Breit-Wheeler process [264, 265]. These fields can be treated as a spectrum of equivalent photons, with the photon flux being proportional to the square of the particle’s charge (Z²), resulting in dilepton production scaling with Z⁴. Initially, dilepton production from the two-photon process was studied in ultra-peripheral collisions, where the impact parameter is large enough to avoid hadronic interactions. However, recent observations have shown that such photo-production also occurs in hadronic heavy ion collisions [266, 267], prompting theoretical advancements to describe these processes [268-270]. In events with hadronic overlap, dilepton photo-production occurs alongside hadronic interactions, offering a new method to probe the QGP, especially its electromagnetic properties. Data from peripheral collisions show discrepancies in the $p_{\text{T}}^2$ distribution between experimental results and theoretical models that do not consider the impact parameter dependence of photon kinematics, suggesting potential alternative origins of $p_{\text{T}}^2$ broadening, possibly linked to a postulated, trapped magnetic field or Coulomb scattering in the hot and dense medium [266, 267]. However, theoretical calculations that account for impact parameter dependence can explain the observed broadening [271-275], indicating the significant influence of the initial electromagnetic field strength, which were subsequently confirmed by the CMS measurements [276]. Future precision measurements at a toroidal LHC apparatus (ATLAS), STAR, and CMS will further explore these effects, potentially revealing medium induced modifications in dilepton kinematics.