Quantifying the Effect of Initial Fluctuations on Isospin-Sensitive Observables from Heavy-Ion Collisions at Intermediate Energies

Initial fluctuation is one of the ingredients that washes fingerprints of the nuclear symmetry energy on observables in heavy-ion collisions. By artificially using the same initial nuclei in all collision events, the effect of the initial fluctuation on isospin-sensitive observables, e.g., the yield ratio of free neutrons with respect to protons Nn/Np, 3H/3He yield ratio, the yield ratio between charged pions π−/π+, and the elliptic flow ratio or difference between free neutrons and protons v2n/v2p (v2n-v2p), are studied within the ultrarelativistic quantum molecular dynamics (UrQMD) model. In practice, Au + Au collisions with impact parameter b = 5 fm and beam energy Elab = 400 MeV/nucleon are calculated. It is found that the effect of the initialization on the yields of free protons and neutrons is small, while for the yield of pions, the directed and elliptic flows are found to be apparently influenced by the choice of initialization because of the strong memory effects. Regarding the isospin-sensitive observables, the effect of the initialization on Nn/Np and 3H/3He is negligible, while π−/π+ and v2n/v2p (v2n-v2p) display a distinct difference among different initializations. The fingerprints of symmetry energy on π−/π+ and v2n/v2p can be either enhanced or reduced when different initializations are utilized.


Introduction
Nuclear symmetry energy, which describes the energy difference between pure neutron matter and isospin symmetric (with equal numbers of protons and neutrons) nuclear matter, is one of the crucial quantities for studying the structures and the properties of nuclei and neutron stars, the dynamics of heavy-ion collision, supernovae explosions, as well as neutron star mergers [1][2][3][4][5][6][7][8][9][10][11][12][13]. Exploring nuclear symmetry energy at various densities is one of the important scientific goals for intermediate-energy heavy-ion collision (HIC) studies in terrestrial laboratories [14][15][16][17][18][19][20][21][22][23]. Quite a few observables have been found or predicted to be sensitive to the nuclear symmetry energy as, e.g., neutron and proton yields and flow ratios, double ratios, or differences, 3 H/ 3 He yield ratio, π − /π + and K 0 /K + meson production ratios, the Σ − /Σ + ratio [24][25][26][27][28][29][30][31][32][33][34][35][36]. In spite of the progress made, a tight constraint on the density-dependent nuclear symmetry energy E sym (ρ) by using HIC observables is still very difficult to achieve [37][38][39][40][41]. Extracting E sym (ρ) with HIC should rely on both transport model simulations and experimental measurements, while various transport models have different philosophies/assumptions/parameters, and as a result constraints on E sym (ρ) with different transport models are usually different to some extent. In view of this situation, a Transport Model Evaluation Project (TMEP) which intends to understand the source of the discrepancies and evaluate transport models was begun several years ago [42][43][44][45]. In addition, the initial fluctuations (the preparation of initial nuclei, i.e., target and projectile) and dynamical fluctuations (i.e., stochastic nucleon collisions) might wash the effects of E sym (ρ) on observables, which make the constraint of E sym (ρ) challenging. The aim of this work is to quantify the effect of initial fluctuation on several observables, especially its influence on isospin-sensitive observables.
The paper is arranged as follows. In the next section the treatments of the initialization process of the ultrarelativistic quantum molecular dynamics (UrQMD) model is presented. In Section 3, the influence of different initializations on various observables is shown and discussed. Finally, a summary and outlook are given in Section 4.

Initialization Process
In the default UrQMD model [46,47], there are two different choices for the initialization, hard-sphere and Wood-Saxon. According to the philosophy of the UrQMD model, each nucleon is respected by a Gaussian with a finite width σ r . Usually σ r depends on the mass number of target (projectile), σ r = 1.44 fm is used for collision with Au. By considering the width effect, hard-sphere is an economical way to sample nucleons in the coordinate space, thus hard-sphere is used in the present work. The initialization is a two-step process in the default UrQMD model. Firstly, calculating the radius of the target (projectile) according to an empirical formula [46] Here, ρ 0 = 0.16 fm −3 is the saturation density, A is the mass number. We note here that the radius calculated with this formula is smaller than the usually used one R = 1.12 * A 1/3 fm. The centroids of the Gaussians (nucleons) are randomly distributed within a sphere of the radius R, but if the minimum distance between two nucleons in a sampled nucleus is smaller than d min = 1.6 fm, it will be re-sampled. d min is set to 1.6 fm in consideration of the size of the nucleon. Indeed, d min is a free parameter of QMD-like models, different values are used in different codes [42][43][44][45]. This minimum distance limitation is useful to prevent from getting unphysical nucleus. Secondly, the momenta of all nucleons in the projectile or target are randomly sampled between 0 and the local Fermi momentum p F . In the default version of UrQMD, p Fp =h(3π 2 ρ p ) 1/3 for protons and p Fn =h(3π 2 ρ n ) 1/3 for neutrons are used. To more closely mimic the real nucleus, the binding energy per nucleon B of the sampled nucleus will be calculated, the phase space of each nucleon might be re-sampled until the absolute difference between the calculated B and the corresponding experimental value is smaller than 1 MeV. The above treatment is named as Normal-IN in the present work.
It is known, that the Wood-Saxon density distribution is a widely used assumption for nuclear density distribution in a nucleus, which reads, Here a = 0.54 and R = 6.52 fm are used for Au. We propose a density constraint method in sampling the coordinate of each nucleon. After the coordinates of each nucleons are sampled, the density at r = 0, 4, 5, 6, 6.5, and 7 fm will be calculated, if the relative difference between the calculated density and the one calculated with Equation (2) is larger than 10% (for r = 0, 4, or 5 fm) or 25% (for r = 6, 6.5, or 7 fm), the coordinate of all nucleons will be re-sampled. With this treatment, fluctuation in coordinate space will be partly reduced. We name this method as DC in the present work. Alternatively, by neglecting the isospin effect on Fermi momentum, one can use p F =h(3π 2 ρ 2 ) 1/3 for all nucleons instead of p Fp and p Fn , which is named as DCPF in the present work. In this mode, fluctuations in the momentum space can be reduced to some extent.
We randomly select nine initializations obtained from the Normal-IN mode, and the nucleon distribution in coordinate and momentum spaces are displayed in Figure 1. Fluctuations in both coordinate and momentum spaces are very large. Relatively, fluctuations in coordinates are smaller because of the minimum distance limitation in the sampling. Figure 2 displays the momentum distribution and its standard deviation obtained with Normal-IN and DCPF modes. As can be seen, the standard deviation in DCPF mode is slightly smaller than that in Normal-IN, which means smaller fluctuations in the momentum space. The reduction on the standard deviation in large momentum region is more visible because samples with larger densities are resampled in the DCPF method. This may have some effect on pion production as initial nucleons with higher momentum have a higher probability to be converted to ∆ resonances so as to produce pions from their decay.

Influence on Time Evolution
In the present work, we artificially use the same initial nuclei for all collision events. Calculations with the nine initializations are performed. In each case, more than 300,000 Au + Au collision events with impact parameter b = 5 fm and beam energy E lab = 400 MeV/nucleon are simulated to ensure enough statistics. In addition, calculations with the Normal-IN mode, as well as DC and DCPF modes are also performed. In the presently used UrQMD model, the symmetry potential is derived from the Skyrme potential energy density function in the same manner as the improved quantum molecular dynamics (ImQMD) model, see e.g., [48,49]. Together with the introduction of the in-medium nucleon-nucleon cross-sections and the improvement on the cluster recognition criteria, many experimental data in HICs at intermediate energies can be reproduced fairly well [50][51][52]. In the presently used UrQMD version, at the end of the reaction, clusters are recognized with an isospin-dependent minimum spanning tree (iso-MST) method [18,53]. In this method, if the relative distances and momenta between two neutrons (protons) or between one neutron and one proton are smaller 3.8 fm (2.8 fm) and 0.25 Gev/c, respectively, they are considered to belong to the same cluster. Throughout, nucleons that are not bound in clusters are called free nucleons. To quantitatively compare the influences of the initialization and the nuclear symmetry energy, results calculated with two Skyrme interactions (Skz4 and SkI1) which yield very soft (with the slope parameter L = 5.8 MeV) and stiff (L = 159 MeV) symmetry energies are also presented, respectively. Overall, 24 calculations with different initializations and symmetry energies are performed. Figure 3 illustrates that even with the same initialization, the result of each collision event still displays distinct fluctuations because of the following stochastic two-particle collisions. Particle distributions among different events become different at about t = 5 fm/c because the target and projectile start to touch each other at that time. At about t = 20 fm/c, the maximum compression is achieved. By artificially using the same initial nuclei in every collision event, one can get the result of each observable from each event and the averaged one over all events. The averaged nuclear density distribution in the reaction plane at t = 0 and 20 fm/c are displayed in Figure 4. The small difference in density distribution at an initial time may result in different evolution afterward. For example, the maximum density reached in calculations with SameIN3 is about 0.35 fm −3 , while it is about 0.32 fm −3 for SameIN9. Moreover, the shapes of the compressed region also display visible differences among different initializations, indicating that the information of the initial nuclei are memorized during the collision. The strong memory effect is one of the features of HICs at intermediate energies [54,55]. We have checked that not only the total densities are different in calculations with different initializations, the neutron and proton densities also display large differences in different calculations. As a result, one may expect that the isospin-sensitive observables will be influenced to some extent. It is known that both the shape and the density distribution of the compressed region are closely related to the elliptic flow. It is easy to imagine that the elliptic flow will be different with different initializations. In addition, ∆ resonances, which are parents of the produced pions, are created in the compressed region through the inelastic scattering of two particles [56,57] so that the pion related observables should be also influenced to some extent by the initialization.

Influence on Rapidity Distribution
The rapidity distributions of the yield, the directed flow v 1 , and the elliptic flow v 2 of free protons are displayed in Figures 5-7, respectively. In the present work, v 1 and v 2 are deduced from the Fourier expansion of the azimuthal distribution of detected particles, in which p x and p y are the two components of the transverse momentum p t = p 2 x + p 2 y . Additionally, the angle brackets in Equations (3) and (4) indicate an average over all considered particles from all events. For analysis flow in heavy-ion collisions at high energies, one should take care of non-flow effects caused by resonance decays and jet production and so on, while it is negligible for intermediate energy heavy-ion collisions. It is known, that for mass symmetric collisions, i.e., the target and projectile nuclei are the same, the yield and v 2 as a function of rapidity y 0 exhibit even functions, and v 1 being an odd function. While the results obtained with the nine initializations are more or less asymmetric, and large difference appears around target/projectile rapidities (y 0 = ±1) as the nucleons in these region are more likely to inherit more information of the initialization. Further, in Figure 7, v 2 obtained with different initializations is not a smooth function of the rapidity. This cannot be fixed with increasing the number of collision events, as the statistical error has been already smaller than the symbol size. The un-smooth behavior of v 2 around target/projectile rapidities can be understood from its definition Equation (4), which enlarges the effect of anisotropic momentum distributions from the initialization. In addition, it can be seen that the results obtained with Normal-IN, DC, and DCPF are almost overlapping. Only a tiny difference can be found on v 2 around target/projectile rapidities, implying v 2 around y 0 ≈ ±1 is a good candidate for probing the initial state of nuclei.     To quantitatively manifest the effect of initialization on these observables, the yields of free protons, free neutrons, produced π − and π + , as well as the v 1 slope at midrapidity and v 2 of free protons within rapidity window |y 0 | < 0.5 calculated using the soft symmetry energy (Skz4) are listed in the left side of Table 1. The influence on the yield of free protons and neutrons is rather small because it is mainly affected by the selection of the nuclear mean-field potential. However, the others display significant differences with different initializations, as they are more strongly related to the status of the compressed region. For example, the yield of π − for SameIN8 is 0.994 while it is 1.217 with SameIN3. v 2 obtained with SameIN2 and SameIN9 are −0.034 and −0.027, which are the two weakest elliptic flows among all calculations as well, respectively. While v 2 obtained with SameIN1 and SameIN3 are −0.049 and −0.045, respectively, which are the two strongest elliptic flows among all calculations. It can be partly understood from the density distribution in the compressed region, as shown in Figure 4, the maximum densities in SameIN2 and SameIN9 cases are the smallest while they are the largest in SameIN1 and SameIN3. Table 1. Column 2-7: The total yields of free neutrons, free protons, π − and π + , as well as the v 1 slope at midrapidity and v 2 of free protons within rapidity window |y 0 | < 0.5, as calculated with the soft symmetry energy. Column 8-15: The total yield ratios N n /N p , the elliptic flow ratio v n 2 /v p 2 within rapidity window |y 0 | < 0.5, π − /π + , and 3 H/ 3 He. Results calculated with both the soft and stiff symmetry energies are displayed. Statistical error is small and not shown here.

Influence on Isospin-Sensitive Observables
As usually the effect of nuclear symmetry energy on various observables is relatively small, the yield or flow of particles with different isospin are presumed to be sensitive to the nuclear symmetry energy. For example, the yield ratio of free neutrons respect to protons N n /N p , pion yield ratio π − /π + , the elliptic flow ratio (difference) between free neutrons and protons v n , the yield ratio of 3 H with respect to 3 He, have been used as sensitive probes to constrain the nuclear symmetry energy. To quantitatively compare the influences of initialization and the nuclear symmetry energy, N n /N p , v n 2 /v p 2 , π − /π + , and 3 H/ 3 He calculated with different initializations together with the soft and stiff symmetry energies are listed in the right of Table 1 and plotted in Figure 8. The relative difference in N n /N p among different initializations is less than 5% while it is about 50% for different nuclear symmetry energies, indicating the influence of initialization on this observable can be neglected. The yields of 3 H and 3 He and their ratio are also found to be insensitive to the initialization. This is because the production of cluster is mainly determined by the behavior of the nuclear mean-field potential at lower densities. Both v n 2 /v p 2 (v n 2 -v p 2 ) and π − /π + display a distinct difference among different initializations. The value of v n 2 /v p 2 calculated with the soft symmetry energy varies from 0.52 to 1, while it varies between 0.98 and 1.84 for the stiff symmetry energy. The largest difference in v n 2 /v p 2 obtained with the soft and stiff symmetry energies is about 1.09 appeared in SameIN9, while the smallest one is only about 0.02 appeared in SameIN6. The value of π − /π + also varies largely, from 2.17 in SameIN6 to 3.72 in SameIN5 and SameIN7, for the soft symmetry energy. In most of the cases, π − /π + obtained with the soft symmetry energy is larger than that with the stiff one, but two exceptions can be found in calculations with SameIN4 and SameIN6 which show the value of π − /π + is equal to or less than that obtained with stiff symmetry energy. The present calculation illustrates that fingerprint of nuclear symmetry energy on v n 2 /v p 2 and π − /π + can be either enhanced or reduced when different initializations are considered. Figure 8a,b display the elliptic flow difference and ratio between free neutrons and protons at mid-rapidity (|y 0 | < 0.1). Again, one finds that the nuclear symmetry energy effects can be either enhanced or reduced by different initializations. With SameIN6, the initialization effect cancels or even reverses the effect of symmetry energy on both v n 2 /v If one compares the v n 2 /v p 2 and π − /π + obtained with Normal-IN and DC, it is found that the introduction of density constraints alone in the initialization displays a negligible effect. While π − /π + calculated with DCPF is suppressed compared to the ratio obtained with Normal-IN. This suppression is a result of the reduced initial momentum fluctuation in DCPF mode. Moreover, the symmetry energy effect on v n 2 /v p 2 and v n 2 -v p 2 is slightly increased in the calculation with DCPF mode. It implies that the sampling of nucleon momentum in the initialization process should be more carefully and consistently treated when studying isospin-sensitive observables. The importance of initializations of transport model on these observables also has been found within different models [58][59][60][61][62].

Summary and Outlook
By artificially using the same initial nuclei in all collision events, the effect of the initial fluctuation on various observables is studied within the UrQMD model. It is found that for pion yield, the directed and elliptic flows of nucleons calculated with different initializations are different, because of the different density distribution in the compressed region, which is the result of the tiny difference in density distribution at the initial time. The yield ratio of free neutrons with respect to free protons N n /N p and the yield ratio between 3 H and 3 He clusters are found to be insensitive to different initializations, while the yield ratio between pions π − /π + and the elliptic flow ratio or difference between free neutrons and protons v n 2 /v p 2 (v n 2 -v p 2 ) displayed distinct differences among different initializations. It means that a detailed treatment on the initialization of the transport model is necessary when using these observables to constrain the density-dependent nuclear symmetry energy.
To further extend the present study to lower (where the nuclear liquid-gas phase transition may occur) or higher (where a phase transition from hadronic matter to the quark-gluon plasma (QGP) may occur) energies is also of particular interest. At lower energies, fluctuations and nuclear symmetry energy in nuclear fragmentation dynamics are correlated [63]. In relativistic heavy-ion collisions, observables like anisotropic flow, flow fluctuations and correlations, are expected to be sensitive to the initial fluctuations [64][65][66][67][68][69][70]. Because of the rapid development of both nuclear theory and experimental technology, our understanding of nuclei has become deeper, more and more experimental data about nuclear structures, e.g., charge density distribution, neutron skin and neutron halo, nuclear deformation, alpha-clustering, and high momentum tail, are available. It would be of particular interest to incorporate all these structures effects in the transport model, and to investigate whether isospin-sensitive observables can be affected or not.
Author Contributions: Investigation, Y.W. and Z.G.; writing-original draft preparation, Y.W.; writing-review and editing, Z.G. and Q.L. All authors have read and agreed to the published version of the manuscript.