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

Zi-Wei Lin· Liang Zheng

Abstract 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.

Keywords QGP · Transport model · Heavy-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)twocomponent 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 (pT) spectra of identified particles observed in heavy ion collisions at SPS and RHIC. However, it significantly underestimates the elliptic flow (v2) 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/3He+Au collisions at RHIC [45, 46]. Hydrodynamic calculations seem to describe the experimental data well [47, 48]. The AMPTSM 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, we review in Sect. 2 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 Sect. 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 Sect. 2.2,the new quark coalescence model in Sect. 2.3, incorporation of the finite nuclear thickness along beam directions in Sect. 2.4, incorporation of modern parton distribution functions of nuclei in Sect. 2.5,improved treatment of heavy quark productions in Sect. 2.6, the introduction of local nuclear scaling of key input parameters to describe the system size dependence in Sect. 2.7, incorporation of PYTHIA8 and nucleon substructure in the initial condition in Sect. 2.8, and benchmark and improvement of the parton cascade algorithm in Sect. 2.9.We then briefly review other developments of the AMPT model in Sect. 3.Finally,in Sect. 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, Sects. 2.2, 2.4, 2.5, 2.6, 2.7, 2.8 are about the initial condition, Sect. 2.9 is about the parton cascade,Sect. 2.3 is about the hadronization, while Sect. 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 Sects. 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 that describe the data, 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]:

Experimentally, the cross sections for both the reaction pp →dπ+[78-80]and the reaction π+d →pp[81-83,87]have been extensively measured, and the former is given by

where pNand 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 pp →π+d, is parameterized as

These parameterizations are compared with the experimental data in Fig. 1. The cross sections for the isospin averaged reactions NN →dπ and πd →NN can then be obtained as σ(NN →dπ)=3σ(pp →dπ+)/4 and σ(dπ →NN)=σ(dπ+→pp).

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 pTspectra 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 pTand 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 BB′?Md and pp ?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 πnp ?π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.

Fig.1 (Color online)Experimental data on the total cross sections of pp →dπ+ [78-80] (filled symbols) and dπ+ →pp [81-83] (open symbols) in comparison with our parameterizations (solid curves)

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 AMPTSM 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]

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

When using the HIJING values [23, 24] for the key Lund string fragmentation parameters, aL=0.5 and bL=0.9 GeV-2, the default AMPT model works well for particle yields and pTspectra 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 ELAB=158A GeV [11, 22]. Instead, modified values of aL=2.2 and bL=0.5 GeV-2were needed to fit the charged particle yield and pTspectra 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 pTspectra in heavy ion collisions,although it underestimates the elliptic flow[26].

On the other hand,the AMPT-SM model[13,26],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 pTspectra (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 pTspectra[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 pTspectrum (except at low pT).

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[13,22].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(-πm2⊥/κ)at transverse mass m⊥. As a result, the strange quark suppression relative to light quarks, exp[-π(m2s-m2u)/κ], is reduced for a higher string tension.It is found that an upper limit of 0.40 on the relative production of strange to nonstrange quarks is needed for the AMPT-SM model [27].

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 midcentral collisions,pT spectra dN/(2πpTdpTdy)in the unit of c2/GeV2 of b π+ and e K+ at mid-rapidity in central collisions, and elliptic flow v2{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

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 pTspectra at mid-rapidity for the 0-5% centrality in panels (b) and (e),and v2{EP} at mid-rapidity for the 20-30% centrality in panels (c) and (f). We see good agreements between the model results and the dN/dy data in both central and midcentral events at RHIC and LHC energies. The value of 0.55 is used for aLat 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 pTspectra at mid-rapidity below ～2 GeV/c. In addition, the AMPT-SM model roughly describes the pion and kaon elliptic flow data on v2{EP}[102,103]at low pT.

This choice of settings for the AMPT-SM model [27]reasonably and simultaneously reproduces the particle yield, pTspectra and elliptic flow of the bulk matter in central and semi-central AA collisions at high energies.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].

Fig.3 (Color online)AMPT-SM results on the η dependence of v2 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 r2(ηa，ηb) as functions of ηa in comparison with the CMS data for different centralities

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 r2(ηa，ηb) are shown in panels (c) to(f) for four centralities. We see that the AMPT-SM model reasonably reproduces the observed v2(η) 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 r2(ηa，ηb) as a function of ηaat 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 vnwill then develop longitudinal correlations.

We note that the AMPT model may not be reliable at higher pT, as indicated by Fig. 2, since it lacks inelastic parton collisions [13, 109] and consequently the radiative parton energy loss that is important for high pTpartons. 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 pTpartons 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 an 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 threequark 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]:

In the above,rBMis 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 rBM→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 rBM→∞,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.

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 yH 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 AMPTSM 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 towards one 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+/Kratio. 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 midrapidity 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 ˉ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 rBMvalue 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 rBMvalue 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.

Not only is the new quark coalescence able to describe the dN/dy yields,pTspectra,and elliptic flows of pions and kaons at low pT, but it also better describes the baryon observables in general, especially the pTspectra 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 rBMfactors 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 centerof-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:

In the above, ATrepresents the full transverse area of the overlap volume, and dET/dy is the initial rapidity density of the transverse energy at mid-rapidity, which is often approximated with the experimental dET/dy value in the final state.Because the Bjorken energy density diverges as t →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, we have calculated the initial energy density ∈(t) averaged over the transverse area of the overlap region as a function of time, including its maximum value ∈max[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 t1to t2:

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,t1=0.29dt&t2=0.71dtare used since they give the same mean and standard deviation of time as the beta profile at n=4.

We then considered both the finite time duration and longitudinal extension of the initial energy production[125]. When τFis not too much smaller than the crossing time of the two nuclei, results from this later study [125]are similar to those from the earlier study [124]. On the other hand, there is a qualitative difference in that the maximum energy density ∈maxat τF=0 is finite after considering both the finite duration time and longitudinal extension [125], while the Bjorken formula diverges as 1/τFand the method that only considered the finite duration time [124] diverges as ln(1/τF) at low energies but as 1/τFat 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 τFthan 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 z0and time t0of each excited string, which is then converted into the initial partons via string melting. Note that in the normal AMPT-SM model [13, 26, 27, 32], the longitudinal coordinate z0and time t0of 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. 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 effect of finite nuclear thickness on the energy density can be mostly neglected at high-enough energies.

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 in 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].

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 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) at τF =0.1 &0.3 fm/c are also shown for comparison

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

In the above, pTis the transverse momentum of the produced minijet parton,y1and y2are the rapidities of the two produced partons c and d,the factor K accounts for higherorder corrections to the leading order jet cross section, x1and x2are, respectively, the momentum fraction x carried by the two initial partons, fa(x1，Q2) is the parton distribution function (PDF) of parton type a at x=x1and factorization scale Q2, σabis the parton-parton cross section for partons a and b, and ^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 p0and considering all the possible combinations of final state parton flavors[23]:

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 implement 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 p0and σsoftparameters[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 of the incoming nuclei. Analogous to the free nucleon case, global analyses of 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 is 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. Note that the HIJING 1.0 shadowing is spatially dependent but independent of Q2or 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

where TA(s) represents the nuclear thickness function at transverse position s, and rAi(x，Q2，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.

As mentioned earlier, p0and σsoftare the two key parameters in the HIJING model that determine the total and inelastic cross sections of pp and pˉp collisions within the Eikonal formalism. In the HIJING 1.0 model that uses the Duke-Owens PDFs,constant values of p0=2.0 GeV/c and σsoft=57 mb are found to reasonably describe the experimental cross sections of pp and pˉ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 p0(s) and σsoft(s) values are needed.

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

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)

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, besides the update with modern parton distributions for proton and nuclei as discussed in Sect. 2.5,we have made several significant improvements on 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 twocomponent model. Second, we remove the minimum transverse momentum requirement (p0) 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]:

Fig.9 (Color online)Identified particle rapidity distributions in 0-5%central Pb+Pb collisions at =2.76 TeV from the AMPT-SM model using the minijet cutoff a (s) or b (s) in comparison with the ALICE data [99, 101]

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 →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],which leads to a smaller σjetthan Eq.(19)(at the same p0).As a result, an increase in the p0value as shown below is needed[122]for Eq.(19)to describe the experimental data on total and inelastic cross sections of pp and pˉp collisions shown in Fig. 8:

The total cˉ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 p0cut, 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 p0(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

Fig.11 (Color online)Charm quark a rapidity distributions and b pT spectra around mid-rapidity in pp collisions at=200 GeV(red curves)and 7 TeV (black curves) from the AMPT model. The shaded band represents statistical errors

Fig. 12 (Color online) dN/dy of charm quark pairs around midrapidity 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 rBMparameter as shown in Eq.(7). In our earlier work that updated the AMPT model with modern PDFs [119], the rBMvalue 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 rBMvalue 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 rBMvalue 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 rBMvalue 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].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 pT spectra of open charm hadrons in central Au+Au collisions at =200 GeV and b the Λc/D ratios versus pT for two centralities of Au+Au collisions at 200 GeV from the AMPT model in comparison with the STAR data [198, 201]

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 bLparameter 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 bLcan not describe the centrality dependence of〈pT〉 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 p0for 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 bLand p0parameters 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 selfconsistently.

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 aLand bLas [13]

Consequently,the〈pT〉of both 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 bLvalue 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 in bLwith the system size is consistent with the expectation of a stronger color field and thus a higher string tension κ since κ ∝1/bL[13] as shown in Eq.(6).

We propose that bLdepends on the local transverse position of the corresponding excited string inside the nucleus in each event [160]. Specifically, we assume that bLscales with the local nuclear thickness functions in a general AB collision as

where E0=200 GeV and Θ(x) is the unit step function.The fitted β(s) function is shown in Fig. 14a (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 bL∝1/TA(sA)for central AA collisions, where TA(sA) is proportional to the local number of participant nucleons or excited strings integrated over the longitudinal length.

Figure 14b shows the bLvalue 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〈bL〉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 〈bL〉 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〈bL〉is higher than that in Pb+Pb collisions at small b and grows faster with b due to its smaller system size.

Fig. 14 (Color online)a Fitted exponent functions α(s), β(s), and 3q(s)versus the center-of-mass energy per nucleon pairb Average p0 and bL values versus the impact parameter of pPb, Au+Au and Pb+Pb collisions at several energies

The minimum transverse momentum cutoff p0for 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 p0is 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 p0cutoff [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 p0as a local variable that depends on the transverse position of the corresponding hard process in each event. As p0is expected to increase with the system size, we related its value to the nuclear thickness functions in a general AB collision as [160]

Figure 14b also shows the average p0value 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 〈p0〉 decreases with the impact parameter and that〈p0〉at the lower LHC energy is smaller than that at the higher LHC energy due to the smaller α(s)value. Also, 〈p0〉 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〈p0〉with the impact parameter is seen to be much weaker than that of 〈bL〉 since α(s)?β(s) for the exponents in the local nuclear scaling relations.

Fig. 15 (Color online) a dNch/dη and b 〈pT〉 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 pT range used for the 〈pT〉calculation is [0.15, 2] GeV/c at 5.02A TeV and [0.2, 2] GeV/c at 200A GeV

The local nuclear scaling relations also predict how observables depend on the system size going from large to small systems. Figures 16a, b show, respectively, the dNch/dη and charged particle 〈pT〉 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 dNch/dη and 〈pT〉 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. In Fig. 16 we also see that the changes of the charge particle yield and〈pT〉from Cu+Cu to Au+Au collisions at 200A GeV are well accounted for by the local nuclear scaling. For example, the 〈pT〉 in Cu+Cu is generally smaller than that in Au+Au due to the larger bLvalue for Cu+Cu collisions. Note however that our calculations here have not considered the deformation of the Xe nucleus [227].

2.8 PYTHIA8 initial condition with sub-nucleon structure

Fig. 16 (Color online) a dNch/dη and b 〈pT〉 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 pT range used for the 〈pT〉 calculation is [0.15, 2]GeV/c at LHC energies and [0.2, 2] GeV/c at 200A GeV

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 〈pT〉 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 〈pT〉, in pp collisions. On the other hand, PYTHIA8 [228] is quite successful in describing the particle production 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 PYTHIA. With the proton charge distribution given by

where z is along the beam directions. In the second way,the initial transverse spatial condition including event-byevent 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), and then the transverse coordinates of the excited strings are randomly assigned to the binary collision center of each interacting constituent pair.

Fig.17 (Color online)〈pT〉of π(black),K(blue)and proton(red)at mid-rapidity within 0＜pT＜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 18a 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 18b 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.

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 subnucleon structure of the proton is implemented. We analyze the projected correlation function of two charged hadrons with a large pseudorapidity gap:

Both trigger and associate hadrons are required to be within 1＜pT＜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 Nsel, the number of selected charge tracks with pT＞0.4 GeV/c and ｜η｜＜2.4. High multiplicity events are defined as those with Nsel＞80, while low multiplicity events are defined as those with Nsel＜20.

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)

Figure 18c shows the multiplicity dependence of the C(Δφ) function from the two ways of generating the subnucleon 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 and 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 AMPTSM 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 above Eqs.(29-30)represent forward-angle scatterings.For isotropic scatterings, dσp/d^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].

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 pTdistribution, while panels (b) and (c) show the time evolution of parton 〈pT〉(scaled by T) and variance of pT(scaled by T2), respectively.The gluon system is initialized in a box with an offequilibrium 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. 19a 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 min(ct1，ct2)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 min(ct1，ct2) instead of(ct1+ct2)/2 enhances the collision rates and thus suppresses the causality violation.

We use the parton subdivision method to obtain the‘‘exact’’ time evolutions of 〈pT〉 and pTvariance (dashed curves) in Figs. 19b, c. We see that the time evolution of the pTvariance from the default scheme deviates significantly from the‘‘exact’’parton subdivision result,although the time evolutions of 〈pT〉 are close to each other (mostly due to the conservation of total momentum). In contrast,the time evolution of the pTvariance 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 →l×N，Vunchanged， (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 l2,the subdivision method is very expensive in terms of the computation time, which roughly scales with l2per 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:

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 106(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]

Fig.19 (Color online)a The final pT distribution,b time evolution of〈pT〉/T,and c time evolution of〈pT variance〉/T2 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=106 (dashed curves) for gluons in a box at T =0.5 GeV and σp =2.6 mb

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]

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:

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 for box calculations 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 is modified to specify the initial proton and neutron spatial distributions in the96Ru or96Zr 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 couples 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 K0or ˉK0.As a result,we change K0to K+and change ˉK0to 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 K0and change half of the final K-to ˉK0.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. Primary nucleon-nucleon collisions are then 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 nonequilibrium 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 stringdominated default version. The AMPT model should be improved to dynamically determine whether the QGP should be formed in the initial state; it would then selfconsistently 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 pTphysics once the radiative energy loss of high pTpartons 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.

Open AccessThis article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,adaptation,distribution and reproduction in any medium or format,as long as you give appropriate credit to the original author(s) and the source,provide a link to the Creative Commons licence,and indicate if changes were made.The images or other third party material in this article are included in the article’s Creative Commons licence,unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.