Determine the neutron skin type by relativistic isobaric collisions

The effects of neutron skin on the multiplicity ($N_{\rm ch}$) and eccentricity($\epsilon_2$) in relativistic $^{96}_{44}$Ru+$^{96}_{44}$Ru and $^{96}_{40}$Zr+$^{96}_{40}$Zr collisions at $\sqrt{s_{_{\rm NN}}}=200$ GeV are investigated with the Trento model. It is found that the Ru+Ru/Zr+Zr ratios of the $N_{\rm ch}$ distributions and $\epsilon_{2}$ in mid-central collisions are exquisitely sensitive to the neutron skin type (skin vs.~halo). The state-of-the-art calculations by energy density functional theory (DFT) favor the halo-type neutron skin and can soon be confronted by experimental data. It is demonstrated that the halo-type density can serve as a good surrogate for the DFT density, and thus can be efficiently employed to probe nuclear deformities by using elliptic flow data in central collisions. We provide hereby a proof-of-principle venue to simultaneously determine the neutron skin type, thickness, and nuclear deformity.


Introduction
Nuclear densities are often parameterized by Woods-Saxon (WS) distributions. Proton and neutron distributions are usually not distinguished and the former, which is well measured, is substituted for the latter. For most studies in heavy ion collisions, this is sufficient. However, it is well known that proton and neutron distributions differ because of Coulomb interactions and the symmetry energy. This subtle difference becomes important when one compares collisions of similar species, such as isobar 96 44 Ru+ 96 44 Ru and 96 40 Zr+ 96 40 Zr collisions [1][2][3][4], recently conducted at the Relativistic Heavy Ion Collider (RHIC) [5]. Because Ru has more protons than Zr, the Ru charge distribution is fatter than the Zr's. Simply substituting the charge density for the mass density, one predicts a lower multiplicity in central Ru+Ru than Zr+Zr collisions [6]. Energy density functional theory (DFT) calculations indicate a significantly thicker neutron skin in 96 Zr than 96 Ru and that the overall size of 96 Zr is bigger than 96 Ru. Because of this, DFT instead predicts a higher multiplicity in central Ru+Ru than ZrZr collisions [1,2]. It is further found that, even incorporating the smaller Ru than Zr radius, the WS skin-type density parameterization cannot reproduce the DFT prediction [2]. In this paper, we demonstrate that the DFT densities can be significantly better modeled by the halo-type WS parameterization. The skinand halo-type densities predict quite different multiplicity distributions and elliptic flow, which can be easily discriminated by experimental data. If the data, soon to be available, favor the halo-type density (and therefore verifying the DFT calculation), then our study suggests that the halo-type WS density approximation can be used as a surrogate for DFT density for most applications in heavy ion collisions, including studies of nuclear deformity. This would provide an efficient way for such studies as DFT calculations of deformed nuclei are challenging.
The rest of the paper is organized as follows. Section 2 gives a brief description of the trento model used in this work. Section 3 discusses the effects of nuclear densities, focusing on the difference of neutron skin types and effect of the DFTcalculated densities. In Sec. 4, the effect of nuclear deformity is investigated. A summary is given in Sec. 5.

Model Setup
In relativistic heavy ion collisions, the charged hadron multiplicity (N ch ) distribution and elliptic flow (v 2 ) are two basic observables. For v 2 , one usually relies on macroscopic hydrodynamic or microscopic transport model calculations; however, those calculations are time-consuming. In this study, we use the trento model [7,8] to simulate Ru+Ru collisions and Zr+Zr collisions at nucleon-nucleon center-of-mass energy √ s NN = 200 GeV [4]. In trento the transverse area entropy density is proportional to the reduced thickness, where T A = ρ(r)dz is the nuclear thickness function. Particle production yield is proportional to the total entropy, N ch ∝ s(x, y)dxdy, and a Poisson-type fluctuation is used to generate the event-by-event N ch . We use the parameter p = 0 (i.e., s ∝ √ T A T B ), a gamma fluctuation parameter k = 1.4, and a Gaussian nucleon size of 0.6 fm, which were found to well describe the multiplicity data in heavy ion collisions [7,8]. The v 2 is approximately proportional to the initial eccentricity [9], 2 , at small amplitude ( 2 < 0.5) [10]. The conversion (proportionality) factor depends on dynamics and is expected to be the same for the highly similar 96  {2}, and substitute that for the v 2 ratio. The 2 can be readily obtained from the initial geometry in the trento model by [11] 2 = r 2 e inφ s(x, y)dxdy where φ is the azimuthal angle of position (x, y) and · · · denotes the event average. We note that, while we focus on the trento model, there are several initial geometry models on market. The most traditional one is the two-component Glauber model where the N ch depends on the numbers of participants and binary collisions [12][13][14][15][16]. The essential difference of all those models lies in how N ch or s depends on the thickness function (Eq. (1)). The treatment of fluctuations in those models can also affect the magnitude of 2 in each collision system, as well as the relation between N ch and centrality [15,[17][18][19][20]. Since we are primarily interested in the relative difference between the isobar systems, these model dependencies are weak. Transport models are more complex, but the essential ingredient relevant for our study is, again, the effective relationship of Eq. (1). Studies have previously been carried out by a multi-phase transport (AMPT) model [2], and consistent results have been obtained as we will refer to along with our results.

Effects of Nuclear densities
Besides some model dependence in N ch , the calculation boils down to the thickness function T A . Different nuclear densities yield different T A ; thus the N ch distribution and v 2 can be used to determine the nuclear density. Using the relative measures of isobar systems reduces the large uncertainties in theoretical modeling of the collision dynamics. Our starting point is the measured charge distributions of Ru and Zr from experimental data [21].

WS charge densities
The charge densities of Ru and Zr were parameterized by WS distribution (also known as two-parameter Fermi distribution), with R = 5.085 fm and 5.021 fm being the radius parameter for Ru and Zr, respectively, and a = 2.3/(4 ln(3)) = 0.523 fm for both [21] (the numbers are also tabulated in Table 1). The simplest approach is to take the proton and neutron WS parameters (R, a) to be as same as those for the charge density; the neutron skin thickness is implicitly zero. Figure 1 shows the charge density distribution (the dotted curve in panel a), the Ru+Ru/Zr+Zr ratio of the N ch distributions (the open squares in panel b) and the 2 ratio (the open squares in panel c) in Ru+Ru over Zr+Zr collisions. It is found Table 1: Skin-and halo-type WS parameterizations (radius parameter R and diffuseness parameter a) of neutron densities based on the measured charge (proton) densities [21] together with neutron skin thicknesses (taken to be ∆r np =0 and 0.12 fm [22,23] that the effect of the finite R difference between Ru and Zr is mostly small. The large effect (bending from unity down to zero) in the N ch distribution ratio at high N ch is due to the fact that the effective radius of Zr is smaller than that of Ru, which makes the overlap region denser (larger T A ) yielding a larger probability of high N ch events in central Zr+Zr than Ru+Ru collisions. Similar result was found by Ref. [6] where the charge WS densities were used. The effect of the R difference on the v 2 ratio is less than 0.5%.

Effect of neutron skin
Intermediate to heavy nuclei have sizeable neutron skin thicknesses, ∆r np = r 2 n 1/2 − r 2 p 1/2 , where r 2 1/2 is root-meansquare radius of either the proton or neutron distribution, There are two extreme cases to constrain the WS parameters (R n , a n ) for the neutron density with finite ∆r np [22]: (i) skintype density where R n > R p and a n = a p ; and (ii) halo-type density where R n = R p and a n > a p . Experimental data suggest a neutron skin thickness of ∆r np = 0.12 ± 0.03 fm for 96 Zr and it is of the halo-type [22,23]; there may be additional significant theoretical uncertainty on its magnitude [24]. We take ∆r np =0.12 fm for 96 Zr and keep it zero for Ru, and investigate whether relativistic heavy ion collisions can distinguish between the two types. If so, then it would add significant strength to our nuclear structure knowledge as relativistic isobar collisions are distinctively different from traditional methods. The corresponding Zr neutron density distributions are shown in Fig. 1(a) by the dashed blue and solid red curves. The N ch and 2 ratios from trento simulations using the skin-type neutron densities are shown by the blue triangles in Fig. 1(b) and (c), respectively. These ratios flip about unity from those with the charge WS densities in the previous section. This is easy to understand because the neutron skin has made Zr larger than Ru, opposite to that from the charge WS density (see Table 1). The now smaller Ru+Ru than Zr+Zr collisions produce more particles in central collisions, making the Ru+Ru/Zr+Zr ratio larger than unity. As previously discussed the strongest effect is on the tail of the N ch distribution. The effects are minor on 2 and the N ch distribution in non-central collisions. Similar results have also been found in our previous study with AMPT [2].   The halo-type neutron density, on the other hand, has profound effects on the N ch and 2 in mid-central isobar collisions. This is shown in Fig. 1(b) and (c) by the red circles. The increasing N ch ratio in central collisions comes from the smaller effective radius of Ru than Zr, similar to the case of skin-type neutron density. The ratios in mid-central collisions are distinct, with an arch shape reaching as large as 2% on the top and dipping below unity in peripheral collisions. This feature will be discussed further in the next section.

DFT-calculated densities
The skin-and halo-type neutron densities are two extreme parameterizations; the truth is likely in-between. Nowadays, DFT is the state-of-the-art in calculating nuclear structures [25][26][27][28]. Figure 2(a) shows the DFT densities of Zr in dashed curves, calculated with the well known skyrme parameter set SLy4 [25] (the corresponding individual neutron and proton densities have already been shown in Ref. [1]). Figure 2(b) and (c) show the Ru+Ru/Zr+Zr ratios of the N ch distributions and 2 , respectively, simulated by the trento model taking the DFT densities as input. It is immediately clear that the ratios are similar to those using halo-type neutron densities in Fig. 1. Similar features have also been observed with Glauber as well as AMPT calculations taking DFT densities as input [1,2]. This suggests, model-independently, that the DFT densities may have halotype neutron skin. In Ref. [2] we have tested a skin-type WS density keeping a fixed and adjusting R to the DFT-calculated density via r 2 of Eq. 5. The DFT-adjusted skin-type WS density yielded results similar to the skin-type neutron density results in Fig. 1(b) and (c). This indicates that the DFT density cannot be substituted by a skin-type density even if the effective radius is forced to be the same.
The DFT density contains more detailed information than the WS parameterization. Many of those details may not be important for a given observable of interest. The question arises whether a proper WS approximation of the DFT density would be sufficient for our N ch and 2 studies. Clearly, skin-and halo- Table 2: The WS parameterizations (radius parameter R and diffuseness parameter a) of proton and neutron (and nucleon) density distributions for 96 Ru and 96 Zr, matching to the corresponding r and r 2 from the DFT-calculated spherical densities with SLy4 skyrme parameter set [1,25]. The WS parameterization of nucleon density assuming a quadrupole deformity parameter β 2 = 0.16 and matching to the spherical DFT density is also listed. All quoted numbers are in fm. 96  type neutron densities do make significant differences as shown in Fig. 1, so not any parameterization would work. We can determine both WS parameters R and a by matching the r and r 2 from the DFT proton and neutron densities [29]. The obtained parameters are tabulated in Table 2. The parameters do indicate that the DFT density is closer to the halo-type WS density (R n ≈ R p but a n > a p ). This is in line with the findings by low-energy nuclear experiments [22,23]. We note that, because of the isospin symmetry of the strong interaction, most observable in relativistic heavy ion collisions, including the N ch and v 2 we study here, are sensitive only to the total nucleon density, the sum of the proton and neutron ones. Since what is experimentally measured is typically the charge (proton) density, how the neutron density is implemented is important. In our trento model simulations, we use the resultant total nucleon densities; their WS parameters are also tabulated in Table 2. Figure 2 shows the Ru+Ru/Zr+Zr ratios of the N ch distributions and 2 calculated by the trento model with WS densities parameterized from the DFT-calculated ones, and compares them to those calculated directly from the DFT densities. It is   found that the DFT densities and the WS parameterized ones give essentially the same Ru+Ru/Zr+Zr ratios. The two densities do yield slightly different results for each individual Ru+Ru or Zr+Zr system, up to a couple of percent (see the insets), but these differences cancel in the Ru+Ru/Zr+Zr ratios. This indicates that the properly parameterized WS densities are indeed sufficient to study the N ch and 2 . Note that the ratios in Fig. 2 are somewhat larger than the ones from the halo-type densities in Fig. 1. This is because the neutron skin thickness of Zr from the DFT calculation is ∆r np =0.16 fm, larger than the 0.12 fm used in Fig. 1. It implies that those ratios in non-central collisions can further be used to quantitatively measure the neutron skin thickness. The halo-type neutron densities from DFT indicate that Zr has a larger diffuseness parameter a and rms radii than Ru. These yield distinctive shapes of Ru+Ru/Zr+Zr ratios of the N ch distributions and 2 .
For this study we have generated 50 million events for each collision species in each density case. Total over 2 billion events have been collected for each isobar system by the STAR experiment at RHIC [5], so the statistical uncertainties would be negligible. Since most of the systematic uncertainties cancel in those ratios between the two isobar systems, our results can be readily compared against experimental data to yield valuable information on the neutron skin type and thickness. If experimental data, on the other hand, do not favor our DFT and halotype neutron density predictions, but rather the skin-type neutron density ones, then an overhaul of our current knowledge of nuclear structure would be called for.

Effect of nuclear deformation
So far, we have focused on the discussion of spherical nuclei. The elliptic flow in most central collisions are extra sensitive to nuclear deformation, as deformed nuclei colliding at zero impact parameter can have large eccentricity depending on the collision orientation [30,31]. To take nuclear deformation into account, the following extended WS formula [32], is used, where Y 0 2 (θ) = 1 4 5 π (3 cos 2 θ − 1) and β 2 is the quadrupole deformity parameter. DFT calculations of deformed nuclei are more challenging and often result in large uncertainties on β 2 . Experimental measurements of β 2 also suffer from large uncertainties. 96 Ru (β Ru 2 = 0.158) could be more deformed than 96 Zr (β Zr 2 = 0.08), or could be less deformed (β Ru 2 = 0.053, β Zr 2 = 0.217) [6]. In the previous section, we have demonstrated that DFT densities can be adequately approximated by halo-type neutron densities for comparative studies of N ch and 2 in isobar collisions. We assume this preserves for deformed Ru and Zr; this assumption is important as it allows us to bypass the challenging DFT calculations to still investigate efficiently the effects of nuclear deformity. With a given β 2 , we determine R and a of the proton and neutron densities, respectively, by matching the r and r 2 (Eq. (5)) of the corresponding spherical DFT densities. We set β 2 = 0.16 for this initial study; the corresponding WS parameters are listed in Table 2. The finite β 2 reduces the a parameter but have negligible effect on R. Figure 3 shows the Ru+Ru/Zr+Zr ratios of the N ch distributions and 2 with various combinations of Ru and Zr deformities. The N ch distribution ratio changes insignificantly due to the finite β 2 ; the change comes primarily from that in a, e.g. a finite β 2 for Zr reduces the a value yielding a larger N ch in Zr+Zr collisions compared to the spherical case and a smaller ratio of the N ch distributions. The 2 ratios, with the same β 2 value for both Ru and Zr, are almost identical for the spherical (β 2 = 0) and deformed (β 2 = 0.16) cases, except perhaps in very central collisions. If the two nuclei differ in β 2 , then the 2 ratio is significantly impacted in central (and semi-central) collisions, as expected [33]. However, the effect of deformation on the 2 ratio in peripheral and semi-peripheral collisions is insignificant.  This feature indicates that the shape of the 2 ratio in peripheral and semi-peripheral collisions can still be used to determine the neutron skin type and thickness as we demonstrated in the previous section with spherical nuclei. As aforementioned, 40 times more statistics have been collected by STAR than used in this simulation [5]. Even though experimental triggers suffer inefficiencies in peripheral collisions, it is significant only in very peripheral collisions (N ch 10), and even there the inefficiency is larger than 10%. In other words, despite the large bias to the v 2 ratio in (semi-)central collisions, our proposal to discriminate the neutron skin by Ru+Ru/Zr+Zr ratios of N ch distributions and v 2 is still feasible even though the Ru and Zr deformities are largely uncertain. Furthermore, the v 2 ratio in central collisions, together with N ch distribution ratio, can be exploited to probe the β 2 parameters of the isobar nuclei; we postpone such a study to further work.

Summary and outlook
By using the trento model, we have investigated the effects of neutron skin on the ratios of the charged hadron multiplicity (N ch ) distributions and eccentricity ( 2 )  found that these Ru+Ru/Zr+Zr ratios are exquisitely sensitive to the neutron skin type (skin vs. halo) and thickness. Taking the Woods-Saxon (WS) neutron densities to be as same as the proton densities for both nuclei (i.e. without any neutron skin), the Ru+Ru/Zr+Zr ratio of the N ch distributions bends down below unity at large N ch because of the smaller effective Zr radius than the Ru's in these charge distributions. Including neutron skin, which is thicker in Zr than in Ru, the ratio bends up above unity at large N ch because the effective nucleon radius of Zr is now larger than the Ru's. The shapes of the Ru+Ru/Zr+Zr ratios of the N ch distributions and 2 in mid-central collisions can further distinguish between skin-type and halo-type neutron densities, both having the same skin thickness (see Fig. 1). The Ru+Ru/Zr+Zr ratios obtained by using nuclear densities from the state-of-the-art calculations by energy density functional theory (DFT) show similar features as those using nuclear WS density distributions with a halo-type neutron skin (see Fig. 2 and Ref. [1,2]). This is because the DFT-calculated densities are more like the halo-type than the skin-type (see Table 2). The soon-to-be-available isobar data should be able to determine whether the neutron skin is skin-type or halo-type, how large the neutron skin thickness (∆r np ) is, and whether or not the DFT densities are a truthful reflection of nuclear structures.
After demonstrating that the halo-type WS densities can substitute the DFT-calculated ones to faithfully describe the Ru+Ru/Zr+Zr ratios, we include a nuclear quadrupole deformity into the WS distribution of Ru and/or Zr to study the effects. We show, as a proof of principle, that the ratios of the N ch distributions and v 2 in isobar collisions can be used to probe the nuclear deformation (see Fig. 3). Corroborating with other means of probing nuclear deformity, e.g. v 2p T correlations [34], should yield valuable information on nuclear structure in the future.