Nuclear modification factor in intermediate-energy heavy-ion collisions

The transverse momentum dependent nuclear modification factors (NMF), namely $R_{CP}$, is investigated for protons produced in Au + Au at 1$A$ GeV within the framework of the isospin-dependent quantum molecular dynamics (IQMD) model. It is found that the radial collective motion during the expansion stage affects the NMF at low transverse momentum a lot. By fitting the transverse mass spectra of protons with the distribution function from the Blast-Wave model, the magnitude of radial flow can be extracted. After removing the contribution from radial flow, the $R_{CP}$ can be regarded as a thermal one and is found to keep unitary at transverse momentum lower than 0.6 GeV/c and enhance at higher transverse momentum, which can be attributed to Cronin effect.


I. INTRODUCTION
Recently, the nuclear modification factor (NMF) has been extensively investigated for different particles at various collision energies in relativistic heavy-ion collisions (HIC) [1][2][3][4]. These studies indicate that the NMF, which can be represented either by the modification factor between nucleus-nucleus (AA) collisions and protonproton (pp) collision R AA or the one between the central collisions and peripheral collisions R CP , is very useful for the study of the quantitative properties of the nuclear medium response when the high speed jet transverses it. In high transverse momentum (p T ) region, NMF is suppressed owing to jet quenching effect in hot-dense matter and thus has become one of the robust evidences on the existence of the Quark-Gluon-Plasma [5,6]. In lower p T region, radial flow boosts or the Cronin Effect [7] competes with the quenching effect and enhances the NMF, which has been also demonstrated by the Relativistic Heavy-Ion Collider (RHIC) beam energy scan (BES) project [8].
Meanwhile, collective motion plays an important role in the time evolution of particles, which has been studied over a wide range of collision energy in heavy-ion collision. Around 1A GeV incident energy in central HIC, the colliding nuclei are expected to be stopped and lead to densities of 2 ∼ 3ρ 0 (ρ 0 is the normal nuclei density) at the largest compression time [9]. At this point, high excitation energy stage is reached and some parts of the excitation energy are converted into collective motion, such as radial flow [10,11]. Thus, in the following expansion stage, the products move outward containing both the collective motion and the thermal motion. It will be of highly interesting to decouple these two parts, * Author to whom all correspondence should be addressed. Email: ygma@sinap.ac.cn because each of them reflects important information of the HIC process. Efforts have been made by the Blast-Wave model [12][13][14], and the collective motion parameter (radial flow velocity) together with the thermal motion parameter (temperature) can be extracted at the same time.
Nevertheless, till now, NMF has not been investigated in intermediate energy HIC yet to our knowledge. In particular, the R CP shape will be strongly affected by the radial flow which plays a significant role in the expansion stage especially for the low p T particles. In order to understand the properties of the nuclear medium response, one may want to know what the behavior of R CP without the contribution from radial flow might be. In the present paper, we address NMF in intermediate energy HIC to study the properties of nuclear medium for the first time. After removing the contribution from the radial flow, the NMF can be regarded as a thermal one, which reflects property of thermal medium produced in intermediate energy HIC. The R CP for emitted protons in Au + Au collision at 1A GeV is investigated systematically.
The article is organized as follows: A brief introduction about the IQMD model is given in Sec. II. In Sec.III, we compare the kinetic energy spectra of light fragments obtained by the IQMD with the one from the EOS experimental result, then we fit the proton transverse mass (m T ) spectra with the Boltzmann function and with the distribution function from Blast-Wave model. In Sec. IV, we show the results of NMFs and compare them with the results extracted from the KaoS experimental data. Finally, the thermal R CP is recalculated by removing the contribution from radial flow. Summary and conclusion are presented in Sec. V.

II. BRIEF DESCRIPTION OF IQMD MODEL
The Quantum Molecular Dynamics model is a transport model which is based on a many body theory to describe heavy ion collisions from low (dozens of MeV) to relativistic energy [15][16][17][18]. The Isospin dependent Quantum Molecular model (IQMD), was extended from the QMD model, with considering the isospin effects. In the past decades, many applications have been successfully performed into nuclear physics studies with the help of IQMD. For instance, IQMD has been successfully applied to treat collective flow, multi-fragmentation, isospin effects in HIC, transport coefficient in HIC, giant resonance, and strangeness production etc [19][20][21][22][23]. In the IQMD model, the wave function of each nucleon is described as a coherent state with the form of Gaussian wave packet, (1) Here r i and p i are the time dependent variables which describe the center of the packet in coordinate and momentum space, respectively. The parameter L, related to the width of wave packet in coordinate space, is determined by the size of reaction system, i.e. L = 1.08 fm 2 for Ca+Ca system and L = 2.16 fm 2 for Au+Au system in this work. The wave function of the system is the direct product of all the nucleon wave functions without considering the Fermion property of nucleon: (2) As a compensation, Pauli blocking is employed in the initializations and collision process to restore some parts of the quantum property of many Fermion system. By applying a generalized variational principle on the action of the many-body system, one can get the equations of motion for p i and r i , which are listed as follows The Hamiltonian H = T + V where T is the kinetic energy, the potential V is expressed by (4) In the above, the Wigner distribution function f i ( r, p, t), which is the phase-space density of the ith nucleon, is obtained by applying the Wigner transformation on the single nucleon wave function: The baryon-potential consists of the real part of the G-Matrix which is supplemented by the Coulomb interaction between the charged particles. The former one can be divided into three parts, the Skyrme-type interaction, the finite-range Yukawa potential, and the momentumdependent interaction (MDI) parts. The two-body interaction potential V ij in Eq. 4 can be expressed as follows: The symmetry potential between protons and neutrons corresponding to the Bethe-Weizsacker mass formula can be taken as with t 6 =100 MeV. By integrating Skyrme part as well as the momentum dependent part of the two-body interaction and introducing the interaction density, one can get the local mean field potential which contains the Skyrme potential and momentum dependent potential where ρ 0 is the saturation density at ground state, (h) 2 is the momentum distribution function, the interaction density ρ = ij ρ ij , and α, β, and γ are the Skyrme parameters, which connect tightly with the EOS of the bulk nuclear matter, as listed in Table I.  [15]. With the help of coalescence mechanism, the information of fragments produced in HICs can be identified in IQMD. A simple coalescence rule to form a fragment is used with the criteria ∆r = 3.5 fm and ∆p = 300 MeV/c between two considered nucleons.

III. TRANSVERSE MASS SPECTRA
The Au+Au collisions at 1A GeV are simulated with the IQMD model for both the soft equation of state with MDI (SM) and the hard equation of state with MDI (HM). In order to test the reliability of the model, the kinetic spectra of proton, deuteron and triton obtained by the IQMD with the above equation of state situations are compared with the EOS experimental data in Fig 1 [12]. The conditions of the chosen fragments in our model calculations are b ≤ 3f m, θ cms = 90 ± 15 • , which are the same as the EOS experimental situation. As shown in Fig 1, the solid lines are the spectra extracted from the IQMD with soft and MDI potential, and the dash lines are the spectra extracted with HM potential. From Fig. 1, we find that both the soft and hard equation of state situation can essentially describe the EOS experimental data nicely. A little higher yield of protons and lower yield of t are observed in the case of the HM potential. The enhanced yield of protons can be attributed to stronger mean field interaction in HM case compared to that in SM case. The solid symbols are the experimental data from the EOS collaboration [12], and the solid line is our IQMD simulation with the soft+MDI potential, and the dash line is for the hard+MDI potential.
Transverse mass m T , is given by m T = p 2 T + m 2 0 , with the rest mass of protons m 0 . Fig. 2 shows our simulation results for the double differential transverse mass spectra ( 2πmT dmT dy ) of protons at different centralities (0−10%, 10−20%, 20−40%, 40−60%, 60−80%) together with two fit results either from Boltzmann distribution function or from Blast-Wave model distribution function in the framework of IQMD model with SM potential. The c.m.s. rapidity cut is |Y /Y proj | < 0.1 where Y proj is the initial projectile rapidity. At the end of this section, a discussion will be given for the fit parameters, i.e. temperature and radial flow.

A. The Boltzmann fit
In heavy-ion collisions, particles collide with each other randomly, which can be described in term of thermal motion [24]. Lots of works demonstrate that the Boltzmann distribution can roughly describe the particle spectra by a thermal source in HICs [25][26][27][28]. If the particles are ejected from a single pure high temperature thermal source, its transverse mass spectra should satisfy the Boltzmann distribution function, i.e. dN/(dm T dy) ≈ f (m T ) · exp(−m T /T ), which is close to a straight line when the y-axis is plotted logarithmically. The results of the m T spectra of protons fitted with the Boltzmann function are shown in Fig. 2(a). The proton spectra at different centralities are fitted by the Boltzmann function in the region of m T >0.165 GeV and extrapolated to the low m T region. It is interesting to find that the spectra in central collisions have an obvious "shoulder"like structure in the low m T region while there is no such structure for peripheral collisions. This means the spectra in central collisions deviate from the Boltzmann distribution in the low m T region, which can be essentially attributed to the radial flow in central collisions. The fit parameters of the slope temperature are listed in the second column of Table II. From peripheral to central collision the slope temperature extracted from the Boltzmann function fit shows a decreasing trend, which is very different from that extracted from the Blast-Wave model fit as discussed in next section. Considering that the Boltzmann function cannot fit well the m T spectra in the low m T region due to the existence of collective radial flow, we adopt a hydrodynamically inspired "Blast-wave" model to describe the m T spectra with two free parameters: collective transverse flow velocity β and kinetic freeze-out temperature T f . In this model, the spectrum is computed by boosting the thermal source both in longitudinal and transverse direction [10,29]. The radial velocity distribution β r in the region 0 ≤ β r ≤ R is described by a self-similar profile, which is parameterized by the surface velocity β S : where R is the freeze-out radius, namely the maximum radius of the expanding source at thermal freeze-out time. β S is the particles' radial velocity on the surface of the freeze out volume when r = R, and the exponent α represents for the radial flow profile which describes the evolution of the radial flow velocity (when α =0, it means uniform velocity; when α =1, the expansion is similar to Hubble's law; when α = 2, it corresponds to a hydromechanical expansion). The particle spectra are the superposition of individual thermal source at different r, each boosted with the angle ρ = tanh −1 β r (r): where K 1 , I 0 are the modified Bessel function, T f is the freeze-out temperature. The spectrum shape is determined by the freeze-out temperature T f , the velocity of the transverse expansion β S , the flow velocity profile α and the mass of the particle m 0 . The average flow velocity is estimated by taking an average over the transverse geometry: β r = β S 2 2+α . The spectra of protons can be well described by the Blast-Wave model as shown in Fig. 2(b) and the fit parameters are listed in Table III. Obviously, the whole spectra are well fitted, even for the "shoulder structure" at low m T in central collisions (for the centralities of 0-10% and 10-20%). From central collisions to peripheral ones, the velocity β S decreases and while the temperature has only a slight increase [30]. It is noted that the similar trend was also seen in the fitting result in RHIC energy region, which might be owing to non-equilibrium effect in peripheral collision. For example, by replacing the Boltzmann statistics with the Tsallis statistics in traditional Blast-Wave model, the Tsallis Blast-Wave model (TBW) can describe the peripheral spectra well and give a more reasonable fitting result (a low value of β S and (q − 1) > 0 at peripheral collisions) [31].
The main results from the fitting of the m T spectra are following: 1. Radial flow ( β r ). In Fig. 2(b), the proton spectra have been described very well with the Blast-Wave model. The average radial flow is about 0.33c in the central collisions (0-10%), and reduces to 0.22c in peripheral collisions (60-80%). This result with a larger radial flow existing in central collisions is consistent with the physics picture of HICs.
In another scenario, the average kinetic energy E kin of protons, deuterons and tritons at mid-rapidity are consistent with the picture described by the Blast-Wave model. Fig. 3 shows the mass number dependence of E kin . They have a linear relationship: E kin = E th + E f low = η · A + 3T /2 where η = 1/ 1 − β 2 r − 1. The slope is the quantity related to radial flow, and the intercept is the quantity related to thermal temperature. We can also extract the radial flow and temperature information by the linear curve about E kin vs A, the result is listed in the lower part of the Table III, which is comparable with the Blast-Wave fit result.
2. Freeze out temperature (T f ) We use Boltzmann distribution to fit the proton spectra, with the fitting parameters at different centralities called "slope temperature" listed in Table II. The temperature shows a decreasing trend from 93.06 MeV in central collisions to 73.87 MeV in peripheral collisions. Comparing with the Boltzmann fit results, the freezeout temperature from the Blast-Wave fit is smaller, and shows a little increasing trend. This difference is owing to that the collective radial flow motion effect has been misunderstood as the thermal part of motion in the Boltzmann description.
3. Flow profile (α) The α exponent describes the evolution of the radial  Table III . The parameters of radial flow βS, βr and the freeze-out temperature T f at various centralities with different fitting methods. The upper one is from proton mT spectra described by the Blast-Wave model. The lower one is from fitting the average kinetic energy E kin of p, d, t, which is displayed in Fig. 3. flow velocity from any radius r (0 < r < R). We set the parameter α free in Blast-Wave fit, which is different from some previous works keeping a fixed value equal to 1 or 2 [29,31,32] or some assuming a common flow (i.e. α = 0) [11,12,33,34]. The α makes a little difference to the fitting value of β S while the choice of R value have no influence on the fitting which chosen R = 40 fm in our work. We get the α value equal to 0.336 on the spectra of 0-10% centrality and fix this value when do the fitting on the other centralities. A similar value of α has been also obtained in Ref. [35] for Au + Au collisions in ultrarelativistic energy.

IV. NUCLEAR MODIFICATION FACTOR
To explore the nuclear medium response, the ratio R CP between the particle yield in central collisions and the particle yield in peripheral collisions, has been introduced. Both yields are normalized by corresponding nucleon-nucleon binary collision numbers N coll (binary scaling): where the yield is the differential invariant yield ( d 2 N 2πpT dpT dy ). If the nucleus-nucleus collision is a mere superposition of N coll independent nucleon-nucleon collisions, R CP would be unit with p T . Thus any departures from R CP =1 indicate nuclear medium effects or other dynamical effects.

A. The normal NMF
The normal R CP versus p T is obtained by dividing the p T spectra in central collisions with the ones in peripheral collisions with taking the respective binary collision number into account. To test the reliability of model, the results from IQMD calculation are compared with that extracted from KaoS experimental data. Following the experimental conditions [26], the central collision events are chosen with charged particle multiplicity M z > 50, and the peripheral collision events with M z < 26. The identified protons are chosen the same as the KaoS experiment, namely, in the polar angle of 40-48 degrees, and the momentum range 0.320-1.440 GeV/c. The spectra of proton in central and peripheral collisions are showed in Fig. 4(a), together with the result from IQMD+SM(HM) representing with different type lines elaborated in figure caption. Results demonstrate that both the soft and hard EOS can describe the experimental data nicely. Shown in Fig. 4(b) is the R CP as a function of total laboratory momentum from IQMD simulation and the KaoS experiment. Two parameter sets of equation of state, namely SM and HM, are employed in IQMD calculation. The R CP of protons from experiment (circle) and from IQMD calculation (solid red line: SM; dash blue line: HM) increase almost linearly with the total momentum in laboratory system (p lab ). This is consistent with those behaviors from lower RHIC Beam-Energy-Scan energies, e.g. the R CP for proton in Au+Au at 7.7A GeV [8], and can be understood by radial boosts and/or the Cronin Effect [7]. The R CP of proton is not so sensitive to the hard or soft equation of state in the IQMD calculations, indicating the multiple collision process dominates here. It's a need to note that the Rcp from model are arbitrarily scaled to match with the experimental result owing to the unknown N coll in KaoS experiment.
In IQMD model, two particles collide if the minimum relative distance of their centroids of the Gaussian wave during their motion, in the CM frame fulfills the requirement: where the cross section σ tot is assumed to be the free cross section of the regarded collision type (N-N, N − ∆, etc.). In this work, the collision numbers in every centralities are counted with the "N collt " counter when each collision occurs. For each collision the phase space densities in the final states are checked in order to assure that the final distribution in phase space is in agreement with the Pauli principle, the Pauli blocking number has been counted by the "N paubl " counter. Then, we can get the real collision number (N coll = N collt − N paubl ). For Au+Au collision at incident energy 1A GeV, the average collision numbers are 242. 5, 473, 795.5, 1136, 1463 in five centralities from peripheral to central collisions in our IQMD model, respectively. The R CP versus p T shows centrality dependence. In Fig. 5(a), four cases at different centralities are investigated, which are 0−10% 60−80% (red dot), 10−20% 60−80% (green square), 20−40% 60−80% (blue up triangle) and 40−60% 60−80% (pink downtriangle). From the most peripheral case in numerator (i.e. 40−60% 60−80% ) to the most central case in numerator (i.e. 0−10% 60−80% ), the R CP becomes more and more p T dependent. For all these cases, R CP rises with p T , with a cross point shows up at p T =0.5 GeV/c, which may suggest a balance between two mechanisms for the p T dependence of R CP .
On the one hand, in the low p T region (p T less than 0.5 GeV/c), radial flow takes a major role in central collisions, which pushes protons to higher p T region and results in the smaller R CP at low p T . On the other hand, in the high p T region (p T greater than 0.5 GeV/c), the Cronin effect (nucleon multiple scattering effect) [36,37] tends to transform the longitudinal momentum into the transverse momentum and increase the p T in central HIC, leading to the larger R CP at high p T .

B. The thermal NMF
In order to investigate the thermal property of the medium in the collision overlapping region, one needs to deduct the contribution of collective flow to focus on thermal effect. The radial flow parameter extracted by the Blast-Wave fit has been displayed in Table III. By subtracting the radial velocity of each particle, we recalculate the R CP of protons without the radial flow contribution.
To determine the magnitude of radial flow, one needs to figure out the size of freeze-out volume. In the transport model IQMD, a density criterion is applied here. During the expansion process, the stage can be chosen as freeze-out stage, when the density at the center of each collision reaches 1/8ρ 0 . At this stage, all the products (including fragments and free nucleons) cease the strong interaction among them and almost fix their momenta. For the case of Au+Au collision at E beam = 1A GeV, this freeze-out time is about 60 fm/c, and the corresponding radius R max of freeze-out volume is about 40 fm (half of the reducing edge of radius distribution of all emitting protons), which is of course model dependent. Fig. 5(b) shows the NMF of protons inside the freeze-out volume. In comparison with Fig. 5(a), freeze-out sphere cuts off most of higher p T protons. Once the freeze-out sphere is fixed, the radial flow β r for the particles inside the volume can be calculated through formula in Eq. (10), here α = 0.336 is taken from the Blast-Wave fit for the p T spectra. The contribution of radial flow on p T is the projection of radial momentum (p r = m 0 · β r · γ r ) on the p T vector, and the thermal transverse momentum is p T = p T − pr * pT | pT | . The thermal R CP , is then obtained by dividing the thermal p T distribution in central collisions to that in the peripheral collisions with taking the respective binary collision number into account, which is shown in the Fig. 5(c). It is found that the thermal R CP becomes almost unitary below 0.6 GeV/c, which illustrates that the original increasing R CP behavior in low p T region essentially stems from the collective radial flow effect, and the thermal motion of nucleus-nucleus collision can be seen as the independent overlap of nucleonnucleon collisions. It is also noted that the thermal R CP is larger than 1 at p T above 0.6 GeV/c, where the Cronin effect plays a dominant role for the R CP increasing behavior versus p T [7,36].

V. CONCLUSION
In summary, the nuclear modification factor has been introduced in intermediate energy HIC in the framework of transport model, and thermal NMF has been proposed. First, we test the reliability of the IQMD model by comparing the kinetic energy spectra of light fragments (p, d, t) obtained by the IQMD model with the EOS experimental data, and we found the model can describe the data nicely. Secondly, the transverse mass spectra of protons are fitted by the Boltzmann distribution function as well as the distribution function from Blast-Wave model. It is found that the latter gives a satisfying description of the transverse mass m T spectra of protons at freeze-out time, which demonstrates that the kinetic energy of protons contains the collective radial motion part and random thermal motion part. The fitting radial flow parameter β S and temperature parameter T f are consistent with the results from the EOS collaboration [12] and the FOPI collaboration [13,14]. While, the Boltzmann description overestimates the proton yield at low m T , because of the missing radial flow contribution. Thus the fit temperature from the Boltzmann description can be interpreted as the apparent temperature of the emission source, which decreases from the central collisions to the peripheral collisions. It is found that both the radial flow effect and Cronin effect play their corresponding roles in shaping the p T dependence of R CP . The radial flow magnitude can be extracted by fitting the m T spectra with the distribution function from Blast-Wave model. The thermal modification factor is then obtained by removing the contribution from the radial flow. It is found that the thermal R CP of protons is close to 1 at lower p T , where protons emitted from Au+Au collisions can be regarded as independent superposition of emitted protons by nucleon-nucleon collisions, and while R CP enhances at higher p T , where the Cronin effect, i.e. multiple production mechanism of protons plays an increasing role. On the other hand, in light of this study, the nuclear modifi-cation factor at very high p T which has been often used to investigate jet-quenching effect at RHIC and LHC energies was actually nearly not affected by the collective radial flow.