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

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.


I. 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][4][5][6][7][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][11][12][13][14], hydrodynamic models [15][16][17][18], and hybrid models [19][20][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 2 ) 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 * Z.-W.L. is supported in part by the National Science Foundation under Grant No. PHY-2012947. L.Z. is supported in part by the National Natural Science Foundation of China under Grant No. 11905188. † Corresponding author, linz@ecu.edu 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][31][32][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][36][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 arXiv:2110.02989v1 [nucl-th] 6 Oct 2021 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/ 3 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, we review in Sect. II 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. II A, 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. II B, the new quark coalescence model in Sect. II C, incorporation of the finite nuclear thickness along beam directions in Sect. II D, incorporation of modern parton distribution functions of nuclei in Sect. II E, improved treatment of heavy quark productions in Sect. II F, the introduction of local nuclear scaling of key input parameters to describe the system size dependence in Sect. II G, incorporation of PYTHIA8 and nucleon substructure in the initial condition in Sect. II H, and benchmark and improvement of the parton cascade algorithm in Sect. II I. We then briefly review other developments of the AMPT model in Sect. III. Finally, in Sect. IV, we summarize and discuss possible directions for further developments of the AMPT model.

II. 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. II B, II D, II E, II F, II G, II H are about the initial condition, Sect. II I is about the parton cascade, Sect. II C is about the hadronization, while Sect. II A 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. II A and II B; changes from the other developments will be released in the future.

A. 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][65][66][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][71][72][73][74], and dynamical models based on the kinetic theory [75][76][77]. We have modified the AMPT model to provide a kinetic theory description of deuteron production and annihilation by adding the following reactions [77]: where M = π, ρ, ω, and η, while B and B stand for baryons N , ∆, P 11 (1440), and S 11 (1535). Note that the hadron cascade component of the AMPT model [13], based on a relativistic transport (ART) model [84][85][86], already includes the interactions of π, K, η, ρ, ω, φ, K * , N , ∆(1232), P 11 (1440), S 11 (1535) as well as their antiparticles. For the cross sections of the reactions BB → 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 → dπ at the same center of mass energy √ s. The cross sections for the inverse reactions M d → BB 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 pp → dπ + [78][79][80] and the reaction π + d → pp [81][82][83]87] have been extensively measured, and the former is given by where p N and p π are, respectively, the magnitude of the threemomenta 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 where √ s is in the unit of GeV and f (s) is in the unit of mb. For the inverse reaction dπ + → pp, its cross section is related to that for pp → dπ + via the detailed balance relation: These parameterizations are compared with the experimental data in Fig. 1. The cross sections for the isospin averaged reactions N N → dπ and πd → N N can then be obtained as σ(N N → dπ) = 3σ(pp → dπ + )/4 and σ(dπ → N N ) = σ(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 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 BB ↔ M d 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.

B. 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] 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 [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 Fig. 2. 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, pT spectra dN/(2πpTdpTdy) in the unit of c 2 /GeV 2 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. 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. The AMPT-SM model can then 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 √ s NN = 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 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 [13,22] κ 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(−πm 2 ⊥ /κ) at transverse mass m ⊥ . As a result, the strange quark suppression relative to light quarks, exp[−π(m 2 s − m 2 u )/κ], 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][100][101] in panels (a) and (d), the p T spectra at mid-rapidity for the 0-5% centrality in panels (b) and (e), and v 2 {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 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 2 {EP} [102,103] at low p T .
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 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].
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 2 (η a , η b ) are shown in panels (c) to (f) for four centralities. We see that the AMPT-SM model reasonably reproduces the observed v 2 (η) 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 2 (η 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 v n will then develop longitudinal correlations.
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 parton energy loss 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].

C. 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 quarkantiquark 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 d M ) 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 d B ) 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]: otherwise : form a meson.
In the above, r BM 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 BM → 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 BM → ∞, 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 √ s NN = 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 poten-tial 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, presumably 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 dotdashed curve in Fig. 4 represents a function that is proportional to cosh y H , which qualitatively agrees with our model results. 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 r BM = 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 + /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.
On the other hand, the AMPT model with the improved quark coalescence [32,119] still underestimates thep/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 r BM 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 r BM 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.
Not only is the new quark coalescence able to describe the dN/dy yields, p T spectra, and elliptic flows of pions and kaons at low p T , but 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 r BM 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.

D. 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-ofmass 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 com-pletely cross each other in the center-of-mass frame: (8) in the hard sphere model of the nucleus. In the above, R A is the hard-sphere radius of the nucleus and y cm is the rapidity of the projectile nucleus. For central Au+Au collisions at √ s NN = 50 GeV, for example, d t ≈ 0.5 fm/c, which is comparable to the typical value of the parton formation time or hydrodynamics initial time when one takes R A = 1.12A 1/3 fm. Therefore, one may expect the effect from finite nuclear thickness to be significant for central Au+Au collisions at √ s NN 50GeV, 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: 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 → 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 t 1 to t 2 : In contrast, the beta time profile assumes the following: 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 1 = 0.29d t & t 2 = 0.71d t are 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 τ F is 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 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 0 and time t 0 of 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 z 0 and time t 0 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. 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 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 semianalytical 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.

E. 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 (12) In the above, p T is the transverse momentum of the produced minijet parton, y 1 and y 2 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 1 and x 2 are respectively the momentum fraction x carried by the two initial partons, f a (x 1 , Q 2 ) is the parton distribution function (PDF) of parton type a at x = x 1 and factorization scale Q 2 , σ ab is 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 p 0 and considering all the pos-sible combinations of final state parton flavors [23]: whereŝ is the standard Mandelstam variable for the minijet subprocess. We see that the minijet transverse momentum cutoff p 0 and the parton distribution functions f (x, Q 2 ) are the key factors affecting the jet cross section. The total jet cross section and the σ soft parameter that describes the soft component determine the nucleon-nucleon interaction cross sections in the Eikonal formalism [131,132]. Note that p 0 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 √ s NN < 10 GeV.
An important ingredient needed in Monte Carlo event generators for hadron collisions is the input parton distribution function [133][134][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 smallx 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 0 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 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][145][146][147][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][150][151][152][153][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 Q 2 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 where f p i (x, Q 2 ) is the corresponding PDF in the free proton. Here R A i (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 where T A (s) represents the nuclear thickness function at transverse position s, and r A i (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.   As mentioned earlier, p 0 and σ soft are the two key parameters in the HIJING model that determine the total and inelastic cross sections of pp and pp collisions within the Eikonal formalism. In the HIJING 1.0 model that uses the Duke-Owens PDFs, constant values of p 0 = 2.0 GeV/c and σ soft = 57 mb are found to reasonably describe the experimental cross sections of pp and pp 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 0 (s) and σ soft (s) values are needed. In our work that implements the CTEQ6.1M PDF in the AMPT model [119], the energy dependent parameters p 0 (s) and σ soft (s) are determined via fitting the experimental total and inelastic cross sections of pp and pp collisions within the energy range 4 < √ s < 10 5 GeV, as shown in Fig. 8. We then obtain the following p 0 (s) and σ soft (s) functions for pp collisions: In the above, p pp 0 and the center-of-mass colliding energy √ s are in the unit of GeV/c and GeV, respectively; while σ soft is in the unit of mb. Note that √ s will be replaced with √ s NN 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 pp collisions from √ s ∼ 4 GeV to 13 TeV [119].
When we apply the above p 0 (s), σ soft (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 pp 0 (s) value overestimates the final state hadron multiplicities in central Pb+Pb collisions at √ s NN = 2.76 TeV. Since a larger p 0 value would suppress the total jet cross section and reduce the particle yields, we introduce a global scaling of the minijet cutoff p 0 to make its value in central AA collisions, p AA 0 (s), nuclear size dependent [119]: 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 √ s NN ≤ 200 GeV since the p pp 0 (s) value works reasonably well there for central Au+Au collisions, while its value approaches 0.16 at √ s NN ∼ 10 7 GeV. This nuclear scaling of the minijet momentum cutoff scale p 0 is motivated by the physics of the color glass condensate [159], where the saturation momentum scale Q s depends on the nuclear size as Q s ∝ 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 AA 0 (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 very small value for the Lund b L parameter, b L = 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 0 and the Lund b L parameter with a local nuclear scaling [160], as shall be discussed in Sect. II G, which would help explain why a bigger p 0 value but a smaller b L value are needed for high energy AA collisions than pp 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. II E, 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 two-component model. Second, we remove the minimum transverse momentum requirement (p 0 ) 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]: 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 σ jet than Eq.(19) (at the same p 0 ). As a result, an increase (GeV) s   [191][192][193][194][195][196][197] as functions of the colliding energy. of the p 0 value as shown below is needed [122] for Eq. (19) to describe the experimental data on total and inelastic cross sections of pp and pp collisions shown in Fig. 8: with p 0 in GeV/c and √ s in GeV. The total cc 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 0 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 0 (dashed curve). Figure 11 shows the charm quark rapidity and transverse momentum distributions from the AMPT model for pp collisions at √ 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 0 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 0 cut (dashed curve) are similar, partly because a p 0 cut (2 GeV/c) on the charm quark production is also used in the original AMPT model.
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 √ s NN = 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. 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 pseudoscalar 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 in-clude 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 r BM parameter as shown in Eq. (7). In our earlier work that updated the AMPT model with modern PDFs [119], the r BM 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 r BM 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 r BM 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 r BM 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 * , D s 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 0 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 √ s NN = 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. 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 √ s NN = 200 GeV or below p T ∼ 15 GeV/c in Pb+Pb collisions at √ s NN = 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 → g + g cross section for scatterings of all parton flavors, where flavordependent 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. Fur-  thermore, hadronic scatterings of heavy flavor hadrons [202][203][204][205][206] have not been included in the AMPT model except for the decays of heavy hadron resonances. Future developments of the AMPT model are expected to improve its description of heavy flavor productions in AA collisions.

G. 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][208][209][210][211]. Second, we have found in earlier developments of the AMPT model [119,122] that the minijet transverse momentum cutoff p 0 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 0 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] .
Consequently, the p T of both partons after string melting and the final hadrons are significantly affected by the value of b L . Since the mean transverse momentum of initial partons in heavy ion collisions is expected to be higher in larger systems  [198,201].
due to the higher initial temperature, we expect the b L 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). We propose that b L depends on the local transverse position of the corresponding excited string inside the nucleus in each event [160]. Specifically, we assume that b L scales with the local nuclear thickness functions in a general AB collision as In the above, b pp L 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 ) = ρ A (s A , z)dz is the nuclear thickness function at the transverse distance s A from the center of nucleus A determined with the Woods-Saxon nuclear density profiles [156], and T p is the average value of the effective thickness function of the proton (taken as 0.22 fm −2 ). Note that T p is used instead of T A (s A ) or T B (s B ) when the projectile or the target is proton or when T A (s A ) or T B (s B ) 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 form (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 where E 0 = 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 ∝ 1/T A (s A ) for central AA collisions, where T A (s A ) is proportional to the local number of participant nucleons or excited strings integrated over the longitudinal length. 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 0 for light flavor minijet productions is another key parameter in the HI-JING 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 0 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 0 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 0 as a local variable that depends on the  transverse position of the corresponding hard process in each event. As p 0 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 pp 0 (s) As T A (s) ∝ A 1/3 , Eq. Since p pp 0 (s) works for charged particle yields in central Au+Au collisions at and below 200A GeV, we assume that the need to modify p 0 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] α(s) = 0.0918 ln with α(s) = 0 for √ 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 0 is closely related to the saturation momentum Q s in the color glass condensate [159], where Q s ∝ A 1/6 in the saturation regime. Figure 14(b) also shows the average p 0 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 0 decreases with the impact parameter and that p 0 at the lower LHC energy is smaller than that at the higher LHC energy due to the smaller α(s) value. Also, p 0 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 0 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][217][218][219][220][221][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 0 and b L but instead use constant b L = 0.15 GeV −2 and p 0 (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 AA 0 (s) and the dotted curves when using p 0 (s) = p pp 0 (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.
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 pp 0 (s), where the global nuclear scaling p 0 (s) = p AA 0 (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 0 in peripheral collisions is much smaller than that in central collisions, because the binary scaling of minijet productions makes p 0 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 AA  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][224][225][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. In Fig. 16 we also 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]. H. 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 HI-JING 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 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][231][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 ap- <3 GeV/c T 0<p π ALICE ALICE K ALICE p Fig. 17. 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.
proach, the fluctuating initial condition of AMPT originally provided by the HIJING model is replaced by the PYTHI-A/Angantyr model [229]. In addition, the sub-nucleon structure, which could be important for collectivity observables in small systems [234][235][236][237][238], can be modeled when implementing the space-time structure of the string system generated by PYTHIA. With the proton charge distribution given by 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) = ρ(x − b/2, y, z)ρ(x + b/2, y, z)dz, (27) 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][240][241][242][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 √ 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. 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.
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: 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 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][249][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.

I. 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 σ 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] where µ is a screening mass to regular the total cross section. This way the total cross section has no explicit dependence onŝ: The above Eqs. (29)(30) represent forward-angle scatterings. For isotropic scatterings, dσ p /dt 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][254][255]. These numerical artifacts due to the causality violation can be reduced or removed by the parton subdivision method [12,43,252,254,[256][257][258][259][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 σ 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; this collision time in the two-parton center of mass frame thus becomes two different collision times in the global frame (named here as ct 1 and ct 2 respectively for the two colliding partons) after the Lorentz transformation. The default collision scheme of ZPC [109] uses (ct 1 + ct 2 )/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 2 ), 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 min(ct 1 , ct 2 ) 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(ct 1 , ct 2 ) instead of (ct 1 + ct 2 )/2 enhances the collision rates and thus suppresses the causality violation.
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: 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 2 , the subdivision method is very expensive in terms of the computation time, which roughly scales with l 2 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: 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 6 (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 This is especially relevant for forward-angle scattering. For example, when parton subdivision requires the decrease of the forward-angle cross section of Eq.(30), one should not do that by increasing the screening mass µ by a factor of √ l because that would change the angular distribution of the scatterings in Eq. (29). Instead, one can decrease the α s parameter by a factor of √ 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][266][267][268][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 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.

III. 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 the 96 Ru or 96 Zr nucleus according to the density functional theory (DFT) calculations [271][272][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 K 0 orK 0 . As a result, we change K 0 to K + and changeK 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 0 and change half of the final K − toK 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 initialstate 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].

IV. 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 incorpora-tion 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 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.