Parameterization of Deformed Nuclei for Glauber Modeling in Relativistic Heavy Ion Collisions

The density distributions of large nuclei are typically modeled with a Woods-Saxon distribution characterized by a radius $R_{0}$ and skin depth $a$. Deformation parameters $\beta$ are then introduced to describe non-spherical nuclei using an expansion in spherical harmonics $R_{0}(1+\beta_2Y^0_2+\beta_4Y^0_4)$. But when a nucleus is non-spherical, the $R_{0}$ and $a$ inferred from electron scattering experiments that integrate over all nuclear orientations cannot be used directly as the parameters in the Woods-Saxon distribution. In addition, the $\beta_2$ values typically derived from the reduced electric quadrupole transition probability B(E2)$\uparrow$ are not directly related to the $\beta_2$ values used in the spherical harmonic expansion. B(E2)$\uparrow$ is more accurately related to the intrinsic quadrupole moment $Q_{0}$ than to $\beta_2$. One can however calculate $Q_0$ for a given $\beta_2$ and then derive B(E2)$\uparrow$ from $Q_0$. In this paper we calculate and tabulate the $R_0$, $a$, and $\beta_2$ values that when used in a Woods-Saxon distribution, will give results consistent with electron scattering data. We then present calculations of the eccentricity $\varepsilon_2$ and $\varepsilon_3$ with the new and old parameters. We demonstrate that $\varepsilon_3$ is particularly sensitive to $a$ and argue that using the incorrect value of $a$ has important implications for the extraction of $\eta/s$ from the QGP created in Heavy Ion collisions.


Introduction
In relativistic nucleus-nucleus collisions, the geometry of the initial overlap region is reflected in the final momentum space distributions of produced particles [1][2][3]. How much that geometry is translated into the final state distributions is used to infer information about the properties of the matter created in the collision fireball like its viscosity [4]. The initial geometry plays a particularly important role in interpreting the data and in extracting the viscosity to entropy ratio η/s. It is important therefore to understand the initial conditions including the exact shape of the colliding nuclei. The inference of the properties of the fireball from data is hindered by uncertainties in the characteristics of the initial state [5]. Recently collisions between Uranium nuclei ( 238 U) which have an intrinsic prolate shape [6], have been used as a way to manipulate this initial geometry in order better test our understanding of the initial state of heavy ion collisions and the subsequent fireball [7,8].
An important part of describing the initial conditions is to correctly model the geometry of the incoming nuclei. For many years in simulations for heavy ion collisions, nuclei were approximated as smooth density dis-tributions and the only anisotropies considered in the initial state were the intrinsic almond shape caused by the overlap of two spherical nuclei. As the accumulation of RHIC data gradually demonstrated that final state anisotropies were sensitive to the initial geometry and its fluctuations, it became necessary to take into account the lumpiness of the colliding nuclei [9][10][11][12][13][14]. This is done through Monte-Carlo simulations (M-C) where each nucleus is generated with a finite number of nucleons distributed with a density ρ described by a Woods-Saxon distribution [15]: where ρ 0 is the density at the center of the nucleus. The nuclear radius R 0 and skin depth a are commonly taken from high-energy electron scattering measurements [16]. For non-spherical nuclei, this description was extended by introducing spherical harmonics in the Woods-Saxon distribution to describe the modulation of the nuclear radius with θ [17]: where β 2 and β 4 are the deformation parameters. β 2 is often derived from measurements of the reduced electric quadrupole transition probability B(E2)↑ from the 0 + ground state to the first 2 + state [6] according to the formula: R 0 is taken to be 1.2A 1/3 . The B(E2)↑ values are directly measured experimental quantities, but the derivation of β 2 from B(E2)↑ is usually done with model dependent assumptions including the assumption that the nuclear charge distribution is a hard edged, step-function rather than a Woods-Saxon distribution. In order to check for consistency with the measured B(E2)↑ values, one can vary the β 2 used in Eq. (2) and then calculate the intrinsic quadrupole moment of the resulting nucleus: and check to ensure that Q 20 is consistent with the measured B(E2)↑ according to the approximation [18]: This procedure ensures that the deformation of the simulated nucleus is consistent with the measured B(E2)↑ values for the given model describing the density profile of the nucleus. Another difficulty in characterizing deformed nuclei comes from a mismatch between the parameters inferred from electron scattering experiments and the Woods-Saxon parameters used in the simulation of the nuclear density profile. Electron scattering experiments probe the spherical part of the density distribution (characterized by radius R 0 and diffuseness a), which is averaged over all orientations of the nucleus. If the nucleus is not spherical, then the R 0 and a used in Eq. (2) will not necessarily correspond to the average radius and average diffuseness inferred from the electron scattering experiments [16]. For this reason, it is important to calculate the R 0 , a, and β 2 values that when used in conjunction with Eq. (2) yield a nucleus that is consistent with the experimental measurements of B(E2)↑, R 0 and a. In this paper we present updated R 0 and a for nuclei commonly used in collisions at RHIC or the LHC, 238 U, 208 Pb and 197 Au, as well as β 2 value for 238 U. We then show calculations of the second and third harmonic participant eccentricity ε 2 and ε 3 for the new and old parameters. We find that in addition to the obvious dependence of ε 2 on β 2 , the ratio of ε 3 /ε 2 is very sensitive to the diffuseness parameter a: ε 3 increases with increasing a while ε 2 decreases so that ε 3 /ε 2 is overestimated if a is overestimated. Since viscous damping decreases the ratio of triangular flow over elliptic flow (v 3 /v 2 ), an overestimate of ε 3 /ε 2 will lead to an over-estimate of the amount of viscous damping needed to match the experimental data. We find therefore that the new values of a presented in this work should lead to a decrease in the value of η/s inferred from model-to-data comparisons.

Parameterization of Deformed Nuclei
When using a Woods-Saxon distribution with parameters R 0 , a, and β 2 to characterize the density distribution of a deformed nucleus, it is obvious that after allowing for all possible rotations of the nucleus, the averaged density distribution will still be accurately described by a Woods-Saxon distribution. We find however, that for the range of R 0 , a, and β 2 values describing typical nuclei, the final density distribution after averaging over all appropriate rotations is indeed well described by a Woods-Saxon Distribution but with a different radius R 0 and diffuseness a . To accurately model a nucleus, R 0 and a should match the parameters measured in electron scattering experiments. The larger β 2 is, however, the more R 0 and a deviate from the R 0 and a used in Eq. (2).
Taking ρ(r, θ) to be the density distribution of a deformed nucleus centered at the origin of the spherical coordinate system, then with n-times randomly rotated and overlapping ellipsoids, the total nucleon numbers of these n ellipsoids reads: n ρ(r, θ)r 2 sin(θ) dr dθ dφ.
On the other hand, such total nucleon numbers can also be given by directly n-times integrating a density of the space, f , which only depends on r: Hence, letting Eq. (6) equal Eq. (7), one can easily identify f (r) as: Namely, integrating ρ(r, θ)sin(θ)/2 over [0, π] will give the nucleon density as a function of r, as shown in Fig. (1) (a.). In order to correctly reproduce the experimentally measured B(E2)↑ (Q 20 ) via Eq. (4) and (5), one has to use the charge density distribution to represent the charge profile of the nucleus as realistically as possible. In this work, we assume that the charge density of nuclei can be calculated by folding the proton (neutron) density with the charge distribution of the proton (neutron): where ρ p (ρ n ) represents proton (neutron) density, while n p (n n ) is the charge distribution of single proton (neutron). It's worth mentioning that the M-C Glauber simulation for the collision of two nuclei does not typically distinguish neutrons and protons, so we assume both neutrons and protons follow the same Woods-Saxon distribution in the form of Eq. (2) for simplicity, however, with different normalization factors: where Z, N and A denote proton number, neutron number and nucleon number, respectively. Based on this assumption, one can utilize B(E2)↑ (Q 20 ) as a reference to optimize the Woods-Saxon parameters and then apply them to generate nuclei via a M-C simulation. Often, in Eq. (2), all nucleons are considered as point-like, viz. n p and n n in Eq. (9) are described by delta functions. In a realistic case however, the finite size of the nucleon should be taken into account [19][20][21]. For protons, we adopt the experimentally supported dipole form of the charge profile [22], which reads: with m = √ 12/0.877 reproducing the proton RMSradius 0.877 fm. For neutrons, we use the following form [23]: with r 2 n = −0.113 fm 2 and r 1 = √ 2/5·0.71 fm. B(E2)↑ (Q 20 ) can then be obtained by using the above euqations. However, our calculation shows that the contribution to Q 20 from neutrons is very small (much less than 1%), so we neglect the charge form factor of neutrons and only take protons into account. In this case, Eq. (9) can be simplified to: where r represents the distance from the center of the nucleus (coordinate origin) to any possible charge points in the proton, while r represents the distance to the center of the proton. After integrating Eq. (13) numerically and employing Eq. (8), an integrated charge distribution of deformed nuclei that only depends on r can be obtained.
In the steps above, directly adopting parameters from e-A scattering is no longer valid. Fig. (1) (b.) shows the charge density distribution of uranium and gold nuclei with a point-like and a finite-size profile. The same parameters of R 0 and a from Ref. [16] are taken for both profiles for each nucleus. Deviations between the two profiles can been seen from this figure implying the parameter sets for finite-size profiles need to be adjusted. A relatively larger radius R 0 and a smaller diffuseness a for the finite-size scenario is needed to reproduce the Woods-Saxon distribution inferred from electron scattering measurements.  Now taking into account the finite size of nucleons and the presence of a deformation β 2 , we want to find the set of parameters β 2 , a, and R 0 that when used in Eq.
(2) will yield a nucleus consistent with data on B(E2)↑ and electron scattering. During the implementation, we loop over all possible combinations of R 0 , a and β 2 . For each set, the deformed Woods-Saxon function with finite size profile is first rotated and then averaged. After that, the charge distribution with respect to the radius is extracted. Meanwhile, Q 20 is directly calculated from Eq. (4) and (13) in order to check for consistency with B(E2)↑. An optimal set should guarantee that (i) the charge distribution from electron scattering data is reproduced, (ii) the B(E2)↑ (Q 20 ) value from experiments can be obtained. For the first criterion, the sum of squares of the residual (S S R) is used to quantify the difference between the two distributions. In this work, only parameter sets providing minimum S S R and the nearest B(E2)↑ (Q 20 ) are considered in the final answer. We didn't tune β 4 because with its typical value, including β 4 makes little difference to our result.

Results
Fig. (2) shows an example of scanning optimal parameter sets for 238 U. According to Ref. [6], the B(E2)↑ value of 238 U is 12.09±0.2 (e 2 b 2 ). Within this experimental uncertainty, we calculate all possible combinations of R 0 , a and β 2 to find the minimum S S R. Fig.  (2) (a.) shows the two dimension scan of the parameter sets. The anti-correlation between a and R 0 is due to the requirement that the calculated B(E2)↑ from this procedure has to match previous measurement of B(E2)↑, of which the uncertainty is reflected as finite band width in the plot. The minimum SSR value visible in the plot corresponds to the purple region centered around the best parameters R 0 =6.86 fm and a = 0.42 fm. The corresponding β 2 value that gives a nucleus with a Q 20 consistent with the B(E2)↑ is 0.265. Fig. (2) (b.) shows that the Woods-Saxon distribution is recreated by the above parameters with nice consistency with the one extracted from electron scattering. We summarize the parameter sets for different nuclear species in Tab. (1). For 208 Pb, the ground state nucleus is considered spherical [6] although its excited state may show a deformed configuration, thus the β 2 ( 208 Pb) is taken to be 0 in this study. Namely, the only contributing step for 208 Pb in the aforementioned procedure is the finite-size correction. For 197 Au, B(E2)↑ has not been reported from experiments, therefore β 2 remains a major uncertainty in the initial conditions for Au+Au collisions. Therefore we fix the commonly used value β 2 ( 197 Au)=-0.13, which comes from a mix of data and model calculations [24], and only provide the updated R 0 and a values with the caveat that these values will change if the β 2 is changed.

Finite size + Rotation
Two-parameter Fermi

Eccentricity
To test how the corrections to the Woods-Saxon parameters may affect observables, we study the multiplicity distribution, the second and third harmonic participant eccentricity (ε 2 and ε 3 ) in U+U collisions from a Glauber model similar to the one discussed in Ref. [25] where N part is the number of struck nucleons, N bin is the number of binary nucleon-nucleon collisions, and x hard is a fractional contribution of N bin to the multiplicity [15,26]. The multiplicity is then generated by sampling a negative binomial distribution n AA times with negative binomial parameters (n pp and k) taken from p+p collisions at the same energy and in the same |η| window [27]. For this calculation we use the parameters n pp = 2.43 and k=2 for the negative binomial and x hard = 0.13.  Glauber model calculations. The maximum dN/dη with the new parameters is slightly higher than that with the original parameters. The increase comes from the smaller skin depth in the new parameters. As the skin depth becomes larger, the nucleus becomes more diffuse leading to a smaller N bin and smaller estimate for dN/dη. To study the impact of the new Woods-Saxon parameters on the initial geometry, we calculate the second and third harmonic participant eccentricity. Since different input parameters generate different multiplicity distributions, we present the eccentricities as a function of centrality intervals based on the percentage of the to-tal multiplicity. Fig. (4) upper and middle panels show ε 2 and ε 3 for the two parameters sets. The lower panel shows the ratio of the results with the parameters over the results with old parameters for ε 2 and ε 3 . For the most central collisions, both ε 2 and ε 3 ratios are below one. For the most central bin, the initial geometry is most sensitive to β 2 so the smaller β 2 values in the new parameter set lead to a smaller ε 2 and ε 3 . In mid-central collisions, however, the new parameter set produces larger ε 2 values but smaller ε 3 values. This behavior can be traced to the smaller value of a in the new parameter set. In non-central collisions, the value of 3 is enhanced by the probability of a nucleon fluctuating out on the edge of one nucleus and impinging on the center of the other nucleus where it encounters a relatively large number of nucleons. This effect will be largest when the nucleon from nucleus A fluctuates in the reaction plane towards nucleus B. That configuration enhances ε 3 but decreases ε 2 . It is this effect that explains the rise and fall of the ridge correlation in heavy ion collisions [13] and it also causes a correlation between the second and third harmonic event planes [28]. To further confirm this explanation, in  From this study, it is clear that the centrality dependence of ε 3 is strongly dependent on the diffuseness. For that reason, using the correct value of the diffuseness parameter when modeling heavy ion collision is crucial especially when calculating v 3 = cos(3φ) where φ is the azimuth angle of each particle relative to the major axis of the third harmonic event anisotropy. Since the relationship of v 3 and v 2 similarly defined are often used to estimate the viscosity to entropy ratio η/s, estimates of η/s from model-to-data comparisons can be adversely affected if the model does not use the correct Woods-Saxon parameters.

Summary
We find that when modeling the density distribution of deformed nuclei with a Woods-Saxon distribution, the radius R 0 and skin depth a used in the model will be substantially different than the average radius R 0 and skin depth a that would be observed in electron scattering experiments. The more deformed the nucleus, the larger the discrepancy. For this reason, one must modify the parameters used in the Woods-Saxon model so that after appropriately averaging over all orientations of the axis-of-symmetry for the nucleus, the average radius R 0 and skin depth a match the values reported from electron scattering experiments. In addition, one should also take into account the radius of the nucleon when carrying out a Monte-Carlo simulation of the positions of nucleons inside the nucleus, otherwise, the effective radius will be larger than the simulated radius by roughly the size of the nucleon. We also find that the model dependent β 2 parameters estimated from B(E2)↑ measurements, lead to an overestimate of the deformation of nuclei when they are used in a deformed Woods-Saxon distribution. We presented a procedure to calculate the correct values of the Woods-Saxon distribution for deformed nuclei (R 0 , a, and β 2 ) so that the resulting nucleus is consistent with electron scattering data and B(E2)↑ measurements. Our calculations show that, for 238 U, 208 Pb and 197 Au, the new value of R 0 is slightly larger than the old value while the new value of a is significantly smaller. For 238 U, the obtained β 2 is also different than the ones that are commonly used. We also presented the results from Glauber-Model Calculations for U+U collisions to study the effect of the new parameters. The decrease in the skin depth has a large impact on eccentricity, increasing ε 2 but decreasing ε 3 . Since the amount of viscous damping needed for a model to match v 3 /v 2 data will depend strongly on ε 3 /ε 2 , overestimates of a will lead to over-estimates of η/s. Predictions of models using the old parameter values may need to be revisited [29].