Nuclear Science and Techniques

NUCLEAR PHYSICS AND INTERDISCIPLINARY RESEARCH

Further developments of a multi-phase transport model for relativistic nuclear collisions

Zi-Wei Lin，

Liang Zheng

Vol.32, No.10

Article number 113

Published in print 01 Oct 2021

Available online 11 Oct 2021

DOI：10.1007/s41365-021-00944-5

8900

A multi-phase transport (AMPT) model was constructed as a self-contained kinetic theory-based description of relativistic nuclear collisions as it contains four main components: the fluctuating initial condition, a parton cascade, hadronization, and a hadron cascade. Here we review the main developments after the first public release of the AMPT source code in 2004 and the corresponding publication that described the physics details of the model at that time. We also discuss possible directions for future developments of the AMPT model to better study the properties of the dense matter created in relativistic collisions of small or large systems.

QGPTransport modelHeavy-ion collisions

1　Introduction

In high energy heavy ion collisions [1], a hot and dense matter made of parton degrees of freedom, the quark-gluon plasma (QGP), has been expected to be created [2]. Experimental data from the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [3-8] strongly indicate that the QGP is indeed created in heavy ion collisions at high energies [9]. Comprehensive comparisons beween the experimental data and theoretical models are essential for the extraction of key properties of the high density matter, including the structure of the QCD phase diagram at high temperature and/or high net-baryon density. Many theoretical models including transport models [10-14], hydrodynamic models [15-18], and hybrid models [19-21] have been constructed to simulate and study the phase space evolution of the QGP.

A multi-phase transport (AMPT) model [13] is one such model. The AMPT model aims to apply the kinetic theory approach to describe the evolution of heavy-ion collisions as it contains four main components: the fluctuating initial condition, partonic interactions, hadronization, and hadronic interactions. The default version of the AMPT model [11, 22] was first constructed. Its initial condition is based on the Heavy Ion Jet INteraction Generator (HIJING) two-component model [23, 24], then minijet partons enter the parton cascade and eventually recombine with their parent strings to hadronize via the Lund string fragmentation [25]. The default AMPT model can well describe the rapidity distributions and transverse momentum (p_T) spectra of identified particles observed in heavy ion collisions at SPS and RHIC. However, it significantly underestimates the elliptic flow (v₂) at RHIC.

Since the matter created in the early stage of high energy heavy ion collisions is expected to have a very high energy density and thus should be in parton degrees of freedom, the string melting version of the AMPT (AMPT-SM) model [26] was then constructed, where all the excited strings from a heavy ion collision are converted into partons and a spatial quark coalescence model is invented to describe the hadronization process. String melting increases the parton density and produces an over-populated partonic matter [27], while quark coalescence further enhances the elliptic flow of hadrons [26, 28]. As a result, the string melting AMPT model is able to describe the large elliptic flow in Au+Au collisions at RHIC energies with a rather small parton cross section [26, 29].

The source code of the AMPT model was first publicly released online around April 2004, and a subsequent publication [13] provided detailed descriptions of the model such as the included physics processes and modeling assumptions. The AMPT model has since been widely used to simulate the evolution of the dense matter created in high energy nuclear collisions. In particular, the string melting version of the AMPT model [13, 26] can well describe the anisotropic flows and particle correlations in collisions of small or large systems at both RHIC and LHC energies [13, 26, 30-33]. The AMPT model is also a useful test bed of different ideas. For example, the connection between the triangular flow and initial geometrical fluctuations was discovered with the help of AMPT simulations [34], and the model has also been applied to studies of vorticity and polarization in heavy ion collisions [35-37].

Experimental data from heavy ion collisions fit with hydrodynamics-inspired models suggest that particles are locally thermalized and possess a common radial flow velocity [38]. Large momentum anisotropies such as the elliptic flow [39] have been measured in large collision systems, as large as the hydrodynamics predictions [7, 40]. This suggests that the collision system is strongly interacting and close to local thermal equilibrium [9]. Transport models can also generate large anisotropic flows. The string melting AMPT model [13, 26] can describe the large anisotropic flows with a rather small parton cross section of ˜3 mb [26] and the flow enhancement from quark coalescence [26, 28, 29, 41, 42]. Without the quark coalescence, a pure parton transport for minijet gluons requires an unusually large parton cross section of ˜40- -50 mb [29, 43] for the freezeout gluons to have a similar magnitude of elliptic flow as charged hadrons in the experiments. This minijet gluon system, despite a factor of ˜2.5 lower parton multiplicity at mid-rapidity, has a factor of ˜6 smaller mean free path than the string melting AMPT model for 200A GeV Au+Au collisions at impact parameter b=8 fm [29]. In general, for large systems at high energies transport models tend to approach hydrodynamics since the average number of collisions per particle is large and thus the bulk matter is close to local equilibrium. Hydrodynamics models and transport models are also complementary to each other. For example, hydrodynamics models provide a direct access to the equation of state and transport coefficients, while transport models can address non-equilibrium dynamics and provide a microscopic picture of the interactions.

Recent data from small systems, however, hint at significant anisotropic flows in high multiplicity pp and pPb collisions at the LHC [44] and p/d/³He+Au collisions at RHIC [45, 46]. Hydrodynamic calculations seem to describe the experimental data well [47, 48]. The AMPT-SM model also seems to describe the measured correlations [30]. This suggests that the collision of these small systems might create a QGP as well, in contrast to naive expectations. On the other hand, it is natural to expect hydrodynamic models and transport models to be different for small colliding systems due to non-equilibrium dynamics. Indeed, recently it has been realized that parton transport can convert the initial spatial anisotropies into significant anisotropic flows in the momentum space through the parton escape mechanism [49, 50], especially in small systems where the average number of collisions per particle is small. Kinetic theory studies also show that a single scattering is very efficient in changing the particle momentum distribution [51]. There are also many studies on whether and how hydrodynamics could be applicable to small systems [52, 53]. In addition, there are active debates on whether the momentum anisotropies in small systems mainly come from initial state correlations [54, 55] or final state interactions [49-51, 56, 57]. Furthermore, the differences between the anisotropic flow data of small systems from different collaborations still need to be fully resolved [46, 58, 59]. Therefore, the system size dependence of various observables, particularly the anisotropic flows from small to large systems, could provide key information on the origin of collectivity.

The paper is organized as follows. After the introduction, in Sec. 2 we review the main developments of the AMPT model after the first public release of its source code in 2004 [13, 60, 61]. They include the addition of deuteron productions in Sec. 2.1, the string melting model that can simultaneously reproduce the yield, transverse momentum spectra and elliptic flow of the bulk matter in heavy ion collisions in Sec. 2.2, the new quark coalescence model in Sec. 2.3, incorporation of the finite nuclear thickness along beam directions in Sec. 2.4, incorporation of modern parton distribution functions of nuclei in Sec. 2.5, improved treatment of heavy quark productions in Sec. 2.6, the introduction of local nuclear scaling of key input parameters to describe the system size dependence in Sec. 2.7, incorporation of PYTHIA8 and nucleon substructure in the initial condition in Sec. 2.8, and benchmark and improvement of the parton cascade algorithm in Sec. 2.9. We then briefly review other developments of the AMPT model in Sec. 3. Finally, in Sec. 4 we summarize and discuss possible directions for further developments of the AMPT model.

2　Main developments

We now review the main developments of the AMPT model after the first public release of the AMPT source code in 2004 [60, 61] and the corresponding publication that described the physics details of the model at that time [13]. These developments are listed mostly in chronological order. In terms of the four main components of the AMPT model, Sections 2.2, 2.4, 2.5, 2.6, 2.7, 2.8 are about the initial condition, Sec. 2.9 is about the parton cascade, Sec. 2.3 is about the hadronization, while Sec. 2.1 is about the hadron cascade. Currently the public versions of the AMPT model since v1.26t5/v2.26t5 [61] have incorporated the changes made in the developments described in Sections 2.1 and 2.2; changes from the other developments will be released in the future.

2.1　 Deuteron productions in the hadron cascade

Light nuclei such as deuteron (d) and triton (t) are produced and observed in high energy nuclear collisions at RHIC and LHC [62, 63]. They have been proposed to be important for the search of the QCD critical point [64-67] and thus the study of light nuclei has become more active recently. Currently the production mechanism of light nuclei is still under debate, as there are several different models including the statistical model [68, 69], the nucleon coalescence model [70-74], and dynamical models based on the kinetic theory [75-77].

We have modified the AMPT model to provide a kinetic theory description of deuteron production and annihilation by adding the following reactions [77]:

B B^{'} \leftrightarrow M d,

(1)

where M = π, ρ, ω, and η, while B and B’ stand for baryons N, Δ, P₁₁(1440), and S₁₁(1535). Note that the hadron cascade component of the AMPT model [13], based on a relativistic transport (ART) model [84-86], already includes the interactions of π, K, η, ρ, ω, ϕ, K^*, N, Δ(1232), P₁₁(1440), S₁₁(1535) as well as their antiparticles. For the cross sections of the reactions $B B^{'} \to M d$ , we assume that their angular integrated mean squared matrix elements that are averaged over initial and summed over final spins and isospins are the same as that for the reaction $N N \to d π$ at same center of mass energy $\sqrt{s}$ . The cross sections for the inverse reactions $M d \to B B^{'}$ are then determined from the detailed balance. In addition to the production and annihilation processes for deuterons, we also include their elastic scatterings with mesons M and baryons B [77].

Experimentally, the cross sections for both the reaction $p p \to d π^{+}$ [78-80] and the reaction $π^{+} d \to p p$ [81-83, 87] have been extensively measured, and the former is given by

σ (p p \to d π^{+}) = \frac{1}{4} \frac{p_{π}}{p_{N}} f (s),

(2)

where p_N and p_π are, respectively, the magnitude of the three-momenta of initial and final particles in the center of mass frame. The function f(s), which is proportional to the angular integrated mean squared matrix elements that are summed over initial and final spins for the reaction $p p \to π^{+} d$ , is parameterized as

(3)

where $\sqrt{s}$ is in the unit of GeV and f(s) is in the unit of mb. For the inverse reaction $d π^{+}$ , its cross section is related to that for $p p \to d π^{+}$ via the detailed balance relation:

σ (d π^{+} \to p p) = \frac{2 p_{N}^{2}}{3 p_{π}^{2}} σ (p p \to d π^{+}) .

(4)

These parameterizations are compared with the experimental data in Fig. 1. The cross sections for the isospin averaged reactions $N N \to d π$ and $π d$ can then be obtained as $σ (N N \to d π) = 3 σ (p p^{+}) / 4$ and $σ (d π \to N N) = σ (d π^{+} \to p p)$ .

Fig. 1

(Color online) Experimental data on the total cross sections of

p p \to d π^{+}

[78-80] (filled symbols) and

d π^{+} \to p p

[81-83] (open symbols) in comparison with our parameterizations (solid curves).

We have coupled the above deuteron transport with an initial hadron distribution after hadronization as parameterized by a blast wave model [77], where a nucleon coalescence model using the deuteron Wigner function [88] was also applied for comparison. We find that the transport model gives very similar deuteron p_T spectra as the coalescence model; however the elliptic flows from the two models are different. In particular, the transport model gives a deviation of the elliptic flow from the exact nucleon number scaling at relatively high p_T and agrees better with the measured data.

On the other hand, the deuteron yield obtained directly from the AMPT-SM model is typically much lower than the experimental data. This could be due to the assumed relation between the $B B^{'} \leftrightarrow M d$ and $p p \leftrightarrow d π$ cross sections, which can be further constrained by using the measured total π d cross section, or the lack of additional production channels such as $π n p \leftrightarrow π d$ [89]. The low yield could also be partly due to the assumption of no primordial deuteron formation from quark coalescence. There are also studies [73, 90] that applied the nucleon coalescence model to the kinetic freezeout nucleon distributions from the AMPT-SM model. It has been found that the resultant light nuclei yields depend sensitively on the freezeout surface, which is affected by both the partonic expansion and the hadronization (quark coalescence) criterion. The yields also depend on the coalescence function used for the light nuclei [90], especially for small collision systems where the suppression due to the light nuclei size [91] could be significant. Further improvements of the AMPT model regarding the deuteron cross sections, the parton phase, and the hadronization criterion will benefit the studies of light nuclei.

2.2　 String melting model to describe the bulk matter

The Lund string model [25] is used in both the default and string melting versions of the AMPT model. In the default AMPT model, minijet partons recombine with their parent strings after the parton cascade to hadronize via the Lund string model into primordial hadrons. In the AMPT-SM model, the primordial hadrons that would be produced from the excited Lund strings in the HIJING model are “melt” into primordial quarks and antiquarks. Therefore, the parameters in the Lund string model affect the AMPT model results. In the Lund model, one assumes that a string fragments into quark-antiquark pairs with a Gaussian distribution in transverse momentum. Hadrons are formed from these quarks and antiquarks, with the longitudinal momentum given by the Lund symmetric fragmentation function [92, 93]

f (z) \propto z^{- 1} {(1 - z)}^{a_{L}} e^{- b_{L} m_{T}^{2} / z} .

(5)

In the above, z represents the light-cone momentum fraction of the produced hadron with respect to that of the fragmenting string and m_T is the transverse mass of the hadron.

When using the HIJING values [23, 24] for the key Lund string fragmentation parameters, a_L=0.5 and b_L=0.9 GeV^-2, the default AMPT model works well for particle yields and p_T spectra in pp collisions at various energies. However, it gives too small a charged particle yield in central Pb+Pb collisions at the SPS energy of E_LAB=158A GeV [11, 22]. Instead, modified values of a_L=2.2 and b_L=0.5 GeV^-2 were needed to fit the charged particle yield and p_T spectra in Pb+Pb collisions at SPS. For heavy ion collisions at higher energies such as RHIC energies, the default version of the AMPT model with these parameter values also reasonably describes hadron dN/dη, dN/dy and the p_T spectra in heavy ion collisions, although it underestimates the elliptic flow [26].

On the other hand, the AMPT-SM model [26, 13], due to its dense parton phase and quark coalescence, reasonably describes the elliptic flow [26] and two-pion interferometry [94] in heavy ion collisions. However, the versions before 2015 [61] (i.e., before v2.26t5) could not reproduce well the hadron dN/dη, dN/dy and p_T spectra (when using the same Lund parameters as the default version). For example, they overestimated the charged particle yield and significantly underestimated the slopes of the p_T spectra [13]. In an earlier attempt to reproduce data in Pb+Pb collisions at LHC energies with the AMPT-SM model, the default HIJING values for the Lund string fragmentation parameters were used [95] together with the strong coupling constant αs=0.33 (instead of 0.47); there the model reasonably reproduced the yield and elliptic flow of charged particles but underestimated the p_T spectrum (except at low p_T).

It was later realized that this problem of the AMPT-SM model can be solved [27] by using a very small value for the Lund fragmentation parameter b_L together with an upper limit on strange quark productions. Then the AMPT-SM model can reasonably reproduce the pion and kaon yields, p_T spectra, and elliptic flows at low p_T (below ˜1.5 GeV/c) in central and semi-central Au+Au collisions at the RHIC energy of $\sqrt{s_{N N}} = 200$ GeV and Pb+Pb collisions at the LHC energy of 2.76 TeV [27]. In particular, we found that b_L=0.15 GeV^-2 is needed [27], which is much lower than in the value used in previous studies [11, 13, 22, 26, 95]. Note that, for a smaller b_L value, the effective string tension κ, as given by [22, 13]

κ \propto \frac{1}{b_{L} (2 + a_{L})},

(6)

is higher and thus gives a larger mean transverse momentum for the initial quarks after string melting. In addition, the AMPT model assumes that the relative production of strange to non-strange quarks increases with the effective string tension [22, 13]. This is because the quark-antiquark pair production from string fragmentation in the Lund model is based on the Schwinger mechanism [96], where the production probability is proportional to $\exp (- π m_{⊥}^{2} / κ)$ at transverse mass $m_{⊥}$ . As a result, the strange quark suppression relative to light quarks, $\exp [- π (m_{s}^{2} - m_{u}^{2}) / κ]$ , is reduced for a higher string tension. It is found that an upper limit of 0.40 on the relative production of strange to non-strange quarks is needed for the AMPT-SM model [27].

Figure 2 shows the AMPT-SM results of pions and kaons for central (b<3 fm) and mid-central (b=7.3 fm) [97] Au+Au events at 200A GeV as well as central (b<3.5 fm) and mid-central (b=7.8 fm) [98] Pb+Pb events at 2.76A TeV. Also plotted for comparisons are the corresponding data for 0-5% and 20-30% centralities on dN/dy [99-101] in panels (a) and (d), the p_T spectra at mid-rapidity for the 0-5% centrality in panels (c) and (e), and v₂{EP} at mid-rapidity for the 20-30% centrality in panels (d) and (f). We see good agreements between the model results and the dN/dy data in both central and mid-central events at RHIC and LHC energies. The value of 0.55 is used for a_L at the top RHIC energy, while the value of 0.30 is used at the LHC energy since it gives a slightly better fit of the ALICE data [101] than the value of 0.55. We also see that the model roughly reproduces the observed p_T spectra at mid-rapidity below ˜2 GeV/c. In addition, the AMPT-SM model roughly describes the pion and kaon elliptic flow data on v₂{EP} [102, 103] at low p_T.

Fig. 2

(Color online) AMPT-SM results for pions (upper panels) and kaons (lower panels) on dN/dy of (a) π⁺ and (d) K⁺ in central and mid-central collisions, p_T spectra dN/(2πp_Tdp_Tdy) in the unit of c²/GeV² of (b) π⁺ and (e) K⁺ at mid-rapidity in central collisions, and elliptic flow v₂{EP} of (c) charged pions and (f) charged kaons at mid-rapidity in mid-central collisions in comparison with the experimental data for central (0-5%) and/or mid-central (20-30%) Au+Au collisions at 200A GeV and Pb+Pb collisions at 2.76A TeV.

This choice of settings for the AMPT-SM model [27] reasonably and simultaneously reproduces the particle yield, p_T spectra and elliptic flow of the bulk matter in central and semi-central AA collisions. Therefore, it enables us to make more reliable studies, such as the calculation of the evolution of energy density, effective temperatures, and transverse flow of the parton phase [27], and comprehensive predictions for Pb+Pb collisions at the top LHC energy of 5.02 TeV [31].

An example of the 5.02 TeV predictions from the AMPT-SM version v2.26t5 [31] is shown in Fig. 3, where the results on the η dependence of elliptic flow are shown in panels (a) and (b) for two centralities and the results on the factorization ratio r₂(ηa,ηb) are shown in panels (c) to (f) for four centralities. We see that the AMPT-SM model reasonably reproduces the observed v₂(η) magnitudes and shapes at 15-25% and 25-50% centralities from CMS [104] (filled circles) and ATLAS [105] (open circles) for Pb+Pb collisions at 2.76 TeV and from PHOBOS [106] (open diamonds) for Au+Au collisions at 200 GeV. We also see that the AMPT results on the factorization ratio r₂(ηa,ηb) as a function of ηa at 2.76 TeV are rather consistent with the corresponding CMS data [107], similar to a study [108] that used the AMPT-SM model as the initial condition for an ideal (3+1)D hydrodynamics. Furthermore, the AMPT-SM results show that the longitudinal correlation is much suppressed in Au+Au collisions at 200 GeV but slightly enhanced in Pb+Pb collisions at 5.02 TeV. Note that the longitudinal correlation comes naturally in the AMPT-SM model since each excited string typically produces many initial partons over a finite η range. Therefore, the initial transverse spatial geometry of the parton matter including the event plane has a strong correlation over a finite η range, and through partonic and hadronic interactions, the azimuthal anisotropies vn will then develop longitudinal correlations.

Fig. 3

(Color online) AMPT-SM results on the η dependence of v₂ in comparison with data for (a) the 15-25% centrality and (b) the 25-50% centrality, and (c)-(f) AMPT-SM results on the factorization ratio r₂(ηa,ηb) as functions of ηa in comparison with the CMS data for different centralities.

We note that the AMPT model may not be reliable at higher p_T, as indicated by Fig. 2, since it lacks inelastic parton collisions [13, 109] and consequently the radiative energy loss of partons that is important for high p_T partons. In addition, the string melting AMPT model up to now uses quark coalescence to model the hadronization of all partons, while the hadronization of high p_T partons and partons far away from their coalescing partners should be treated differently, e.g., with independent fragmentation [110] or string fragmentation [111].

2.3　 Improved quark coalescence

After parton scatterings, a spatial quark coalescence model is used to describe the hadronization process in the AMPT-SM model. It combines a quark with a nearby antiquark to form a meson and combines three nearby quarks (or antiquarks) into a baryon (or an antibaryon). For quarks and antiquarks in an event, the original quark coalescence model in AMPT [13, 26, 27, 31] searches for a meson partner before searching for baryon or antibaryon partners. Specifically, each quark (or antiquark) has its default coalescence partner(s), which are just the one or two valence parton(s) from the decomposition of the quark’s parent hadron from the string melting process. Then for any available (i.e., not-yet-coalesced) quark (or antiquark) that originally comes from the decomposition of a meson, the quark coalescence model searches all available antiquarks (or quarks) and selects the closest one in distance (in the rest frame of the quark-antiquark system) as the new coalescence partner to form a meson. After these meson coalescences are all finished, for each remaining quark (or antiquark) the model searches all available quarks (or antiquarks) and selects the closest two in distance as the new coalescence partners to form a baryon (or an antibaryon). As a result, the total number of baryons in any event after quark coalescence is the same as the total number before. Similarly, the quark coalescence process also conserves the number of antibaryons and the number of mesons in an event.

However, this separate conservation of the numbers of baryons, antibaryons, and mesons through the quark coalescence for each event is unnecessary, because only conserved charges such as the number of net-baryons and the number of net-strangeness need to be conserved. Therefore, we improved the coalescence method [32, 112] by removing the constraint that forced the separate conservations. Specifically, for any available quark, the new coalescence model searches all available antiquarks and records the closest one in relative distance (denoted as dM) as the potential coalescence partner to form a meson. The model also searches all available quarks and records the closest one in distance as a potential coalescence partner to form a baryon, and then searches all other available quarks again and records the one that gives the smallest average distance (i.e. the average of the three relative distances among these three quarks in the rest frame of the three-quark system, denoted as dB) as the other potential coalescence partner to form a baryon.

In the general case where both the meson partner and baryon partners are available, the quark will form a meson or a baryon according to the following criteria [32]:

\begin{array}{l} d_{B} < d_{M} r_{B M} : form a baryon; \\ otherwise:form a meson . \end{array}

(7)

In the above, rBM is the new coalescence parameter, which controls the relative probability of a quark forming a baryon instead of forming a meson. Note that the same coalescence procedure is also applied to all antiquarks, and the above criteria are not needed when only the meson partner or baryon partners can be found for a parton. In the limit of $r_{B M} \to 0$ , there would be no antibaryon formation at all while the minimum number of baryons would be formed as a result of the conservation of the (positive) net-baryon number. On the other hand, in the limit of $r_{B M} \to \infty$ , there would be almost no meson formation; more specifically, only 0, 1, or 2 mesons would be formed depending on the remainder when dividing the total quark number in the event by three. As a result, the new quark coalescence allows a (anti)quark the freedom to form a meson or a (anti)baryon depending on the distance from the coalescence partner(s). This is a more physical picture; for example, if a subvolume of the dense matter is only made of quarks which total number is a multiple of three, it would hadronize to only baryons (with no mesons) as one would expect.

We take central Au+Au collisions at $\sqrt{s_{N N}} = 200$ GeV from the AMPT-SM model as an example to compare the old and new quark coalescence [32]. The same parton cross section is used so that the parton phase-space configuration just before quark coalescence is statistically the same for the old and new quark coalescence. Figure 4 shows the average coalescence time of partons in mesons and (anti)baryons as functions of the hadron rapidity y_H. We see that baryons and antibaryons in the new quark coalescence (curve with open circles) are now formed much earlier than before. This is because the old quark coalescence tends to form (anti)baryons late, since it searches for meson partners before (anti)baryon partners and a parton will be unavailable for (anti)baryon formation when it is already used for meson formation. In contrast, the new quark coalescence searches for the potential meson partner and (anti)baryon partners concurrently and then determines the hadron type to be formed, making the coalescence process more physical as well as more efficient. In addition, we see that mesons in the new quark coalescence (curve with filled circles) are also formed earlier than before, possibly because of the improved efficiency after giving partons the freedom to form either a meson or a (anti)baryon. Since the plotted coalescence time is in the center-of-mass frame of the AA collision, we would expect a cosh y_H dependence if the dense matter were boost-invariant. The dot-dashed curve in Fig. 4 represents a function that is proportional to cosh y_H, which qualitatively agrees with our model results.

Fig. 4

(Color online) The average coalescence time of partons in mesons and (anti)baryons as functions of the hadron rapidity from the new (curves with circles) and old (dashed curves) quark coalescence for central Au+Au collisions at 200 GeV; the dot-dashed curve represents a cosh y_H curve for comparison.

Therefore, the new quark coalescence is more efficient, especially for the formation of (anti)baryons, due to the freedom of a parton to form either a meson or a (anti)baryon. As a result, it leads to improvements in the descriptions of (anti)baryon observables from the AMPT-SM model [32, 33, 113]. Figure 5 shows the AMPT results (with rBM=0.61) on various antiparticle-to-particle ratios around mid-rapidity for central Au+Au collisions at 200 GeV [38, 114, 115] and Pb+Pb collisions at 2.76 TeV [101, 116, 117] in comparison with the experimental data at mid-rapidity. Both the data and model results here are for the 0-5% centrality except that Ω at 200 GeV corresponds to the 0-10% centrality. We see that the results from the new quark coalescence (solid curves) are generally consistent with the experimental data, while results from the old quark coalescence (dashed curves) severely overestimate the ratios for Ξ and Ω. In addition, the antibaryon-to-baryon ratios generally increase with the strangeness content in both the AMPT model and the data. This is consistent with models such as the ALCOR model [118], which predict that these ratios are sequentially higher by a multiplicative factor, the K⁺/K^- ratio. Since the K⁺/K^- ratio is usually slightly larger than one at high energies, we see that our results from the improved quark coalescence agree rather well with this expectation and with the experimental data.

Fig. 5

(Color online) Antiparticle-to-particle ratios around mid-rapidity for (a) central Au+Au collisions at 200 GeV and (b) central Pb+Pb collisions at 2.76 TeV from the new (solid curves) and old (dashed curves) quark coalescence in comparison with the experimental data.

On the other hand, the AMPT model with the improved quark coalescence [32, 119] still underestimates the $\bar{p} / p$ ratio in central Au+Au collisions at and below 200 GeV. We note that quark coalescence should be augmented with other hadronization mechanisms such as fragmentation [110, 111] for partons that cannot find nearby partners. This will also help avoid the potential violation of the second law of thermodynamics during the hadronization process [120], where whether the entropy decreases during a phase-space quark coalescence has been found to depend on details such as the duration of the mixed phase, volume expansion, and resonance decays [121]. Also note that the rBM value of 0.61 is found to reasonably reproduce the proton and antiproton yields of AA collisions in the AMPT model with the original parton distribution function and HIJING’s nuclear shadowing [32], while the preferred rBM value is 0.53 for light (u/d/s) hadrons [119, 122] and 1.0 for charm hadrons [122] in the AMPT model with modern parton distribution functions of nuclei.

The new quark coalescence is not only able to describe the dN/dy yields, p_T spectra, and elliptic flows of pions and kaons at low p_T, it also better describes the baryon observables in general, especially the p_T spectra of (anti)baryons and antibaryon-to-baryon ratios for Ξ and Ω. It has also been shown to qualitatively describe the near-side anticorrelation feature of baryon-baryon azimuthal correlations observed in small systems at the LHC [33, 113]. In addition, it can be easily extended to include individual rBM factors specific to given hadron species, e.g., to describe the enhanced multi-strange baryon productions in nuclear collisions [123]. The string melting AMPT model with the new quark coalescence thus provides a better overall description of the bulk matter in high-energy nuclear collisions.

2.4　 Importance of finite nuclear thickness at lower energies

For heavy ion collisions at lower energies, the thickness of the incoming projectile and target nuclei in the center-of-mass frame becomes larger due to the finite Lorentz contraction along the beam directions. Therefore, one needs to consider the finite nuclear thickness in dynamical models of heavy ion collisions at lower energies, which correspond to higher net-baryon densities. The finite nuclear thickness increases the longitudinal width of the created matter and thus will obviously affect the initial energy and net-baryon densities [124, 125]. Furthermore, it will lead to a significant time duration of the initial particle and energy production; therefore, one cannot use a fixed proper time to describe the initial condition for hydrodynamic-based models but use a dynamical initialization scheme [126, 127].

For a central collision of two identical nuclei of mass number A, it takes the following time for the two nuclei to completely cross each other in the center-of-mass frame:

d_{t} = \frac{2 R_{A}}{\sinh y_{c m}}

(8)

in the hard sphere model of the nucleus. In the above, RA is the hard-sphere radius of the nucleus and ycm is the rapidity of the projectile nucleus. For central Au+Au collisions at $\sqrt{s_{N N}} = 50$ GeV, for example, dt≈0.5 fm/c, which is comparable to the typical value of the parton formation time or hydrodynamics initial time when one takes RA=1.12 A^1/3 fm. Therefore, one may expect the effect from finite nuclear thickness to be significant for central Au+Au collisions at $\sqrt{s_{N N}} ≲ 50$ GeV, which is the focus energy range of the RHIC Beam Energy Scan (BES) program.

We have developed semi-analytical methods [124, 125] to include the finite nuclear thickness in the calculation of the initial energy density, which is crucial in determining the initial temperature (and net-baryon chemical potential at low energies) of the produced QGP. Traditionally, the Bjorken formula [128] has been the standard semi-analytical tool in estimating the initial energy density in the central rapidity region right after the two nuclei pass each other:

ϵ_{B j} (t) = \frac{1}{A_{T} t} \frac{d E_{T}}{d y} .

(9)

In the above, A_T represents the full transverse area of the overlap volume, and dE_T/dy is the initial rapidity density of the transverse energy at mid-rapidity, which is often approximated with the experimental dE_T/dy value in the final state. Because the Bjorken energy density diverges as $t \to 0$ , a finite value is needed for the initial time, which is often taken as the proper formation time of the produced quanta τ_f. However, a serious limitation of the Bjorken formula results from the fact that it neglects the finite thickness of the colliding nuclei. Therefore, one expects that the Bjorken formula may break down when the crossing time is not small compared to the formation time [6].

Using the semi-analytical methods that include the finite nuclear thickness, the initial energy density ϵ(t) averaged over the transverse area of the overlap region as a function of time, including its maximum value ϵ^max, have been calculated [124, 125]. We first considered the finite time duration of the initial energy production but neglected the finite longitudinal extension [124], which enabled us to obtain explicit analytical solutions of ϵ(t). Both the uniform time profile and beta time profile have been considered, where in the uniform time profile one assumes that the initial transverse energy at y ≈ 0 is produced uniformly in time (x) from t₂ to t₂:

\frac{d^{2} E_{_{T}}}{d y d x} = 1 t_{2} - t_{1} \frac{d E_{_{T}}}{d y}, i f x \in [t_{1}, t_{2}] .

(10)

In contrast, the beta time profile assumes the following:

\frac{d^{2} E_{_{T}}}{d y d x} \propto {[x (d_{t} - x)]}^{n} \frac{d E_{_{T}}}{d y}, i f x \in [0, d_{t}] .

(11)

Note that n=4 is chosen [124] from the comparison to the time profile of partons within mid-spacetime-rapidity in central Au+Au collisions from the AMPT-SM model. In addition, for the uniform profile shown here, t₁=0.29dt & t₂=0.71 dt are used since they give the same mean and standard deviation of time as the beta profile at n=4.

Then we considered both the finite time duration and longitudinal extension of the initial energy production [125]. When τ_f is not too much smaller than the crossing time of the two nuclei, results from the latter study [125] are similar to those from the former study [124]. On the other hand, there is a qualitative difference in that the maximum energy density ϵ^max at τ_f=0 is finite after considering both the finite duration time and longitudinal extension [125], while the Bjorken formula diverges as 1/τ_f and the method that only considered the finite duration time [124] diverges as ln (1/τ_f) at low energies but as 1/τ_f at high energies. Overall, these studies have yielded the following qualitative conclusions: the initial energy density after considering the finite nuclear thickness approaches the Bjorken formula at high colliding energies and/or large formation time τ_f. At low colliding energies, however, the initial energy density has a much lower maximum, evolves much longer, and is much less sensitive to τ_f than the Bjorken formula. Note that we have written a web interface [129] that performs the semi-analytical calculation [125] for central AA collisions, where the user can input the colliding system, energy and the proper formation time.

To include the finite nuclear thickness, we have modified the initial condition of the AMPT-SM model [124] to specify in each heavy ion event the longitudinal coordinate z₀ and time t₀ of each excited string, which is then converted into the initial partons via string melting. Note that in the normal AMPT-SM model [26, 13, 27, 32] the longitudinal coordinate z₀ and time t₀ of each excited string in the initial state are both set to zero, which would be correct only at very high energies. Figure 6 shows the results on the time evolution of the average energy density at ηs ≈ 0 for central Au+Au collisions at four different energies for formation time τ_f=0.1 and 0.3 fm/c. At the high energy of 200 GeV, the AMPT-SM results with (curves with filled circles) and without (curves with open circles) the finite nuclear thickness are essentially the same. This is consistent with the fact that the Bjorken result and our semi-analytical result are also very similar (after shifting the results in time); it also confirms the expectation that the finite nuclear thickness can be mostly neglected at high-enough energies.

Fig. 6

(Color online) Average parton energy densities at central spacetime-rapidity from AMPT with (filled circles) and without (open circles) the finite nuclear width for central Au+Au collisions and τ_f=0.1 & 0.3 fm/c at (a) 4.84, (b) 11.5, (c) 27, and (d) 200 GeV; corresponding results for the uniform time profile (dashed curves), the beta time profile (solid curves), and the Bjorken formula (dot-dashed curves) are also plotted for comparison.

At lower energies, however, the AMPT results after including the finite nuclear thickness are very different: the maximum energy density is lower, and the time evolution of the energy density (e.g., the time spent above half the maximum energy density) is longer. These key features agree with the semi-analytical results [124, 125], where the results from the uniform time profile and the more realistic beta time profile are close to each other after the uniform profile is set to the same mean and standard deviation of time as the beta profile [124]. We also see from Fig. 6 that the increase of the maximum initial energy density with the colliding energy is much faster after including the finite nuclear thickness, which is consistent with the analytical finding that the Bjorken formula overestimates the maximum energy density more at lower energies [124, 125]. In addition, we see in Fig. 6 that the AMPT results are generally wider in time; partly because the parton proper formation time in AMPT is not set as a constant but is inversely proportional to the parent hadron transverse mass [13]. Secondary parton scatterings and the transverse expansion of the dense matter in AMPT can also cause differences from the semi-analytical results, which do not consider such effects. Overall, we see that the AMPT results without considering the finite nuclear thickness are similar to the Bjorken results, while the AMPT results including the finite thickness are similar to our semi-analytical results. These results suggest that it is important to include the finite nuclear thickness in dynamical models of relativistic heavy ion collisions, especially at lower energies.

2.5　 Modern parton distribution functions in nuclei

The initial condition of the AMPT model is based on the HIJING two-component model [23]. The primary interactions between the two incoming nuclei are divided into two parts: the soft component described by the Lund string fragmentation model [25, 92, 93], and the hard component with minijet productions described by perturbative QCD through the PYTHIA5 program [25].

The minijet differential cross sections in HIJING model can be computed using the factorization theorem in the perturbative QCD framework [130] as

\frac{d σ_{j e t}^{c d}}{d p_{T}^{2} d y_{1} d y_{2}} = K \sum_{a, b} x_{1} f_{a} (x_{1}, Q^{2}) x_{2} f_{b} (x_{2}, Q^{2}) \frac{d σ^{a b \to c d}}{d \hat{t}} .

(12)

In the above, p_T is the transverse momentum of the produced minijet parton, y₁ and y₂ are the rapidities of the two produced partons c and d, the factor K accounts for higher-order corrections to the leading order jet cross section, x₁ and x₂ are respectively the momentum fraction x carried by the two initial partons, f_a(x₁, Q₂) is the parton distribution function (PDF) of parton type a at x=x₁ and factorization scale Q², σab is the parton-parton cross section for partons a and b, and $\hat{t}$ is the standard Mandelstam variable for the minijet production subprocess.

The total inclusive jet cross section (for the production of minijet gluons and u/d/s quarks and antiquarks) is then obtained by integrating the above differential cross section with a transverse momentum cutoff p₀ and considering all the possible combinations of final state parton flavors [23]:

σ_{j e t} = \frac{1}{2} \sum_{c, d} \int_{p_{0}^{2}}^{\hat{s} / 4} d p_{T}^{2} d y_{1} d y_{2} \frac{d σ_{j e t}^{c d}}{d p_{T}^{2} d y_{1} d y_{2}},

(13)

where $\hat{s}$ is the standard Mandelstam variable for the minijet subprocess. We see that the minijet transverse momentum cutoff p₀ and the parton distribution function fa(x,Q²) are the key factors affecting the jet cross section. The total jet cross section together with the σ_soft parameter that describes the soft component determine the nucleon-nucleon interaction cross sections in the Eikonal formalism [131, 132]. Note that p₀ is only relevant when the center-of-mass energy per nucleon pair is higher than 10 GeV, because the jet production in the HIJING model is switched off at $\sqrt{s_{N N}} < 10$ GeV.

An important ingredient needed in Monte Carlo event generators for hadron collisions is the input parton distribution function [133-135]. Efforts have been made to understand the subtilties of implementing various parton distributions for phenomenological studies based on event generators[136, 137]. The impacts of different parton distributions in the event generators for pp collisions are found to be sizable and the key parameters in the generators usually depend on the details of the input PDF [138]. Specifically, the parton distribution function in the AMPT model affects the initial state radiation and the minijet production within the two-component model framework. Using modern parton distributions along with fine tuned model parameters is required to generate reliable exclusive final states in the AMPT model.

The HIJING 1.0 model [23, 24] that generates the initial condition of the original AMPT model employs the Duke-Owens parton distribution function set 1 [139] for the free proton. However, it is well known that the Duke-Owens PDFs were obtained at a time when a large array of experimental data used in the global fittings for modern PDFs were not yet available [135]. The parton densities at small-x relevant for minijet and heavy flavor productions at high energies are generally underestimated by the Duke-Owens PDFs [140]. Therefore, we have updated the AMPT model with a modern parton distribution function of the free nucleon (the CTEQ6.1M set [141]) and retuning of the relevant p₀ and σ_soft parameters [119]. Note that this update is based on the AMPT model with the new quark coalescence [32]. Also note that the HIJING 2.0 model [142] is a similar update, which replaces the Duke-Owens PDFs in the HIJING 1.0 model with the GRV PDFs [143].

For nuclear collisions at sufficiently high energies, results from event generators depend on the parton distribution functions inside the incoming nuclei. Analogous to the free nucleon case, global analyses of the parameterizations for the modifications of the nuclear PDFs relative to the free nucleon PDFs have been performed by several groups [144-148]. In addition, it is natural to expect the nuclear modification to depend on a nucleon’s position inside a nucleus. Therefore, the spatial dependence of nuclear parton densities are considered [23, 149-154], and a global analysis to extract the nuclear PDFs with spatial dependence is carried out [155] based on the EKS98 [156] and EPS09 [157] fits. In a recent study [119], we have included the spatially dependent EPS09s nuclear modifications [155] in the AMPT model to replace the original HIJING 1.0 nuclear shadowing parameterization. Note that the HIJING 1.0 shadowing is spatially dependent but independent of Q² or the parton flavor [13, 23], similar to the HIJING 2.0 nuclear shadowing [142].

For a proton bound in a nucleus, its parton distribution function of flavor i can be written as

f_{i}^{p / A} (x, Q^{2}) \equiv R_{i}^{A} (x, Q^{2}) f_{i}^{p} (x, Q^{2}),

(14)

where $f_{i}^{p} (x, Q^{2})$ is the corresponding PDF in the free proton. Here $R_{i}^{A} (x, Q^{2})$ represents the spatially averaged nuclear modification function, which typically depends on the x range: the shadowing region at small x, an anti-shadowing region at x ˜0.1, and the EMC effect region at x close to 1. The spatial dependence of the nuclear modification function can be formulated as

R_{i}^{A} (x, Q^{2}) \equiv \frac{1}{A} \int d^{2} s T_{A} (s) r_{i}^{A} (x, Q^{2}, s),

(15)

where TA(s) represents the nuclear thickness function at transverse position s, and $r_{i}^{A} (x, Q^{2}, s)$ represents the spatially dependent nuclear modification function.

Solid curves in Fig. 7 show the gluon density distributions (multiplied by x) in the free proton from the original and the updated AMPT model. The gluon density distributions of a bound proton at the center of a lead nucleus from the EPS09s and HIJING nuclear modifications are also shown in Fig. 7. We see that in the updated AMPT model with the CTEQ6.1M set the gluon densities are quite different from the old Duke-Owens set and much higher at small x. We also see that the gluon shadowing in EPS09s is much weaker than that in the HIJING 1.0 model.

Fig. 7

(Color online) Gluon density distributions (multiplied by x) in free proton (solid curves) and proton inside the lead nucleus (dashed curves) at Q²=10 GeV² from the original (black) and updated (red) AMPT.

As mentioned earlier, p₀ and σ_soft are the two key parameters in the HIJING model that determine the total and inelastic cross sections of pp and $p \bar{p}$ collisions within the Eikonal formalism. In the HIJING 1.0 model that uses the Duke-Owens PDFs, constant values of p₀ = 2.0 GeV/c and σ_soft = 57 mb are found to reasonably describe the experimental cross sections of pp and $p \bar{p}$ collisions over a wide energy range [23, 24, 158]. On the other hand, when the PDFs are updated in the HIJING 2.0 model [119, 142], energy-dependent p₀(s) and (s) values are needed.

In our work that implements the CTEQ6.1M PDF in the AMPT model [119], the energy dependent parameters p₀(s) and (s) are determined via fitting the experimental total and inelastic cross sections of pp and $p \bar{p}$ collisions within the energy range $4 < \sqrt{s} {< 10}^{5}$ GeV, as shown in Fig. 8. We then obtain the following p₀(s) and (s) functions for pp collisions:

Fig. 8

(Color online) Total and elastic cross sections versus the colliding energy of pp collisions from the updated AMPT model (solid and dot-dashed curves) in comparison with the experimental data (symbols); σ_jet from the model is also shown (dotted curve).

\begin{array}{l} p_{0}^{p p} (s) = - 1.71 + 1.63 \ln (\sqrt{s}) - 0.256 \ln^{2} (\sqrt{s}) \\ + 0.0167 \ln^{3} (\sqrt{s}), \end{array}

(16)

\begin{array}{l} σ_{s o f t} (s) = 45.1 + 0.718 \ln (\sqrt{s}) + 0.144 \ln^{2} (\sqrt{s}) \\ + 0.0185 \ln^{3} (\sqrt{s}) . \end{array}

(17)

In the above, $p_{0}^{p p}$ and the center-of-mass colliding energy $\sqrt{s}$ are in the unit of GeV/c and GeV, respectively; while σ_soft is in the unit of mb. Note that $\sqrt{s}$ will be replaced with $\sqrt{s_{N N}}$ for nuclear collisions. We also find that the updated AMPT-SM model reasonably describes the charged particle yield and p_T spectrum in pp and/or $p \bar{p}$ collisions from $\sqrt{s} \sim$ 4 GeV to 13 TeV [119].

When we apply the above p₀(s), (s) and the EPS09s nuclear shadowing to central AA collisions at LHC energies, however, we find rather surprisingly that the hadron yields from the AMPT-SM model are significantly higher than the experimental data. As shown in Fig. 9(a), the AMPT-SM model that uses the $p_{0}^{p p} (s)$ value overestimates the final state hadron multiplicities in central Pb+Pb collisions at $\sqrt{s_{N N}} = 2.76$ TeV. Since a larger p₀ value would suppress the total jet cross section and reduce the particle yields, we introduce a global scaling of the minijet cutoff p₀ to make its value in central AA collisions, $p_{0}^{A A} (s)$ , nuclear size dependent [119]:

Fig. 9

(Color online) Identified particle rapidity distributions in 0-5% central Pb+Pb collisions at

\sqrt{s_{N N}} = 2.76

TeV from the AMPT-SM model using the minijet cutoff (a)

p_{0}^{p p} (s)

or (b)

p_{0}^{A A} (s)

in comparison with the ALICE data [99, 101].

\begin{array}{l} p_{0}^{A A} (s) = p_{0}^{p p} (s) A^{q (s)}, \\ q (s) = 0.0334 \ln (\sqrt{s} E_{0}) - 0.00232 \ln^{2} (\sqrt{s} E_{0}) \\ + 0.0000541 \ln^{3} (\sqrt{s} E_{0}), f o r \sqrt{s} \geq E_{0} = 200 G e V . \end{array}

(18)

The above q(s) function is determined from the fit to the overall particle yields of central Au+Au collisions at the RHIC energies and central Pb+Pb collisions at LHC energies (see [119] for details). Its value is zero at $\sqrt{s_{N N}} \leq 200$ GeV since the $p_{0}^{p p} (s)$ value works reasonably well there for central Au+Au collisions, while its value approaches 0.16 at $\sqrt{s_{N N}} \sim 10^{7}$ GeV. This nuclear scaling of the minijet momentum cutoff scale p₀ is motivated by the physics of the color glass condensate [159], where the saturation momentum scale Qs depends on the nuclear size as $Q_{s} \propto A^{1 / 6}$ in the saturation regime for small-x gluons in AA collisions at high-enough energies.

As shown in Fig. 9(b), the updated AMPT-SM model using $p_{0}^{A A} (s)$ from the global nuclear scaling well reproduces the identified particle yields in central Pb+Pb collisions at the LHC energy. In addition, we find that a significantly smaller value for the Lund bL parameter, bL=0.15 GeV^-2, is needed to describe the particle p_T spectra in central AA collisions [119], similar to an earlier study [27]. Note that recently we have generalized the minijet cutoff p₀ and the Lund bL parameter with a local nuclear scaling [160], as shall be discussed in Sec. 2.7, which would help explain why a bigger p₀ value but a smaller bL value are needed for high energy AA collisions than pp collisions.

2.6　 Improvements of heavy flavor productions

Heavy flavors are predominantly produced from the initial hard scatterings at early times in nuclear collisions [161-163]. Therefore, they are powerful observables to probe the strong electromagnetic field created in heavy ion collisions [164-166] and transport properties of the dense matter [167-171]. Multiple theoretical frameworks have been developed for the description of open heavy flavor productions in high energy pp and pA collisions based on the pQCD framework [172-175]. Medium effects, such as those from pQCD calculations of the parton energy loss [176, 177] or the Langevin/Boltzmann equation methods [178-186] can be included for AA collisions.

Study of heavy flavor productions within the AMPT model [187, 188] has the potential to provide a unified model for both light and heavy flavor transport and improve our understanding of the non-equilibrium effects of the QGP evolution [49, 50, 189, 190]. In addition, using parton scatterings to model the interactions between heavy quarks and the evolving medium, the parton cascade approach is able to implement any scattering angular distribution without the need to assume small-angle scatterings. Therefore, after the update with modern parton distributions for proton and nuclei as discussed in Sec. 2.5, we have made several significant improvements of heavy flavor productions in the AMPT model [122]. First, for self consistency we include the heavy flavor cross sections in the total minijet cross section in the HIJING two-component model. Second, we remove the minimum transverse momentum requirement (p₀) for initial heavy quark productions since the heavy quark pair production cross sections from pQCD are already finite due to the heavy quark mass. These changes can be illustrated with the following modified formula for the minijet cross section [122]:

\begin{array}{l} σ_{jet} = ​ ​ ​ ​ \sum_{c, d} \frac{1}{1 + δ_{c d}} \int_{p_{0}^{2}}^{\hat{s} / 4} d p_{T}^{2} d y_{1} d y_{2} \frac{d σ_{light}^{c d}}{d p_{T}^{2} d y_{1} d y_{2}} \\ + \sum_{c, d} \int_{0}^{\hat{s} / 4} d p_{T}^{2} d y_{1} d y_{2} \frac{d σ_{heavy}^{c d}}{d p_{T}^{2} d y_{1} d y_{2}}, \end{array}

(19)

where the first term on the right hand side represents the cross section of light flavor (u/d/s/g) minijets and the second term represents that of heavy flavors such as charm and bottom. Note that the factor 1/(1+δcd) above becomes 1/2 for final states with identical partons, such as $g + g \to g + g$ for minijet gluon productions. In contrast, the original HIJING model uses Eq.(13) and applies the factor of 1/2 to all light flavor minijet production processes [23], and that leads to a smaller σ_jet than Eq.(19) (at the same p₀). As a result, an increase of the p₀ value as shown below is needed [122] to describe the experimental data on total and inelastic cross sections of pp and $p \bar{p}$ collisions shown in Fig. 8:

\begin{array}{l} p_{0}^{p p} (s) = - 1.92 + 1.77 \ln (\sqrt{s}) - 0.274 \ln^{2} (\sqrt{s}) \\ + 0.0176 \ln^{3} (\sqrt{s}) \end{array}

(20)

with p₀ in GeV/c and $\sqrt{s}$ in GeV.

The total $c \bar{c}$ cross section for pp collisions from the updated AMPT model (solid curve) is shown in Fig. 10 versus the colliding energy in comparison with the available world data. We see that the updated AMPT model can well describe the data in pp collisions at various collision energies. The original AMPT model (dotted curve), however, significantly underestimates the charm quark yield, especially at low energies. The enhanced charm quark production in the updated model is largely due to the removal of the p₀ cut, since the charm quark cross section is much lower when charm quarks in the updated AMPT model are required to have a transverse momentum above p₀ (dashed curve).

Fig. 10

(Color online) Total cross sections of charm-anticharm quark pairs from the AMPT model for pp collisions in comparison with the world data [191-197] as functions of the colliding energy.

Figure 11 shows the charm quark rapidity and transverse momentum distributions from the AMPT model for pp collisions at $\sqrt{s} = 200$ GeV and 7 TeV. Note that the charm quark or hadron results shown in this section have been averaged over those for particles and the corresponding antiparticles. As shown in Fig. 11(a), the charm quark yield in the updated AMPT model is found to be significantly higher than that in the original AMPT model over the full rapidity range at both RHIC and LHC energies. From the charm quark p_T spectra at mid-rapidity shown in Fig. 11(b), we see that the removal of the p₀ cut for charm quarks mostly enhances charm quark productions at low p_T. We also see that results from the original AMPT model (dotted curve) and the updated AMPT model that includes the p₀ cut (dashed curve) are similar, partly because a p₀ cut (2 GeV/c) on the charm quark production is also used in the original AMPT model.

Fig. 11

(Color online) Charm-anticharm quark pair (a) rapidity distributions and (b) p_T spectra around mid-rapidity in pp collisions at

\sqrt{s} = 200

GeV (red curves) and 7 TeV (black curves) from the AMPT model. The shaded band represents statistical errors.

In AA collisions, heavy quark production is subject to additional medium-induced initial state and final state effects. Within the AMPT model, initial state effects include the nuclear modification of the parton distribution functions in nuclei, while the final state effects are mostly treated with parton elastic rescatterings in the parton cascade [109]. Figure 12 shows the charm quark yield at mid-rapidity for 0-10% central Au+Au or 0-10% central Pb+Pb collisions at different energies. We see that the EPS09s nuclear modification leads to an enhancement of the charm quark yield in central AA collisions at lower energies but a suppression at high energies. This is expected due to the anti-shadowing feature at large x and the shadowing feature at small x in the nuclear modification functions. We also see that the result from the updated AMPT model (solid curve) is in good agreement with the charm quark data, which is obtained for 0-10% central Au+Au collisions at $\sqrt{s_{N N}} = 200$ GeV by scaling the STAR pp charm quark cross section data with the number of binary collisions [198, 199]. Similar to the results for pp collisions, the updated AMPT model gives a significantly higher charm quark yield at mid-rapidity in central AA collisions compared to the original AMPT model.

Fig. 12

(Color online) dN/dy of charm quark pairs around mid-rapidity from the AMPT model for 0-10% central Au+Au collisions at RHIC energies and Pb+Pb collisions above RHIC energies as functions of the colliding energy in comparison with the STAR data.

The open heavy flavor hadron species formed by quark coalescence include charm and bottom hadrons with all possible charges. To reproduce the observed vector to pseudo-scalar meson ratios of open heavy flavors in pp collisions, we fit the relative probability of forming primordial vector versus pseudo-scalar heavy mesons in the quark coalescence model, e.g., the ratio is set to 1.0 for the primordial D^*/D and B^*/B ratios [122], instead of using only the invariant mass of the coalescing partons in the original AMPT-SM model [13]. Note that only the primordial ratios right after coalescence are specified with these parameters, not the vector to pseudo-scalar meson ratios in the final state which include effects from resonance decays. In addition, in the new quark coalescence model [32] that is used in this heavy flavor work [122], the relative probability for a quark to form a baryon instead of a meson is determined by the rBM parameter as shown in Eq.(7). In our earlier work that updated the AMPT model with modern PDFs [119], the rBM value for light flavor (u/d/s) hadrons is set to 0.53, which value is also used here. On the other hand, we set the rBM value for heavy flavor hadrons to 1.0, because using the light flavor value would lead to too few charm baryons (by a factor of ˜4) compared to the experimental data in pp or AA collisions. In principle, the rBM value for charm hadrons depends on properties such as the number and masses of available charm baryon states versus charm meson states. The necessity of using a higher rBM value for charm is consistent with the assumption that there are more charm baryon states than charm meson states compared to the light flavor sectors [200].

After the improvements on heavy flavor productions, we find that the updated AMPT model [119] can well describe the yields and p_T spectra of open charm hadrons including D, D^*, Ds and Λc in pp collisions at different energies. The updated model also describes the charm data in central AA collisions much better than the original AMPT model. However, the updated AMPT model still does not well describe the charm hadron productions in AA collisions [122]. As shown in Fig. 13(a), the updated AMPT model overestimates the D⁰ yield at low p_T but underestimates it at p_T above 2.5 GeV/c when compared to the STAR data for 0-10% Au+Au collisions at $\sqrt{s_{N N}} = 200$ GeV [198]. Compared to the original AMPT model, the charm hadron yield from the updated AMPT model is significantly enhanced at low p_T, similar to the results at the parton level shown in Fig. 11(b). We also find that the charm baryon to meson ratio (Λc/D) in 0-10% Au+Au collisions at 200 GeV is much larger in the updated AMPT model (dashed curve) than the original AMPT model (dotted curve), as shown in Fig. 13(b). Compared to the STAR data for 10-80% Au+Au collisions at 200 GeV [201], result of the Λc/D ratio versus p_T from the updated AMPT model is somewhat higher. We also see that this ratio from the AMPT model is slightly lower in the 10-80% centrality than the 0-10% centrality for Au+Au collisions at 200 GeV.

Fig. 13

(Color online)(a) p_T spectra of open charm hadrons in central Au+Au collisions at

\sqrt{s_{N N}} = 200

GeV and (b) the Λc/D ratios versus p_T for two centralities of Au+Au collisions at 200 GeV from the AMPT model in comparison with the STAR data [198, 201].

We note that only elastic parton scatterings are included in the AMPT-SM model; therefore, the model is only applicable in the region where the effect of parton radiative energy loss is small. Studies have suggested that the elastic collisional energy loss could be dominant for charm hadrons below p_T ˜5-6 GeV/c in Au+Au collisions at $\sqrt{s_{N N}}$ = 200 GeV or below p_T ˜15 GeV/c in Pb+Pb collisions at $\sqrt{s_{N N}}$ = 2.76TeV [184]. Therefore, the charm results from the AMPT model are not reliable at p_T higher than these values. Also note that the charm p_T spectra are affected by the charm quark scattering cross section and its angular distribution in ZPC. The AMPT model currently uses the $g + g \to g + g$ cross section for scatterings of all parton flavors, where flavor-dependent cross sections and angular distributions should be used for the parton scatterings. In addition, there is still a large uncertainty on the nuclear shadowing of gluons [155], which has not been fully explored in the AMPT model. Furthermore, hadronic scatterings of heavy flavor hadrons [202-206] have not been included in the AMPT model except for the decays of heavy hadron resonances. Further developments of the AMPT model are expected to improve its description of heavy flavor productions in AA collisions.

2.7　 System size dependence under local nuclear scaling

The system size dependence of observables can be useful to uncover the transition of certain phenomena in nuclear collisions, such as the onset of collectivity and whether it comes from initial state momentum correlations [54, 55] or final state interactions [49-51, 56, 57]. It is known from multiple studies that certain key parameters in the initial condition of the AMPT model for AA collisions need to be different from their values for pp collisions to reasonably describe the data [13, 22, 27, 32, 95, 119]. First, the Lund b_L parameter in the symmetric string fragmentation function [92, 93], as shown in Eq.(5), for large collision systems needs to be significantly smaller than its value for pp collisions. An earlier study has also shown that a constant b_L can not describe the centrality dependence of 〈p_T〉 in heavy ion collisions [31], where the system size dependence of the Lund fragmentation parameters was suggested as a possible solution. Note that similar frameworks for the system size dependence have been implemented in the string fragmentation model [207-211]. Second, we have found in earlier developments of the AMPT model [119, 122] that the minijet transverse momentum cutoff p₀ for central Pb+Pb collisions at the LHC energies needs to be significantly higher than that for pp collisions at the same energy. These observations suggest that the above two parameters should be related to the size of the colliding system to provide better initial conditions for the AMPT model.

Therefore, we have recently proposed [160] that the b_L and p₀ parameters in AMPT can be considered as local variables that depend on the nuclear thickness functions of the two incoming nuclei. This prescription allows us to use the parameter values obtained for pp collisions and the local nuclear scaling relations to obtain the values for AA collisions; the model would then describe the system size and centrality dependences of nuclear collisions self-consistently.

In the Lund string model [92, 93], the symmetric fragmentation function is given by Eq.(5). The average squared transverse momentum of massless hadrons is related to the Lund fragmentation parameters a_L and b_L as [13]

〈 p_{T}^{2} 〉 = \frac{1}{b_{L} (2 + a_{L})} .

(21)

Consequently, the 〈p_T〉of partons after string melting and the final hadrons are significantly affected by the value of bL. Since the mean transverse momentum of initial partons in heavy ion collisions is expected to be higher in larger systems due to the higher initial temperature, we expect the bL value to decrease with the system size. Note that the string tension is believed to be larger in a denser matter [209, 210, 212, 213], thus a decrease of b_L with the system size is consistent with the expectation of a stronger color field and thus a higher string tension κ since κα/ 1=b_L[13] as shown in Eq.(6).

b_{L} (s_{A}, s_{B}, s) = \frac{b_{L}^{p p}}{{[\sqrt{T_{A} (s_{A}) T_{B} (s_{B})} / T_{p}]}^{β (s)}} .

(22)

In the above, $b_{L}^{p p}$ is the value for pp collisions (chosen to be 0.7 GeV^-2 based on the fit of 〈p_T〉 data), s represents the square of the center-of-mass collision energy per nucleon pair, $T_{A} (s_{A}) = \int ρ_{A} (s_{A}, z) d z$ is the nuclear thickness function at the transverse distance sA from the center of nucleus A determined with the Woods-Saxon nuclear density profiles [156], and Tp is the average value of the effective thickness function of the proton (taken as 0.22 fm^-2). Note that Tp is used instead of TA(sA) or TB(sB) when the projectile or the target is proton or when TA(sA) or TB(sB) from the nucleus is smaller than the T_p value. Although different mathematical forms from that of Eq.(22) can be used in the local scaling relation, our study [160] shows that the geometric scaling functions (i.e., using the geometric mean of the two nuclear thickness functions) generally works better than the arithmetic form. We note that a systematic Bayesian analysis based on the TRENTo initial condition [214] with a hybrid model found that the geometric form for the local nuclear scaling is preferred by the experimental data [215].

The exponent function β(s) describes the energy dependence of the local nuclear scaling of b_L. From the fits to charged particle 〈p_T〉 data in the most central Au+Au collisions at RHIC energies and most central Pb+Pb collisions at LHC energies, it is parameterized as

β (s) = 0.620 + 0.112 \ln (\frac{\sqrt{s}}{E_{0}}) Θ (\sqrt{s} - E_{0}),

(23)

where E₀=200 GeV and Θ(x) is the unit step function. The fitted β(s) function is shown in Fig. 14(a) (dashed curve), which is a constant at RHIC energies but grows rapidly at LHC energies. Note that β=1 at high energies (dotted line) may be a “natural” limit for Eq.(22) if we imagine that all local strings fully overlap so that the string tension adds up. That would give $b_{L} \propto 1 / T_{A} (s_{A})$ for central AA collisions, where TA(sA) is proportional to the local number of participant nucleons or excited strings integrated over the longitudinal length.

Fig. 14

(Color online)(a) Fitted exponent functions α(s), β(s), and 3q(s) versus the center-of-mass energy per nucleon pair

\sqrt{s}

. (b) Average p₀ and bL values versus the impact parameter of pPb, Au+Au and Pb+Pb collisions at several energies.

Figure 14(b) shows the b_L value averaged over the overlap volume versus the impact parameter for Pb+Pb and pPb collisions at 5.02A TeV and Au+Au collisions at two RHIC energies. We see that 〈b_L〉 for Pb+Pb collisions at the LHC energy is lower than that for Au+Au collisions at RHIC energies, which corresponds to a larger string tension due to the larger value of the exponent β(s) at LHC energies. On the other hand, the impact parameter dependences of 〈b_L〉 at different RHIC energies are essentially the same since β(s) is a constant within that energy range. For pPb collisions at 5.02A TeV, its 〈b_L〉 is higher than that in Pb+Pb collisions at small b and grows faster with b due to its smaller system size.

The minimum transverse momentum cutoff p₀ for light flavor minijet productions is another key parameter in the HIJING model and thus in the initial condition of the AMPT model [23, 119, 142]. In our update of the AMPT model with modern nPDFs [119], the collision energy dependence of p₀ is determined from fitting the pp cross section data. Then motivated by the physics of the color glass condensate [159], a global nuclear scaling of the p₀ cutoff [119] has been introduced for central AA collisions above the top RHIC energy of 200A GeV to describe the experimental data on charged particle yields in central Pb+Pb collisions at LHC energies. Here [160] we go beyond the global nuclear scaling and instead consider p₀ as a local variable that depends on the transverse position of the corresponding hard process in each event. As p₀ is expected to increase with the system size, we related its value to the nuclear thickness functions in a general AB collision as [160]

p_{0} (s_{A}, s_{B}, s) = p_{0}^{p p} (s) {[\sqrt{T_{A} (s_{A}) T_{B} (s_{B})} / T_{p}]}^{α (s)} .

(24)

As $T_{A} (s) \propto A^{1 / 3}$ , Eq.(23) approximately gives $p_{0} \propto A^{α (s) / 3}$ for central AA collisions and thus essentially recovers the global nuclear scaling if α(s)=3q(s). On the other hand, for peripheral collisions where TA(sA) and TB(sB) are very small and thus replaced with the proton value (T_p), Eq.(23) automatically gives the p₀ value for pp collisions.

Since $p_{0}^{p p} (s)$ works for charged particle yields in central Au+Au collisions at and below 200A GeV, we assume that the need to modify p₀ in nuclear collisions starts at the top RHIC energy [119]. From the comparison to charged particle yields in the most central Pb+Pb collisions at 2.76A and 5.02A TeV, we obtain the preferred α(s) values at those two energies. We then fit the α(s) function as [160]

\begin{array}{l} α (s) = 0.0918 \ln (\frac{\sqrt{s}}{E_{0}}) - 0.00602 \ln^{2} (\frac{\sqrt{s}}{E_{0}}) \\ + 0.000134 \ln^{3} (\frac{\sqrt{s}}{E_{0}}), f o r \sqrt{s} \geq E_{0}, \end{array}

(25)

with α(s)=0 for $\sqrt{s} < E_{0} = 200$ GeV. We see in Fig. 14(a) that α(s) ≈ 3q(s) as expected. We also see that both approach the value of 1/2 at very high energies; this is consistent with our expectation [119] that p₀ is closely related to the saturation momentum Qs in the color glass condensate [159], where $Q_{s} \propto A^{1 / 6}$ in the saturation regime.

Figure 14(b) also shows the average p₀ value as a function of the impact parameter for Pb+Pb collisions at 2.76A TeV and 5.02A TeV as well as pPb collisions at 5.02A TeV. As expected, we see that 〈p₀〉 decreases with the impact parameter and that 〈p₀〉 at the lower LHC energy is smaller than that at the higher LHC energy due to the smaller α(s) value. Also, 〈p₀〉 in pPb collisions is smaller than that in Pb+Pb collisions at the same colliding energy due to its smaller size. In addition, the relative variation of 〈p₀〉 with the impact parameter is seen to be much weaker than that of 〈b_L〉 since $α (s) ≪ β (s)$ for the exponents in the local nuclear scaling relations.

We show in Fig. 15 the dN_ch/dη yield in panel (a) and charged particle 〈p_T〉 in panel (b) around mid-pseudorapidity versus centrality from different AMPT versions in comparison with the experimental data for Au+Au collisions at 200A GeV and Pb+Pb collisions at 5.02A TeV [216-222]. Using the local nuclear scaling, the improved AMPT model (solid curves) reasonably describes these centrality dependence data in AA collisions at both RHIC and the LHC energies, with a significant improvement in the 〈p_T〉 description as shown in Fig. 15(b). When we switch off the local nuclear scaling of p₀ and b_L but instead use constant b_L=0.15 GeV^-2 and p₀(s) (constant at a given energy), we recover the AMPT model developed earlier [122] and obtain the dot-dashed curves when using $p_{0} (s) = p_{0}^{A A} (s)$ and the dotted curves when using $p_{0} (s) = p_{0}^{p p} (s)$ . They both give the wrong centrality dependence of 〈p_T〉, since the model results (dot-dashed or dotted) show a mostly increasing trend from central to peripheral collisions while the data show a mostly decreasing trend. Results from the public AMPT version 2.26t9 [61] are also shown (dashed curves) [27], which also fail to describe the centrality dependence of charged particle 〈p_T〉 data.

Fig. 15

(Color online)(a) dN_ch/dη and (b) 〈p_T〉 around mid-pseudorapidity versus the centrality of 5.02A TeV Pb+Pb collisions (thick curves) and 200A GeV Au+Au collisions (thin curves) from this work (solid curves) and earlier AMPT versions in comparison with the experimental data (symbols). The p_T range used for the 〈p_T〉 calculation is ［0.15, 2］ GeV/c at 5.02A TeV and ［0.2, 2］ GeV/c at 200A GeV.

In Fig. 15(a), we see that the charged particle yield in central Pb+Pb collisions at 5.02A TeV is significantly overestimated when using $p_{0} (s) = p_{0}^{p p} (s)$ , and the global nuclear scaling $p_{0} (s) = p_{0}^{A A} (s)$ is needed to reproduce the particle yield. From the Pb+Pb results, we also see that the effect from the global nuclear scaling of p₀ in peripheral collisions is much smaller than that in central collisions, because the binary scaling of minijet productions makes p₀ less important for peripheral collisions. It is thus not surprising to see that the dN_ch/dη results from the local nuclear scaling are similar to the AMPT results using the constant $p_{0}^{A A}$ for central collisions but close to the AMPT results using the constant $p_{0}^{p p}$ for peripheral collisions.

The local nuclear scaling relations also predict how observables depend on the system size going from large to small systems. Figures 16(a) and 16(b) show respectively the dN_ch/dη and charged particle 〈p_T〉 around mid-pseudorapidity from the AMPT model [160] versus centrality in comparison with the experimental data for Au+Au collisions and several smaller collision systems [221, 223-226]. We see that the improved AMPT model describes these data rather well, further demonstrating the validity of the local nuclear scaling assumption. Note that, although the mid-pseudorapidity dN_ch/dη and 〈p_T〉 data for the most central Au+Au/Pb+Pb collisions have been used in the determination of the parameter functions α(s) and β(s), the data of these smaller systems are not considered in the fitting of the parameters. Results for Au+Au collisions at 200A GeV are also presented in Fig. 16, and we see that the changes of the charge particle yield and 〈p_T〉 from Cu+Cu to Au+Au collisions at 200A GeV are well accounted for by the local nuclear scaling. For example, the 〈p_T〉 in Cu+Cu is generally smaller than that in Au+Au due to the larger b_L value for Cu+Cu collisions. Note however that our calculations here have not considered the deformation of the Xe nucleus [227].

Fig. 16

(Color online) (a) dN_ch/dη and (b) 〈p_T〉 around mid-pseudorapidity versus the centrality of Xe+Xe collisions at 5.44A TeV, Cu+Cu collisions at 200A GeV, and pPb collisions at 5.02A TeV from the AMPT model (curves) in comparison with the experimental data (symbols). The p_T range used for the 〈p_T〉 calculation is ［0.15, 2］ GeV/c at LHC energies and ［0.2, 2］ GeV/c at 200A GeV.

2.8　 PYTHIA8 initial condition with sub-nucleon structure

The modifications of the AMPT initial condition discussed so far have been performed within the framework of the HIJING two-component model that uses the PYTHIA5 program. While the development of local nuclear scaling [160] enables the AMPT model to reproduce the system size dependence and centrality dependence of changed particle yields and 〈p_T〉 in pA and AA collisions using the parameter values for minimum bias pp collisions, we have not directly addressed the multiplicity dependence of these observables, especially the 〈p_T〉, in pp collisions. On the other hand, PYTHIA8 [228] is quite successful in describing the particle yields in pp collisions. It has been extended to treat pA or AA collisions based on the Angantyr framework [229], and PYTHIA8 has been used as the initial condition generator for multiple heavy ion Monte Carlo models [230-232]. Therefore, it is worthwhile to have the option to use PYTHIA8 as the initial condition for the AMPT model.

Recently we have coupled PYTHIA8 with the final state parton and hadron interactions and quark coalescence [32] of the AMPT-SM model to study pp collisions [233]. In this approach, the fluctuating initial condition of AMPT originally provided by the HIJING model is replaced by the PYTHIA/Angantyr model [229]. In addition, the sub-nucleon structure, which could be important for collectivity observables in small systems [234-238], can be modeled when implementing the space-time structure of the string system generated by the PYTHIA multiparton interaction. With the proton charge distribution given by

ρ (r) = \frac{1}{8 π R^{3}} e^{- r / R}

(26)

with R=0.2 fm, the sub-nucleon spatial structure can be related to the transverse positions of the excited strings in two ways. In the first way, the transverse coordinates of the produced string objects are sampled according to the overlap function of a pp collision at a given impact parameter b:

T (x, y, b) = \int ρ (x - b / 2, y, z) ρ (x + b / 2, y, z) d z,

(27)

where z is along the beam directions. In the second way, the initial transverse spatial condition including event-by-event sub-nucleon fluctuations is generated with a Glauber Monte Carlo method based on the constituent quark picture [236, 239-243]. By modeling the proton as three constituent quarks, the interaction of two protons can be interpreted as collisions between the constituent quarks from each incoming proton within the Glauber model framework [241, 244]. The positions of the quark constituents are first sampled with the proton profile ρ(r), then the transverse coordinates of the excited strings are randomly assigned to the binary collision center of each interacting constituent pair.

Figure 17 shows the effect of using PYTHIA8 as the AMPT initial condition on the identified particle 〈p_T〉 versus the charge particle pseudo-rapidity density in pp collisions at $\sqrt{s} = 13$ TeV. Note that only hadrons within 0<p_T<3 GeV/c and |y|<0.5 are included in this comparison, and the central values of ALICE data are obtained with a refit to the data [245]. We see that this AMPT model (solid curves), which uses the PYTHIA8 initial condition and includes both parton and hadron evolutions, roughly reproduces the experimental data. On the other hand, the original AMPT model (dashed curves) reasonably describes the pion 〈p_T〉 but gives a very weak multiplicity dependence for the proton 〈p_T〉. The significant improvement compared to the original AMPT model on the multiplicity dependence of the proton 〈p_T〉 presumably results from multiparton interaction in the PYTHIA8 model.

Fig. 17

(Color online) 〈p_T〉 of π (black), K (blue) and proton (red) at mid-rapidity within 0<p_T<3 GeV/c versus the charged hadron multiplicity density in 13 TeV pp collisions. The AMPT model using the PYTHIA8 initial condition (solid curves) are compared to the original AMPT model (dashed curves) and the ALICE data.

Figure 18(a) shows the average initial spatial eccentricity of partons in the transverse plane right after string melting as a function of the parton multiplicity of each event from the two ways of generating the sub-nucleon spatial structure. Note that only partons with formation time less than 5 fm/c are considered, and eccentricities are calculated with the initial position of each parton at its formation time [246]. When using the overlap function weighting method (black curves), the eccentricity is largely driven by the geometric shape of the transverse overlap area and thus decreases significantly with the parton multiplicity as shown in panel (a) and increases significantly with the impact parameter as shown in panel (b). On the other hand, when using the Monte Carlo method with constituent quarks (red curves), large eccentricities in the initial condition can be generated even in very central collisions or events at high multiplicities. Figure 18(b) actually shows that the initial eccentricity from the constituent quark method is larger for pp collisions at smaller impact parameters, opposite to the behavior from the overlap function method.

Fig. 18

(Color online) The initial eccentricity of partons right after string melting versus (a) the number of partons and (b) the impact parameter, and (c) two-particle long-range angular correlations for events at two different multiplicity classes, for pp collisions at 13 TeV. AMPT results with the sub-nucleon structure are shown for the overlap function method (black curves) and the MC method with constituent quarks (red curves).

The difference in the initial spatial eccentricity could certainly affect final state momentum anisotropies in small collision systems after interactions in the AMPT model convert the spatial anisotropies into momentum anisotropies [26, 49, 50]. Using the AMPT model with PYTHIA8 as the initial condition, we have found [233] that two-particle long-range correlations in high multiplicity pp collisions at the LHC depend sensitively on how the sub-nucleon structure of the proton is implemented. We analyze the projected correlation function of two charged hadrons with a large pseudorapidity gap:

C (Δ ϕ) = \frac{1}{N_{trig}} \frac{d N^{pair}}{d Δ ϕ} .

(28)

Both trigger and associate hadrons are required to be within 1<p_T<3 GeV/c and |η|<2.4 following the analysis procedure of the CMS Collaboration [247], and the two hadrons in each pair must be separated in pseudo-rapidity with a gap |Δη|>2. Events are separated into two categories based on N_sel, the number of selected charge tracks with p_T>0.4 GeV/c and |η|<2.4. High multiplicity events are defined as those with N_sel>80, while low multiplicity events are defined as those with N_sel<20.

Figure 18(c) shows the multiplicity dependence of the C(Δϕ) function from the two ways of generating the sub-nucleon spatial structure for 0.2 mb parton cross section [233]. We see that the AMPT model using PYTHIA8 shows a long-range ridge-like structure for high multiplicity events when the proton geometry is modeled with the constituent quark method (red solid curve), while the overlap function weighting method (black solid curve) does not show this structure. This demonstrates the connection between two-particle long-range correlations with the underlying sub-nucleon structure and fluctuations. Note that a significant near-side ridge structure in the correlation function is found in the experimental data, which has been regarded as an important signature of collectivity in high multiplicity pp events [44, 247].

We note that the original AMPT-SM model also shows the long-range near-side correlations, although it does not include the sub-nucleon structure [233]. In addition, the PYTHIA event generator itself has considered final state hadronic rescatterings [206, 248-250]. Using the AMPT-SM model with PYTHIA8 initial conditions, we can extend the study of pp collisions [233] to pA and AA collisions with the Angantyr model within the PYTHIA8 framework. That would lay a solid foundation for the studies of different mechanisms of collectivity, such as string shoving and parton/hadron evolutions, with the same model.

2.9　 Improved algorithm for the parton cascade

Particle correlations and momentum anisotropies in the AMPT-SM model are usually dominated by parton interactions [13, 26, 41]. We have also found that even a few parton scatterings in a small system is enough to generate significant momentum anisotropies through the parton escape mechanism [49, 50]. It is therefore important to ensure that the parton cascade solution in the AMPT model is accurate.

The ZPC elastic parton cascade [109] in the AMPT model solves the Boltzmann equation by the cascade method, where a scattering happens when the closest distance between two partons is less than the range of interaction $\sqrt{σ_{p} / π}$ with σ_p being the parton scattering cross section. The default differential cross section in ZPC for two-parton scatterings, based on the gluon elastic scattering cross section as calculated with QCD at leading order, is given by [13, 109]

\frac{d σ_{p}}{d \hat{t}} = \frac{9 π α_{s}^{2}}{2} (1 + \frac{μ^{2}}{\hat{s}}) \frac{1}{{(\hat{t} - μ^{2})}^{2}},

(29)

where μ is a screening mass to regular the total cross section. This way the total cross section has no explicit dependence on $\hat{s}$ :

σ_{p} = \frac{9 π α_{s}^{2}}{2 μ^{2}} .

(30)

The above Eqs.(28-29) represent forward-angle scatterings. For isotropic scatterings, $d σ_{p} / d \hat{t}$ is independent of the scattering angle.

It is well known that cascade calculations suffer from the causality violation [251, 252] due to the geometrical interpretation of cross section. This leads to inaccurate numerical results at high densities and/or large scattering cross sections (i.e., large opacities). For example, a recent study [29] has shown that the effect of causality violation on the elliptic flow from the AMPT-SM model [13] is small but non-zero. Causality violation also leads to the fact that different choices of performing collisions and/or the reference frame can lead to different numerical results [253-255]. These numerical artifacts due to the causality violation can be reduced or removed by the parton subdivision method [12, 43, 252, 254, 256–260]. However, parton subdivision usually alters the event-by-event correlations and fluctuations, the importance of which has been more appreciated in recent years [34]; it is also much more computationally expensive. Therefore, it is preferred to improve the parton cascade to yield solutions that are accurate enough without using parton subdivision. We have recently pursued this goal for box calculations [261].

In ZPC, one can take different choices or collision schemes to implement the cascade method [109]. With the closest approach criterion for parton scatterings, the closest approach distance is usually calculated in the two-parton center of mass frame. Two partons may collide when their closest approach distance is smaller than $\sqrt{σ_{p} / π}$ , and at a given global time all such possible collisions in the future are ordered in a collision list with the ordering time of each collision so that they can be carried out sequentially. The collision list is updated continuously after each collision, and for expansion cases the parton system dynamically freezes out when the collision list is empty. For calculations of a parton system in a box, we terminate the parton cascade at a global time that is large enough so that the parton momentum distribution changes little afterwards. When the closest approach distance is calculated in the two-parton center of mass frame, the collision time of a scattering is a well-defined single value. However, because of the finite σ_p the two partons have different spatial coordinates in general; thus this collision time in the two-parton center of mass frame becomes two different collision times in the global frame (named here as ct₁ and ct₂ respectively for the two colliding partons) after the Lorentz transformation. The default collision scheme of ZPC [109] uses (ct₁+ct₂)/2 as both the collision time and ordering time; this is the case for the AMPT model [13]

Results from the default ZPC scheme [261] at σ_p=2.6 mb are shown in Fig. 19 (curves with open circles). Panel (a) shows the final parton p_T distribution, while panels (b) and (c) show the time evolution of parton 〈p_T〉 (scaled by T) and variance of p_T (scaled by T²), respectively. The gluon system is initialized in a box with an off-equilibrium initial momentum distribution as shown by the dot-dashed curve in panel (a), where the gluon density is set the same as that for a thermalized gluon system with the Boltzmann distribution at temperature T=0.5 GeV. We see from Fig. 19(a) that the final distribution from the default ZPC scheme deviates considerably from the expected thermal distribution (dotted curve). On the other hand, we find that a new collision scheme, which uses (ct₁,ct₂) as both the collision time and ordering time, gives a final distribution very close to the thermal distribution [261]. The causality violation usually suppresses collision rates, which is the case for the default ZPC scheme; it is therefore understandable that choosing time (ct₁,ct₂) instead of (ct₁+ct₂)/2 enhances the collision rates and thus suppresses the causality violation.

Fig. 19

(Color online)(a) The final p_T distribution, (b) time evolution of 〈p_T /T〉, and (c) time evolution of 〈p_T variance〉/T² from ZPC from the default collision scheme (open circles) and the new collision scheme (filled squares) in comparison with parton subdivision results at subdivision factor l=10⁶ (dashed curves) for gluons in a box at T=0.5 GeV and σ_p=2.6 mb.

We use the parton subdivision method to obtain the “exact” time evolutions of 〈p_T〉 and p_T variance (dashed curves) in Figs. 19(b) and (c). We see that the time evolution of the p_T variance from the default scheme deviates significantly from the “exact” parton subdivision result, although the time evolutions of 〈p_T〉 are close to each other (mostly due to the conservation of total momentum). In contrast, the time evolution of the p_T variance from the new scheme [261] is very close to the parton subdivision result, which at late times agrees with theoretical expectation (diamond). By examining cases of different parton densities and cross sections [261], we find rather surprisingly that the new scheme for ZPC gives very accurate results (i.e., very close to parton subdivision results and/or theoretical values) even at very large opacities, such as the case of T=0.7 GeV and σ_p=10 mb.

We have used a novel parton subdivision method for the results shown in Fig. 19. In the standard method, one increases the initial parton number per event by factor l while decreasing the cross section by the same factor, which can be schematically represented by the following:

N \to l \times N, V unchanged,

(31)

where N is the initial parton number in an event and V is the initial volume of the parton system. Since the number of possible collisions scales with l², the subdivision method is very expensive in terms of the computation time, which roughly scales with l² per subdivision event or l per simulated parton. However, for box calculations where the density function f(x,p,t) is spatially homogeneous, the following new subdivision method can be used:

N unchanged, V \to V / l,

(32)

where we decrease the volume of the box by factor l while keeping the same parton number and momentum distribution in each event. This subdivision method is much more efficient than the standard subdivision method; we therefore use a huge subdivision factor 10⁶ (instead of the usual value of up to a few hundreds).

We emphasize that the differential cross section must not be changed when performing parton subdivision; as a result, the exact transformation for parton subdivision is [261]

f (x, p, t) \to l \times f (x, p, t), \frac{d σ_{p}}{d \hat{t}} \to \frac{d σ_{p}}{d \hat{t}} / l .

(33)

This is especially relevant for forward-angle scattering. For example, when parton subdivision requires the decrease of the forward-angle cross section of Eq.(29), one should not do that by increasing the screening mass μ by a factor of $\sqrt{l}$ because that would change the angular distribution of the scatterings in Eq.(28). Instead, one can decrease the αs parameter by a factor of $\sqrt{l}$ , which decreases the total scattering cross section while keeping its angular distribution the same.

Transport coefficients such as the shear viscosity η represent important properties of the created matter [262]. Therefore, we have also evaluated the effect of the new collision scheme on the shear viscosity η and its ratio over the entropy density η/s. The Green-Kubo relation [263, 264] has been applied [265-269] to calculate the shear viscosity at or near equilibrium. We thus start with an equilibrium initial condition for shear viscosity calculations according to the Green-Kubo relation [265].

Figure 20 shows our η/s results as functions of the opacity parameter χ, which is defined as [254]

Fig. 20

(Color online) The η/s ratio for different cases of gluon scatterings in a box versus the opacity χ; the solid curve without symbols represents the Navier-Stokes expectation for isotropic scatterings.

χ = \sqrt{\frac{σ_{p}}{π}} / λ = n \sqrt{\frac{σ_{p}^{3}}{π}},

(34)

where n is the parton density and λ is the mean free path. The case shown in Fig. 19 for gluons in a box at T=0.5 GeV and σ_p=2.6 mb then corresponds to χ=2.0, and other χ values shown in Fig. 20 are obtained for the following cases: T=0.2 GeV and σ_p=2.6 mb, T=0.7 GeV and σ_p=5.2 mb, and T=0.7 GeV and σ_p=10 mb [261]. For isotropic scatterings of a massless Maxwell-Boltzmann gluon gas in equilibrium (where s=4n and degeneracy factor dg=16), we have the following Navier-Stokes expectation:

{(\frac{η}{s})}^{N ​ S} ≃ \frac{0.4633}{d_{g}^{1 / 3} χ^{2 / 3}} = \frac{0.1839}{χ^{2 / 3}},

(35)

which only depends on the opacity χ. We see in Fig. 20 that for isotropic scatterings the subdivision result agrees well with the Navier-Stokes expectation (solid curve). On the other hand, the extracted η and η/s values from the default ZPC scheme are significantly different from the Navier-Stokes expectation or the parton subdivision results at large opacities, although they agree at low opacities as expected. We also see that the results from the new collision scheme are very close to the subdivision results for both forward-angle scatterings and isotropic scatterings, even at a huge opacity χ=41. The new ZPC collision scheme is the first step towards the validation and improvement of the ZPC parton cascade for scatterings in 3-dimensional expansion cases.

3　Other developments

There are other developments of the AMPT model that have not been covered in the previous section. Here we gave a brief overview of some of these works.

The AMPT model has been extended to include deformed nuclei as the projectile and/or target. First, deformed uranium nuclei are implemented [270] to study various observables in U+U collisions at 200A GeV and the effect of nuclear deformation. Later, the AMPT model has been modified to specify the initial proton and neutron spatial distributions in the ⁹⁶Ru or ⁹⁶Zr nucleus according to the density functional theory (DFT) calculations [271-273]. The effects of the DFT nuclear density distributions on the backgrounds and possible signals of the chiral magnetic effect (CME) in isobar collisions are then investigated [271]. The extended AMPT model is also used in the study that proposes a novel method to search for the CME in a single heavy ion collision system [272]. Another study [273] uses the model to study multiplicity distributions and elliptic flow in isobar collisions, where the differences between the two isobar systems have the potential to decisively discriminate DFT nuclear distributions from the usual Woods-Saxon density distributions.

The AMPT model has also been extended to include mean field potentials in the hadronic phase in a study of the elliptic flow splitting of particles and antiparticles at the RHIC BES energies [274]. A later study coupled the AMPT model with a parton transport based on the 3-flavor Nambu-Jona-Lasinio model [275] to include the partonic mean field potentials; it shows that a combination of partonic and hadronic mean field potentials can describe the observed splitting of elliptic flows.

The current AMPT model has been known to violate the electric charge conservation because of two reasons [276]. First, the hadron cascade is based on the ART model [84] that has K⁺ and K^- as explicit particles but not K⁰ or ${\bar{K}}^{0}$ . As a result, we change K⁰ to K⁺ and change $^{0}$ to K^- prior to the hadron cascade in order to include hadronic interactions of neutral kaons, and after the hadron cascade we assume the isospin symmetry and thus change half of the final K⁺ into K⁰ and change half of the final K^- to $^{0}$ . The second reason is that many hadron reactions and some resonance decays in AMPT violate the electric charge conservation. Some reaction channels do not consider electric charges of the initial-state hadrons; instead the isospin-averaged cross section is used and the electric charge of each final state hadron is set randomly [276]. We have developed a version of the AMPT model that has corrected these problems and thus satisfies the electric charge conservation [277]. This charge-conserved version of the AMPT model has been shared with some colleagues for their recent studies of charge-dependent CME signals [278, 279].

Recently we have developed a pure hadron cascade version of the AMPT model (AMPT-HC) [280] to study heavy ion collisions at low energies below a few GeVs. Note that the Eikonal formalism, which is a basis of the HIJING model and thus the initial condition of the standard AMPT model, is expected to break down for nuclear collisions at low enough energies. We thus treat a heavy ion collision as individual nucleon-nucleon collisions in the AMPT-HC model. First, we use the Woods-Saxon nucleon density distribution and the local Thomas-Fermi approximation to initialize the position and momentum of each nucleon in the incoming nuclei. Then primary nucleon-nucleon collisions are treated with the hadron cascade component of AMPT, without going through the Lund string fragmentation, the parton cascade, or quark coalescence. In addition to the usual elastic and inelastic collisions, the hadron cascade in the AMPT-HC model also includes hadron mean field potentials for kaons, baryons and antibaryons. This model has been used to study the $Ξ^{-}$ production in low energy Au+Au collisions, which is proposed as a better probe of the nuclear equation of state at high densities than single strangeness (kaon or Λ) productions [280].

4　Summary and outlook

A multi-phase transport model was constructed to provide a self-contained kinetic theory-based description of relativistic nuclear collisions with its four main components: the fluctuating initial condition, partonic interactions, hadronization, and hadronic interactions. Here we review the main developments since the public release of the AMPT source code in 2004 and the 2005 publication that described the details of the model at that time. Several developments have been carried out to improve the initial condition, including the incorporation of finite nuclear thickness relevant for heavy ion collisions below the energy of tens of GeVs, the incorporation of modern parton distribution functions of nuclei for high energy heavy ion collisions, improvement of heavy quark productions, the use of local nuclear scaling of key input parameters for the system size dependence and centrality dependence, and the incorporation of PYTHIA8 and sub-nucleon structure. There are also ongoing efforts to improve the accuracy of the parton cascade without using the parton subdivision method that would alter event-by-event correlations and fluctuations. In addition, the spatial quark coalescence model has been further developed to allow a quark the freedom to form either a meson or a baryon depending on the distance to its coalescing partner(s), which improves baryon and antibaryon productions of the model. Furthermore, deuteron production and annihilation processes have been included in the hadron cascade, an AMPT version that satisfies the electric charge conservation has been developed, and a pure hadron cascade version of the AMPT model is recently developed to study heavy ion collisions at low energies below a few GeVs. For high energy nuclear collisions where the quark-gluon plasma is expected, the string melting version of the AMPT model can now reasonably and simultaneously describe the yield, transverse momentum spectrum and elliptic flow of the bulk matter from small to large collision systems. Consequently, the AMPT model has been applied to the study of various observables in nuclear collisions such as particle yields, particle correlations and anisotropic flows, vorticity and polarization.

Because the transport model approach can address non-equilibrium dynamics, it provides a complementary framework to hydrodynamical models for large systems at high energies, and more importantly it is well suited to study the transition from the dilute limit to the hydrodynamic limit. Therefore, it will be worthwhile to further develop a multi-phase transport as a dynamical model for relativistic nuclear collisions.

There are multiple areas that should be addressed in the future. Regarding the initial condition, at low enough energies the pure hadron cascade version should be applicable while at high enough energies the Eikonal formalism should be valid. It would be desirable to have a unified physics formulation that self-consistently changes from one regime to the other as the colliding energy increases. In addition, for high enough energies and/or large enough collision systems the QGP is expected to be formed, and consequently the string melting version of the AMPT model should be applicable instead the string-dominated default version. The AMPT model should be improved to dynamically determine whether the QGP should be formed in the initial state; it would then self-consistently change from a string-dominated initial condition to a parton-dominated one when the initial energy density is high enough. Another deficiency in the initial condition of the string melting AMPT model is the lack of gluons in the parton phase, and the color-glass-condensate approach would be ideal for including initial gluons once the approach can be generalized to address the quark degrees of freedom such as the nonzero net-baryon number. Regarding the parton phase, the parton cascade should be generalized to perform transport in the presence of an electromagnetic field to enable studies of the electromagnetic field and related observables. Another area of development concerns the study of high net-baryon density physics and the QCD critical point. The AMPT model could be coupled to or improved with effective theories such as the functional renormalization group method or the Nambu-Jona-Lasinio model to treat parton interactions self-consistently including the effective equation of state and effects from the critical point. Regarding the hadronization process, a dynamical parton recombination criterion, e.g., by using the local parton energy density as the recombination criterion instead of starting hadronization at the parton kinetic freezeout, should be developed. Also, additional mechanisms such as independent fragmentation should be included to treat partons that do not find suitable coalescence partners within the local phase space; this would enable the AMPT model to be applicable to studies of high p_T physics once the radiative energy loss of high p_T partons is considered in the parton phase. Regarding the hadron cascade, it can benefit from the inclusion of more resonances for more realistic thermodynamic properties and chemical equilibration of the hadron matter, and modern models such as the SMASH model could be a good choice as the new hadron cascade component. We expect that the AMPT model in the near future, even if only improved in a few focused areas, will enable us to address some key questions in heavy ion physics and also serve as a more reliable open source transport model for the community for various studies of nuclear collisions.

References

[1]

T. D. Lee and G. C. Wick,

Vacuum stability and vacuum excitation in a spin 0 field theory

. Phys. Rev. D 9, 2291-2316 (1974). doi: 10.1103/PhysRevD.9.2291

Further developments of a multi-phase transport model for relativistic nuclear collisions

1 Introduction

2 Main developments

2.1 Deuteron productions in the hadron cascade

2.2 String melting model to describe the bulk matter

2.3 Improved quark coalescence

2.4 Importance of finite nuclear thickness at lower energies

2.5 Modern parton distribution functions in nuclei

2.6 Improvements of heavy flavor productions

2.7 System size dependence under local nuclear scaling

2.8 PYTHIA8 initial condition with sub-nucleon structure

2.9 Improved algorithm for the parton cascade

3 Other developments

4 Summary and outlook