A Study of the Properties of the QCD Phase Diagram in High-Energy Nuclear Collisions

With the aim of understanding the phase structure of nuclear matter created in high-energy nuclear collisions at finite baryon density, a beam energy scan program has been carried out at Relativistic Heavy Ion Collider (RHIC). In this mini-review, most recent experimental results on collectivity, criticality and heavy flavor productions will be discussed. The goal here is to establish the connection between current available data and future heavy-ion collision experiments in a high baryon density region.


Introduction
Most of the visible matter in our universe can be described by the Quantum Chromdynamics (QCD), the standard theory of strong interactions. In the beginning of the century, the new form of matter, the quark-gluon plasma (QGP) in which quarks and gluons are 'freed' in a much larger volume compared to that of nucleon's, was discovered in the largest heavy-ion colliders RHIC and LHC [1][2][3][4] at vanishing baryonic density. Soon after the discovery, a serious question was asked: what is the structure of the nuclear matter at high baryonic density?
Tremendous efforts from both experimental and theoretical sides have launched in order to address the question. Figure 1 summarizes the current status of the studies. At the zero baryonic density, the transition from QGP to hadronic matter is a smooth-crossover at T = 150 − 160 MeV [5][6][7][8][9], see dashed-line in the figure. These results are extracted from the state of the art Lattice gauge theory calculations. At high baryonic density, on the other hand, one would expect a first-order phase transition, shown as a black solid-line. Thermodynamically the first-order phase boundary line must end at finite baryonic density, this is the illusive QCD critical point (CP). Again, recent Lattice calculations have concluded that the QCD critical is 'unfavored' [10,11] when µ B /T < 2.5. The red-line is the chemical freeze-out curve extracted from the measured hadron yields. The collision energies with the corresponding accelerator complex are indicated at the top of the figure. Experimentally extracted chemical freeze-out parameters are shown as the red-line [12]. The dashed-line near µ B = 0 indicates the crossover transition from the Quark-Gluon-Plasma to Hadron Gas [5][6][7][8][9]. The solid-line and the square show the expected 1st-order phase boundary and its end point: critical point. The regions of µ B /T < 2 and < 3 are shown as red and blue dot-dashed lines, respectively. The Lattice QCD calculations have concluded that the QCD critical point is unlikely to exist at µ B /T < 2.5 [10,11]. The red-circle and the yellow-line represent the liquid-gas transition [13]. Regions for physics programs of RHIC beam energy scan and fixed-target (US) [14], NICA (Russia) [15], FAIR (Germany) [16] as well as the HIAF (China) [17] are indicated at the top of the plot.
Experimental status on the beam energy scan (BES) is highlighted in Figure 2. Plot (a) shows the chemical freeze-out temperature T ch as a function of the baryonic chemical potential µ B . Both ALICE at LHC and STAR at RHIC results clearly show that at the vanishing baryon density, i.e., at high collision energy, the data driven freeze-out temperature is consistent with the Lattice calculation, T ch ∼ 160 MeV. Many authors have tried to analyze the chemical freeze-out conditions. Those include the results fluctuations analysis from Lattice QCD [18][19][20][21][22], HRG model [12] and other methods [23,24]. In the low baryon density region, µ B ≤ 400 MeV, the µ B dependence of the freeze-out temperature is quite weak and the value of the freeze-out temperature is around 150-160 MeV. More dramatic drop of the temperature is seen in the high baryon density region. Plot (b) is the kaon over pion yield ratios, extracted from central heavy-ion collisions, as a function of the collision energies. While one observes the smooth increase of the negative kaon over pion ratio with the collision energy, the positive ratio shows a broad peak around √ s NN ∼ 8 GeV and eventually merged with the negative ratios at high collision energy, √ s NN ≥ 100 GeV, where the pair production becomes dominant. Due to the associate channel, N + N → N + Λ + K + , positive kaons carry information on baryon density. The peak in plot (b) implies the maximum freeze-out density reached at 8 GeV. Later in the discussions, we attribute the region of 2 ≤ √ s NN ≤ 8 GeV as the high baryon density region (HBDR) as indicated by the yellow-area in the plot. Collective flow (collectivity) and the critical behavior (criticality) are important aspects in high-energy nuclear collisions. In this short review, we will discuss the experimental status including the results on collectivity and criticality from high-energy nuclear collisions. RHIC has provided most recent data so we will focus on the information. At the end we will address the importance of the future fixed-target experiments such as STAR fixed-target program in BES-II, CBM at FAIR [16] as well the CEE at HIAF [17].

Beam Energy Dependence of the Collectivity
The main goal of high energy heavy-ion collisions, such as collisions at LHC and top energy collisions at RHIC, is to study the properties of new form of matter QGP. QGP is a thermalized (or nearly) system with partonic degree of freedom. The elliptic flow measurement of multi-strange hadrons and φ mesons indicates the partonic collectivity has been built up at the top energy heavy-ion collisions at RHIC [31-37]. The Heavy Flavor Tracker, a high resolution silicon detector system, which was installed in the year of 2013, provides high vertex position resolution. The significance of charmed hadron reconstruction is significantly improved. Thus the precise measurement of D 0 becomes possible. Figure 3 shows the v 2 results for D 0 , Ξ − , Λ and K 0 S [38]. A number of constituent quark (n q ) is tested by scaling both v 2 and m T − m 0 with n q . A simple quark coalescence or recombination model suggests the baryon v 2 would be 1.5 times of the meson v 2 assuming the collectivity has been attained in the partonic stage, as the number of constituent quarks for baryons is 3 where it is 2 for mesons. When discussing the n q scaling, we usually focus on the intermediate p T range in which the v 2 value saturates. The D 0 v 2 follows the n q scalings with selected multi-strange and strange hadrons in Figure 3. It indicates the collectivity of parton level has been built up from light flavor u, d quarks to strange and charm quarks. Since the mass of the charm-quark is much larger than the temperature reached in the system, the observed strong charm-quark collectivity can be interpreted as the thermalization of the medium created in the 200 GeV Au+Au collisions at RHIC [38]. To some extent, this result justified the phase diagram sketched in Figure 1. The main motivation of BES program is to explore the QCD phase boundary and critical point. The logic is straightforward: as the collision energy decreases, the conditions of QGP formation are no longer satisfied at some point. It offers us a unique experimental way to investigate the QCD phase structure. The transverse radial flow velocity β is obtained by fitting the transverse momentum p T spectra with a blast wave model [39]: where I 0 and K 1 are the modified Bessel functions and ρ(r) = tanh −1 β. The model assumes a radially boosted thermalized source with two key parameters, kinetic freeze-out temperature T kin and a transverse collective flow velocity β. Plot (a) of Figure 4 shows the extracted parameter β as a function of collision energy. The data points are taken from E802 [40][41][42][43], E866 [44,45], E877 [46], E895 [47], NA49 [48][49][50][51], STAR [12,52,53] and ALICE experiments [54] and references therein. The p T spectra of π ± , K ± , p andp are fitted simultaneously with the blast wave model. The p T range for simultaneous fitting are similar across all RHIC BES and LHC energies. A rapid increase of β is observed at low energies (<5 GeV), then a steady increase follows up to LHC energy. The six points from RHIC BES (7.7-39 GeV) are almost flat within uncertainties.   [44,45], E877 [46], E895 [47], NA49 [48][49][50][51], STAR [12,52,53] and ALICE experiments [54] and references therein. The data points of RHIC and LHC are from 0-5% central collisions. AGS and SPS energies are mostly from 0-5% and 0-7% central collisions respectively. (b) The p T -integrated v 2 in 20%-30% most central collisions (or similar centrality) from various collision energies. The data points are from E895 [55] for protons, NA49 [56] for pions and FOPI [57], E877 [58], CERES [59], STAR and ALICE [60] for charged hadrons. The RHIC results of inclusive charged particles for 130 and 200 GeV are from refs. [61][62][63][64][65]. The RHIC BES data are from refs. [66,67]. The energy regions of future collider and fix-target experiments, NICA and FAIR, are indicated in both plots.
The β parameter extracted from blast wave model reflects the transverse radial flow built-up in the collision system. The second order coefficient of final azimuth distribution in the momentum space, v 2 , is sensitive to the initial geometry and interactions of early stage of the collisions. It suggests that a non-monotonic variation could be observed around the so-called, "softest point of EOS" [68,69]. The "softest point of EOS" is usually defined as a strong drop of speed of sound (a minimum value) or a reduction in the pressure of the system during the dynamic evolution. Plot (b) of Figure 4 shows the p T -integrated v 2 from 20%-30% or similar centrality as a function of collision energy. The data points are from E895 [55] for protons, NA49 [56] for pions and FOPI [57], E877 [58], CERES [59], STAR and ALICE [60] for charged hadrons. The RHIC results of inclusive charged particles for 130 and 200 GeV are from refs. [61][62][63][64][65]. The RHIC BES data are from refs. [66,67]. The negative v 2 ( √ s NN <3 GeV) is known due to the "squeeze-out" effect [57]. An increasing trend is observed for p T integrated v 2 from AGS to LHC. It appears that the slope of v 2 with collision energy is steeper for 3-7.7 GeV compared to 7.7-2760 GeV, which is consistent with that we observe for β parameter. The v 2 of charged hadrons as a function of p T does not change significantly at RHIC BES and LHC energies. Due to the rise in mean p T which is expected from larger radial flow, the p T -integrated v 2 increases. It is consistent with the collision energy dependence of β parameter, as discussed in panel (a) of Figure 4. Non-monotonic behavior which is predicted by the softening of the equation of state for a system close to the critical temperature [68] is not observed. The first order coefficient of final azimuth distribution in the momentum space, v 1 (rapidity-odd), as a function of rapidity is sensitive to the system expansion during the early stage of collisions. Both hydrodynamic and nuclear transport models indicate that v 1 in the midrapidity region offers sensitivity to details of the expansion of the participant matter during the early collision stages [70][71][72]. Hydrodynamic plus first-order phase transition calculations suggest a minimum of net-baryon v 1 slope (dv 1 /dy) near mid-rapidity is a signal of phase transition between QGP and hadronic matter [73,74]. Net-particle is defined as the excess yield of a particle type over its anti-particle [75,76]. The v 1 of net-particle is defined as: v 1X = r(y)v 1X + [1 − r(y)]v 1net−X , where X represents particle,X represents the corresponding anti-particle, r(y) is the ratio of particle to anti-particle yield. Plot (a) of Figure 5 shows the v 1 slope relative to rapidity for net-proton, net-Λ and net-kaon. Similar energy dependence is observed for net-proton and net-Λ. The non-monotonic behavior is consistent with the hydrodynamical calculations with first-order phase transition [73,74]. Large divergence between dv 1 /dy of net-kaon and net-proton (net-Λ) is observed below √ s NN < 20 GeV, whereas all three agrees well at and above 20 GeV. More theoretical inputs are needed to understand the difference. At the same time the measurements of centrality dependence in the future BES program will further verify the energy dependence and constrain model calculations. In Ref. [76], the dv 1 /dy of φ mesons shows larger magnitude than pions and kaons at and above 14.5 GeV, and more interesting, the φ meson slope seems to increase sharply at 11.5 GeV. Because of the large statistical uncertainties, it is still not conclusive. It opens a new direction for both experimental and theoretical investigation on directed flow [77].  As discussed above, the v 2 of multi-strange hadrons and φ mesons are more sensitive to the parton level collectivity as their hadronic cross sections are smaller than light flavor hadrons [81-83]. The results from RHIC BES I suggest a possible drop of φ meson v 2 compared to other hadrons in Au+Au collisions at √ s NN = 11.5 and 7.7 with ∼2σ effect [78][79][80]. Data of high precision will be available with RHIC BES II. A significant difference in the v 2 values between particles and the corresponding anti-particles is observed at low energy heavy-ion collisions at RHIC. As shown in plot (b) of Figure 5, the difference is more pronounced for v 2 of baryons and anti-baryons when the collision energy is less than 20 GeV [78][79][80]. These differences naturally break the n q scaling discussed previously, as the number of constituent quarks are same for particle and the corresponding anti-particle. Several models try to explain the data [84-88]: the hydro + transport (UrQMD) calculation can reproduce the proton data, but not the meson data [84]; A analytic hydro model can quantitatively reproduce the π, K and proton data, but the flavor dependence (∆v is not consistent with data [85]; A Nambu-Jona-Lasino (NJL) model incorporating partonic and hadronic potentials can describe the data qualitatively, but not quantitatively [86,87]. New data from RHIC BES II, especially data of multi-strange hadrons, could offer more constrains on the model calculation.

Beam Energy Dependence of the Higher-Order Cumulants of Net-Particle Multiplicity Distributions and Light Nuclei Productions
Fluctuations of conserved quantities, such as net-baryon (B), net-charge (Q) and net-strangeness (S), are sensitive observables to search for the QCD critical point in heavy-ion collisions [89-93]. The higher-order cumulants (C n , the n th order cumulants), which can be used to quantify the fluctuations and describe the shape of the event-by-event multiplicity distributions, are predicted to be sensitive to the correlation length (ξ) of the system as C 4 ∝ ξ 7 and C 3 ∝ ξ 4.5 [90,94]. The various order cumulants and cumulant ratios can be expressed in terms of moments as C 2 = σ 2 , C 3 = Sσ 3 , C 4 = κσ 4 and where σ 2 , S and κ are variance, skewness and kurtosis, respectively [95]. The various order cumulants are extensive quantities and are proportional to the system volume, which is difficult to be measured in heavy-ion collisions. By taking the ratio between various order cumulants, the system volume can be cancelled to the first order and are directly related to the ratios of the thermodynamic susceptibilities ( . The red and blue regions in the density plot denote the negative and positive contributions to the fourth order cumulant, respectively. Experimentally, by tuning the beam energy, the T and µ B at chemical freeze-out are varied accordingly. The green dashed line represents the chemical freeze-out points (T, µ B ) passing through the critical region when one varies the beam energies. Figure 6 (right) shows fourth order fluctuation κσ 2 as a function of baryon chemical potential (µ B ). Due to the negative and positive critical contributions near the critical point, the κσ 2 will show a non-monotonic energy or µ B dependence with respect to the non-critical baseline. This is the characteristic experimental signature of the critical point we are looking for in the heavy-ion collision experiment. Theoretically, the properties of QCD phase diagram at finite baryon density and the signatures of conserved charge fluctuations near the QCD critical point have been extensively studied by various model calculations, such as Lattice QCD [10,[18][19][20][21][22]98], NJL, PNJL model [99][100][101][102][103][104][105][106], PQM, FRG model [107][108][109], Dyson-Schwinger Equation (DSE) method [110][111][112][113], chiral hydrodynamics [114] and other effective models [94, [115][116][117][118][119]. However, one should keep in mind that the above results are under the assumption of thermal equilibrium with infinite and static medium. In the real heavy-ion collisions, there exists the effects of finite size/time [120][121][122], non-equilibrium [123][124][125][126][127] and thermal blurring effects [128]. Dynamical modeling of heavy-ion collisions by implementing both the critical and those background effects are ongoing [129][130][131]. Before turning to the experimental status, we would like to stress that there is a long history of using the higher-order cumulants to extract the information on that state created in high-energy nuclear collisions. As extensively discussed in Ref.
[26], hadron yields or their ratios, which are the first order moment of the multiplicity distributions, have been used for determine the nature of thermalization in such collisions. Once the thermalization is established [132], on the other hand, the higher-orders cumulants of the very same distributions can be used to study the fine structures of the QCD matter. For example, the critical point [94], the nature of the crossover transition [5,107] and the phase boundary [133] at vanishing and large net-baryon region, respectively.
To make precise measurements, various corrections and techniques have been applied in the data analysis, those include : (1) Select proper collision centralities to avoid auto-correlations and suppress volume fluctuations [143,144], (2) Detector Efficiency Correction [145][146][147][148][149], (3) Centrality Bin Width Correction (CBWC) [143], (4) Statistical error estimation with Delta theorem and/or Bootstrap [143,150]. Figure 7 shows the measured event-by-event net-charge, net-kaon and net-proton multiplicity distributions of three different centralities in Au+Au collisions at √ s NN = 14.5 GeV. Those are raw distributions and not corrected for detector efficiency and acceptance. One should apply the efficiency correction and CBWC to obtain the final efficiency corrected cumulants. In Ref. [150], we have shown that the statistical uncertainties of the cumulants (C n ) strongly depend on the width of the distributions (error ∝ σ n / √ N). In general, the widths of the distributions in central collisions are wider and with larger mean values than those from peripheral collisions. Further, the widths of the net-charge distributions are much wider than those of net-proton and net-kaon in the same centrality. That's the reason why we observe larger statistical uncertainties in central collisions than those from peripheral. Assuming the measured particles are emitted from many independent sources in the fire ball, the multiplicity distributions in central collisions would be more symmetric and close to Gaussian distribution than the peripheral based on the central limit theorem (CLT). If those sources are identical and uncorrelated, the cumulant ratios are expected to be a constant as a function of collision centralities.   which are plotted as the dashed lines. We found that the σ 2 /M of net-charge, net-kaon and net-proton monotonically increase when increasing the collision energy. The Sσ/Skellam and κσ 2 show weak energy dependence for net-charge and net-kaon measurements. We didn't observe significant deviations of net-charge and net-kaon cumulant ratios Sσ/Skellam and κσ 2 from the Poisson expectations and UrQMD calculations within uncertainties. However, a clear non-monotonic energy dependence of net-proton κσ 2 was observed in top 0-5% central Au+Au collisions. The 0-5% net-proton κσ 2 values are close to unity for energies above 39 GeV and show large deviations below unity around 19.6 and 27 GeV, and then increasing above unity below 19.6 GeV. The UrQMD calculations of net-proton κσ 2 displaying a strong suppression below unity at lower energies is due to the effects of baryon number conservation [154][155][156][157]. However, this suppression is not observed at low energies in the STAR data. Another transport model (JAM model) study further demonstrates that the resonance weak decay and hadronic re-scattering have very small effects on the proton number fluctuations (C 1 -C 4 ) at low energies [158].  In Figure 9, we summarize the energy dependence of κσ 2 of net-charge, net-kaon and net-proton multiplicity distributions in Au+Au collisions measured by the STAR experiment. For comparison, the net-charge results in Au+Au collisions at √ s NN = 7.7, 19.6, 27, 39, 62.4 and 200 GeV measured by the PHENIX experiment [159] are shown in the panel (b). We found that the κσ 2 of the net-charge and net-kaon multiplicity distributions measured by the STAR experiment show larger statistical uncertainties than those of net-proton κσ 2 . This can be understood as the statistical uncertainties of κσ 2 depend on the width (σ) of the multiplicity distributions and the particle detecting efficiencies ( ) in the detector as error(C n /C 2 ) ∝ σ n−2 /( √ N n ) [146]. The width of the net-charge distributions are larger than those of net-proton and net-kaon. Meanwhile, due to decays, the efficiency of kaon (∼40%) is much lower than proton (∼80%). It is the reason why we observe larger statistical uncertainties for net-kaon fluctuations than net-proton. For the net-charge and net-kaon κσ 2 from STAR, we observe weak energy dependence within current statistical uncertainties. The PHENIX net-charge κσ 2 are with much smaller statistical uncertainties than the results from STAR. This is due to smaller acceptance of PHENIX detector than the STAR detector, thus the width of the net-charge multiplicity distributions measured by the PHENIX experiment is much narrower than those measured by STAR. We observe a clear non-monotonic energy dependence for net-proton κσ 2 in the most central (0-5%) Au+Au collisions with a minimum around 19.6 GeV. This non-monotonic behavior cannot be described by various model calculations without the physics of phase transition and critical point.  Figure 10 shows the energy dependence of the fourth-order fluctuations (κσ 2 ) of net-proton from the top 5% central Au+Au collisions measured by STAR experiment [139]. Recent result from the HADES experiment is also shown in the figure. Note that there are differences in the data shown in Figure 10 : while STAR data points are from the top 5% central collisions and |y| < 0.5, 0.4 < p T < 2.0 GeV/c, the HADES data is from the top 10% central Au+Au collisions and |y| < 0.4, 0.4 < p T < 1.6 GeV/c. From 200 GeV to 7.7 GeV, non-monotonic energy dependence is clearly shown in the κσ 2 of net-proton multiplicity distributions and one can observe a strong enhancement at the highest µ B ∼ 420 MeV, corresponding to the Au+Au central collisions at √ s NN = 7.7 GeV. This might indicate attractive correlations between nucleons in nature at the large baryon density region. However, interestingly, the strong enhancement in the fourth-order fluctuation seems disappeared as shown by the HADES result at √ s NN = 2.4 GeV [160].  In Figure 10, the results from the transport model UrQMD (grey band) show a monotonic decrease from low to high baryon density region, which is due to the effect of baryon number conservation in high-energy nuclear collisions. Note that in the Poisson limit, the absence of criticality or other dynamical correlations, the κσ 2 is expected to be unity. Besides the conserved charge fluctuations, the light nuclei production is predicted to be sensitive to the baryon density fluctuations assuming that the light nuclei is formed from the nucleon coalescence. Model calculations show that the yield ratio between deuteron, triton and proton, N t × N p /N 2 d is related to the neutron density fluctuations, thus can be used to search for the QCD critical point in heavy-ion collisions [164,165]. Experimentally, the STAR experiment has measured the production of deuteron As shown in Figure 11, non-monotonic energy dependence is observed for the yield ratio, N t × N p /N 2 d , in 0-10% central Au+Au collisions with a peak around 20-30 GeV [166][167][168]. The yield ratios measured by STAR experiment below 20 GeV are consistent with the results calculated from NA49 experiment [164]. Since there is no critical physics implemented in the JAM model, the results of central (b < 3 fm) Au+Au collisions from JAM model is also plotted as blue band in Figure 11 for comparison [166]. The model results show a flat energy dependence and cannot describe the non-monotonic trend observed in the STAR data. The current STAR results shown in Figure 11 is for 0-10% centrality, it is also worthwhile to perform centrality dependence study on this yield ratio. On the other hand, more theoretical studies and dynamical modeling of heavy-ion collisions with critical physics are needed to understand whether this non-monotonic behavior is related to the QCD critical fluctuations.

Beam Energy Dependence of the Heavy-Flavor Production
Since the masses of heavy flavor quarks are much larger than the temperature of the system created in the high-energy nuclear collisions, they can be used as clean probes of the medium properties at early stage of the collisions. As shown in Figure 12, heavy flavor quark masses are all generated in the electro-weak sector while light quarks (u, d, and s) are dominated by the spontaneous breaking of chiral symmetry in QCD. Thus heavy quarks keep massive when participating in strong interactions. Due to their large masses, these heavy flavor quarks are primarily pair-created in initial hard pQCD processes. These facts making heavy quark hadrons are ideal for studying the medium effects including the thermalization of the system. One example has already discussed in previous section, see Figure 3.  [169]. The plot is taken from Ref. [170].
The QCD calculations can evaluate the charm production cross sections at high energies via a perturbation scheme in p+p collisions [171,172]. Figure 13 shows the charm production cross sections at midrapidity as a function of p T in p+p collisions at √ s = 7 TeV, 1.96 TeV, 500 GeV and 200 GeV from ALICE [173], CDF [174] and STAR [175,176] experiments, respectively. Within uncertainties the Fixed-Order-Next-to-Leading-Logarithm (FONLL) calculations [172] agree with data. The data points are more on top of the upper limit of the theoretical uncertainties for all of the collision energies. In heavy-ion collisions, charm quark interacts with the QGP matter when traversing in the medium. The transverse momentum of charm quark is modified by the medium via energy loss or collective flow. However, the total number of charm quarks may keep conserved since they are produced in initial hard processes before the QGP formation and there is no more charm quark created later via thermal production at RHIC energies. Figure 14 shows the p T -integrated cross section for D 0 production per nucleon-nucleon collision dσ NN /dy| y=0 from different centrality bins in √ s NN = 200 GeV Au+Au collisions for the full p T range (a) and for p T > 4 GeV/c (b), respectively [177]. The result from the p+p measurement at the same collision energy is also shown in both panels [175].  uncertainties. This is consistent with charm quarks lose more energy in more central collisions at high p T . However, for the dσ NN /dy| y=0 integrated over full p T range shows approximately a flat distribution as a function of N part . The values for the full p T range in mid-central to central Au+Au collisions are smaller than that in p+p collisions with ∼1.5σ effect considering the large uncertainties from the p+p measurements. The total charm quark yield in heavy-ion collisions is expected to follow the number-of-binary-collision scaling since charm quarks are conserved at RHIC energies. However, the cold nuclear matter (CNM) effect including shadowing could also play an important role. In addition, hadronization through coalescence could alter the hadrochemistry distributions of charm quark in various charmed-hadron states which may lead to the reduction in the observed D 0 yields in Au+Au collisions [178]. For instance, hadronization through coalescence can lead to an enhancement of the charmed baryon Λ + c over D 0 yield [179][180][181], and together with the strangeness enhancement in the hot QCD medium and sequential hadronization, can also lead to an enhancement in the charmed strange meson D + s yield relative to D 0 [180][181][182]. The STAR Heavy Flavor Tracker (HFT) with a silicon pixel detector achieved ∼30 µm spacial resolution of the track impact parameter to the primary vertex allows a topological reconstruction of the decay vertices of open charm hadrons. Figure 15 left panels show the charmed baryon over meson ratio compared with light and strange baryon over meson ratios [183,184] (a) and various models (b). The Λ c /D 0 ratio is comparable in magnitude to the Λ/K 0 s and p/π ratios and shows a similar p T dependence in the measured region. A significant enhancement is seen compared to the calculations from the latest PYTHIA 8.24 release (Monash tune [185]) without (green solid curve) and without (magenta dot-dashed curve) color reconnections (CR) [186]. The implementation with CR is found to enhance the baryon production with respect to mesons. However, both calculations fail to fully describe the data and its p T dependence. Figure 15b also shows the comparison to various models with coalescence hadronization of charm quarks [179][180][181][182]. The comparisons suggest coalescence hadronization plays an important role in charm-quark hadronization in the presence of QGP. Also, the data can be used to constrain the coalescence model calculations and their model parameters.  Figure 15 right panel shows the D s /D 0 ratio as a function of p T compared to coalescence model calculations for 0-10% (c) and 10%-40% (d) collision centralities. Several models incorporating coalescence hadronization of charm quarks and strangeness enhancement are used to describe the p T dependence of D s /D 0 ratio. Those models assume that D ± s mesons are formed by recombination of charm quarks with equilibrated strange quarks in the QGP [179][180][181][182]. In particular, the sequential coalescence model together with charm quark conservation [180] considers that more charm quarks are hadronized to D ± s mesons than D 0 since the former is created earlier in the QGP, which results in further enhancement of D s /D 0 ratio in Au+Au collisions relative to p + p collisions.
D meson R AA and v 2 have been observed similar as light flavor hadrons in 200 GeV Au+Au collisions [177], which indicates charm is thermalized in the system with T ∼170 MeV. In low energy region, such as the energies in RHIC beam energy scan program, it is of particular interest to measure open charm hadron production in a relative smaller and colder system compared to top energy at 200 GeV. This may provide a chance to tell us in what temperature charm behaves different from light flavors. However, in low energy region, the perturbation algorithm in theoretical calculations of charm production cross section becomes invalid, which may result in large theoretical uncertainties. Meanwhile the charm production cross section drops rapidly when collision energy decreases, it is very challenging to measure open charm production at low energies. The previous measurements at SPS energies are with large uncertainties [188,189]. Since the HFT detector was taken out from STAR for the BES-II runs together with low cross section, it is impossible to reconstruct open charm hadrons via hadronic decay channels, the electron production from heavy flavor semi-leptonic decays becomes the unique way to measure heavy flavor productions at low energy. Figure 16 shows the v 2 of electrons from heavy flavor decays as a function of p T in √ s NN = 54 and 200 GeV Au+Au collisions [190] as solid circles and open stars, respectively. A semi-empirical exponential function [191] is used to fit all the data points and the ratios of data over the fit function are shown in the bottom panel. The result in 54 GeV agrees with that in 200 GeV within uncertainties, which may suggest charm quarks are still thermalized in 54 GeV Au+Au collisions. On the other hand, it is of interest to repeat the same measurement in lower energies, such as 27 GeV from STAR BESII experiment.  [190]. Red curve denotes an empirical reversed exponential fit [191] to all the data points. Bottom panel: The ratios of data over the fit function. Vertical bars and brackets denote statistical and systematic uncertainties, respectively.

Future Upgrades and Physics Program at High Baryon Density Region
The RHIC BES program II and future FAIR and NICA experiments will focus on collision energy below 20 GeV offering us a unique opportunity to explore the QCD phase structure at high baryon density region. In Figure 18, interaction rates from both collider experiments and fixed-target experiments are shown. The region for the high baryon density, largely covered by the fixed-target experiments, is highlighted with yellow. In the following, we discuss few key measurements with the future experimental facilities. Again the discussions are arranged around the headlines of Collectivity, Criticality and Heavy Flavor Productions.
Collectivity: The flow results from top energy heavy-ion collisions at RHIC indicate that the partonic collectivity has been built up from light u, d and s quarks to heavy c quark as well. This is one of the most important experimental evidences for the creation of the QGP in high-energy nuclear collisions [1][2][3][4]. As a function of the collision energy, both radial and elliptic flow show an increasing trend, say above √ s NN ≈ 15 GeV (Figure 4). Above that energy, the v 1 slope of net-particles for both baryons and mesons, is observed to be almost the same, the v 2 difference between particle and anti-particle also becomes similar ( Figure 5). While the dv 1 /dy shows large divergence between net-kaon and net-proton (and net-Λ), the particle and anti-particle v 2 difference splits between baryons and mesons dramatically below √ s NN = 15 GeV, see Figure 5. All of these observations imply that the medium properties created in heavy-ion collisions would be different above/below √ s NN = 15 GeV.
In collectivity, two noticeable observations are the splitting between baryon's and meson's v 1 and v 2 in the low energy. Mesons such as kaons and φ−mesons are important, especially the φ−meson as it has the similar mass of proton. The precise results of φ−meson's v 1 and v 2 will reveal the origin of collectivity at the high baryon density region. In addition, the ratio of N(K − )/N(φ) will shed light on the production mechanism. It could be treated as a micro-laboratory for understanding the quarkonia productions in nucleus-nucleus collisions.
Criticality: One of the main goal of RHIC Beam Energy Scan program is to search for the QCD critical point, which is the end point of the first order phase boundary in the QCD phase diagram. The experimental confirmation of the existence of the CP will be a landmark of exploring the QCD phase structure. Near the QCD critical point, the density fluctuations and correlation length will diverge. The conserved charge fluctuations and light nuclei productions have been proposed as sensitive observables to search for the signature of QCD phase transition and the QCD critical point. Experimentally, the STAR experiment has measured the higher order cumulants of net-particle distributions and light nuclei productions (deuteron and triton) in Au+Au collisions at √ s NN = 7.7 to 200 GeV. In Figures 10 and 11, it has been observed that both fourth order fluctuations of net-proton (κσ 2 ) and light nuclei yield ratio N t × N p /N 2 d show non-monotonic energy dependence in central Au+Au collisions, with a minimum and peak around 20 GeV, respectively. Although the two measurements are of different order of fluctuations, the two observations are consistent with the expectation of model calculations with CP physics and might suggest that the created system skims close by the CP receiving the contributions from critical fluctuations. To confirm the above two observed non-monotonic energy dependence trends in BES-I, the second phase of Beam Energy Scan (BES-II) has been planned at RHIC (2019-2021). It will allow us to have 10-20 times more statistics at energies √ s NN = 7.7-19.6 GeV. In addition, one observes large changes between 19.6 and 14.5 GeV in the energy dependence of net-proton kurtosis ( Figure 10) and light nuclei yield ratio N t × N p /N 2 d (Figure 11) measured in the RHIC BES-I data. This could indicate that the QCD critical point is put by nature between the thermodynamic condition (T, µ B ) of 19.6 and 14.5 GeV. Thus, it is important to conduct a finer beam energy scan between these two energies, i.e., 19.  GeV) as well as the fixed-target projects including HADES (diamonds), HIAF (triangles) and FAIR (filled triangles). STAR fixed-target range is indicated with filled red circles.
On the other hand, it is predicted that the higher order conserved charge fluctuations, such as sixth order (C 6 ) or eighth order (C 8 ) cumulants, should be more sensitive to the phase transition. If the chemical freeze-out temperature in heavy-ion collisions are close enough to the phase boundary, the sixth and eighth order fluctuations could show negative values [98,107]. STAR experiment has measured the centrality dependence of sixth order (C 6 [216,217]. The negative sign of net-proton C 6 /C 2 observed at 200 GeV could be an experimental evidence of smooth crossover at small baryon chemical potential [98,107]. In future fixed target experiments, with much more statistics of low energy data, we can perform precise measurements of those higher order cumulants of conserved charges at the high baryon density region. Heavy Flavor Production: FAIR-CBM and NICA-MPD experiments with advanced fast detector technology under high luminosity beam condition will provide unique chance to measure open and hidden charm hadrons with large statistics close to the production energy threshold [162,218]. It is expected that this measurements will improve the precision of the total charm cross section at low energies and will provide constraints on pQCD calculations, as well as the unknown interactions between charmed particles and cold hadronic medium. Taking the prediction of the HSD model [219], the yield obtained in one week of running of CBM detectors with 10 MHz event rate would be about 300 J/ψ for central Au+Au collisions at 10A GeV, and about 600 J/ψ for central Ni+Ni collisions at 15A GeV. In the latter case, also open charm production can be studied at a rate of 300 kHz with a silicon vertex detector MVD in operation for charmed hadron decay vertex reconstruction. As a result, the expected yield in central Ni+Ni collisions at 15A GeV will be about 30 D mesons per week. This would be sufficient for cross section measurement and an analysis of charmonium propagation and absorption in dense baryonic matter based on the ratio of hidden to open charm at low energy. In order to extend the coverage to even larger baryon density region, STAR has developed a fixed-target (FXT) program. As shown in Figure 18, a gold-target (1% interaction length) is placed at the right entrance of TPC. The end cap time-of-flight wall will be constructed at approximately the other side of the TPC entrance. The time-of-flight detectors are on loan from the CBM experiment at FAIR [220,221]. In addition, the inner-TPC upgrade [222] will extend the rapidity coverage, essential for the search for the QCD critical point measurement. STAR will be setup in such a way that data taking from both colliding and FXT modes will take place concurrently. With this configuration, STAR detector system will measure particle productions and correlations in Au+Au collisions from √ s NN = 3-19.6 GeV, extending its coverage of baryon chemical potential from about µ B = 400 MeV to µ B ∼ 750 MeV. The center of mass energy from the highest energy of the FXT mode is overlay with the lowest colliding mode at √ s NN = 7.7 GeV and the lower part of the FXT energies overlap with the future collision energies provided by CBM at FAIR [16]. These allow systematic crosschecks on many of the observables in STAR experiment, for both colliding and FXT modes, and CBM experiment for the FXT mode. The BES program II and future FAIR and NICA experiments will focus on the high baryon density region (<20 GeV), offer us a unique opportunity to explore the QCD phase structure. In summary, the precise flow measurements of φ mesons and multi-strange hadrons with STAR BES-II and future fixed-target experiments will reveal the degree of freedom originates from partonic or hadronic level at the high baryon density region. The confirmation of non-monotonic energy dependence in central Au+Au collisions for net-proton κσ 2 and/or light nuclei yield ratio N t × N p /N 2 d with STAR BES-II and future fixed-target experiments will provide crucial experimental evidences for establishing the case for the discovery of the QCD critical point. The energy dependence of heavy flavor measurements will provide crucial information on thermalization of the system and provide unique opportunity to study the unknown interactions between heavy quark and the cold nuclear matter. A great deal of new information on the QCD phase diagram will be extracted with current and planned heavy-ion collision programs.
Author Contributions: X.L., S.S., N.X. and Y.Z. contribute equally to this paper. All authors have read and agreed to the published version of the manuscript.