Introduction

The recent isobar data [1, 2] stimulated a wide interest from the relativistic heavy ion community in the physics of nuclear structures of isobars (nuclei with equal number of nucleons but different numbers of protons and neutrons) [3-12]. The isobar collisions of ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ were initially proposed to search for the chiral magnetic effect (CME) [13-16]. The measurements in isobar collisions at $\sqrt{s_{{}_{\text{NN}}}}=200$ GeV by the STAR Collaborations showed sizable differences in the multiplicity distribution, elliptic flow, and triangular flow [1], indicating differences in the density profiles of the colliding nuclei [3, 6]. Nuclear density is usually obtained from theory calculations, such as energy density functional theory (DFT) [17-20], incorporating information from experimental measurements. In relativistic heavy ion collision simulations, it is common to use the parameterized nuclear density in the Woods-Saxon formula, $\rho \left(r\right)=\frac{{\rho }_0}{\exp \left[\left(r-R(1+{\beta }_{2}Y{2}^0+{\beta }_3{Y}_3^0+\mathrm{...}\right)/a\right]+1}\mathrm{,}$ (1) with the radius parameter R, the diffuseness parameter a, and the quadrupole (octupole) deformation parameter β₂ (β₃). A large a value in ⁹⁶Zr is suggested by the STAR data on multiplicity distribution and elliptic flow ratios in non-central collisions [1], consistent with a halo-type neutron skin thickness as predicted by the DFT calculations [21]. The flow ratios observed in the most central collisions suggest that the ⁹⁶Ru has a larger β₂, while ⁹⁶Zr has a larger β₃ [6].

Much progress has been made in quantitatively assessing the differences in size and shape of isobars by using relativistic isobar collisions [5, 7-9, 12, 22-26]. These differences can place stringent constraints on key parameters in theoretical modeling of nuclear interactions, which is the backbone of all theoretical calculations of nuclear structure including DFT. One of such parameters is the slope parameter $L(\rho )$ of the nuclear symmetry energy as a function of nuclear density, which cannot yet be well determined by low-energy nuclear experiments [27-35]. In addition, nuclear deformations are largely uncertain from nuclear structure calculations [36, 37, 15]. On the other hand, the STAR experiment has collected approximately 2 billion minimum-bias events for each species of ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions at $\sqrt{s_{{}_{\text{NN}}}}=200\text{ GeV}$, with good control of systematic uncertainties [1]. These provide an unique opportunity to study the density slope parameter of the symmetry energy and the deformation parameters of the isobars. Previous studies suggest that ratio observables in isobar collisions are less sensitive to final-state interactions [5, 21, 9]. In this work, we investigate the effect of initial-state fluctuations on those ratio observables in relativistic isobar collisions.

We use the Optical Glauber model and Monte Carlo Glauber model [38-40] to investigate the effect of nuclear deformations and initial fluctuations on the relative differences $R_{\langle N_{\text{ch}}\rangle }$, $R_{{ϵ}_2}$, and $R_{{ϵ}_3}$ in relativistic isobar collisions, where $R_X\equiv 2\frac{X_{\text{RuRu}}-X_{\text{ZrZr}}}{X_{\text{RuRu}}+X_{\text{ZrZr}}}\mathrm{.}$ (2) The Glauber model is widely used to determine the centrality in experiments, including the blind analysis of the STAR isobar data [1]. Due to symmetry, the elliptic flow v₂ with impact parameter b=0 and the triangular flow v₃ can not be generated in spherical nucleus-nucleus collisions with smooth initial conditions in Optical Glauber simulation. The collision configuration of deformed nuclei will contribute to spatial anisotropy, while event-by-event fluctuation contributions are also important. We further investigate the effect of initial fluctuations by comparing the results under the scenarios of nucleon participants and quark participants.

The rest of the paper is organized as follows. Section 2 gives a brief description of the models and the nuclear density profiles used in this work. Section 3 discusses the effects of nuclear deformation and initial fluctuations on $R_{\langle N_{\text{ch}}\rangle }$, $R_{{ϵ}_2}$, and $R_{{ϵ}_3}$. A summary is given in Sect. 4.

Glauber model and isobar nuclear densities

Particle production in relativistic heavy ion collisions is usually parameterized by the two-component model, $N_{\text{ch}}=n_{\text{pp}}\left[\left(1-x\right)N_{\text{part}}/2+x{N}_{\text{coll}}\right]\mathrm{,}$ (3) with the number of participants N_part and number of binary collisions N_coll obtained from the Glauber model [39]. Here the x is the relative contribution of hard processes and n_pp=N_ch in nucleon-nucleon collisions. For the case of x=0, the particle production is independent of hard processes and it can alternatively be described by a variant of the Glauber model, the so-called Trento model [41, 42], $N_{\text{ch}}\propto {\displaystyle \int }s\left(r_{\perp }\right)d^2{r}_{\perp }\propto {\displaystyle \int }{\left(\frac{T_A^p\left(r_{\perp }\right)+T_B^p\left(r_{\perp }\right)}{2}\right)}^{1/p}{d}^2{r}_{\perp }$ (4) with p=1 (we note that previous studies favor p~0 [42]). Here s is the entropy density, $r_{\perp }$ is the transverse radius, and T_A(B) is the reduced thickness function for the colliding nuclei A(B) with ${\displaystyle \int }\left(T_A+T_B\right)d^2{r}_{\perp }=N_{\text{part}}$. In this work, we focus on the centrality dependence of the ratio observables in isobar collisions, so the parameter n_pp and the normalization factor hidden in the above equation are not relevant to our study.

The elliptic flow (v₂), triangular flow (v₃) are the Fourier coefficients of the azimuthal angle (ϕ) distribution $v_n=\langle \cos \left[n\left(\varphi -{\Psi }_n\right)\right]\rangle$, where $\langle \mathrm{...}\rangle$ denotes the average over an event and Ψn is the plane angle in momentum phase. The large anisotropic flow can be successfully described by relativistic hydrodynamics with a large spatial eccentricity [43], ${ϵ}_n{e}^{in{\text{Φ}}_n}=\frac{{\displaystyle \int }{r}_{\perp }^n{e}^{in\phi }w\left(r_{\perp }\right)d^2{r}_{\perp }}{{\displaystyle \int }{r}_{\perp }^{n}w\left(r_{\perp }\right)d^2{r}_{\perp }}\mathrm{.}$ (5) Here $w\left(r_{\perp }\right)$ is the weight factor, often taken to be the entropy density s or the energy density e. Φn is the nth order plane angle in the configuration space. The anisotropic flow is approximately proportional to the initial eccentricity for small amplitude [44]. In this work, we use the participant density to calculate ϵ₂ and ϵ₃, however, our main conclusions do not change if alternative w is used.

For spherical nucleus-nucleus collisions, the multiplicity and eccentricities can be calculated analytically in the Optical Glauber model for a given impact parameter b. For deformed nucleus-nucleus collisions, we calculate the ratio observables at a given shape orientation with a specific Euler rotation and then take the average over all orientations $\langle X\left(b\right)\rangle =\frac{1}{{\left(4\pi \right)}^2}{\displaystyle \int }X\left({\Omega }_A\mathrm{,}{\Omega }_B;b\right)\text{d}{\Omega }_A\text{d}{\Omega }_B\mathrm{.}$ (6) Here X is the observable of interest and ${\Omega }_{A\left(B\right)}$ is the orientation solid angle for colliding nuclei A(B).

The event-by-event fluctuations are not included in the Optical Glauber model [39], and such an effect is crucial for the anisotropic flow observables, especially the odd-order ones. We, therefore, use the Monte Carlo Glauber model to include the fluctuation effect, where the finite number of nucleons (96 for the Ru and Zr isobars) are sampled in each colliding nucleus. The event-by-event fluctuations depend on the number of nucleons in the colliding nuclei, so these effects are expected to be larger in the isobar collisions than in Au+Au collisions. After sampling the nucleons with the appropriate nuclear density, the multiplicity is calculated based on Eq. (3), and the eccentricity is calculated from the transverse position of those finite participating nucleons. Here the participants and binary collisions are determined from nucleon-nucleon interactions which happens when the distance between the two nucleons, each from the incoming and outgoing nuclei, is less than $\sqrt{{\sigma }_{\text{NN}}/\pi }$ [39], with the inelastic cross section set to ${\sigma }_{\text{NN}}=42$ mb for $\sqrt{s_{{}_{\text{NN}}}}=200$ GeV. Poisson-type multiplicity fluctuations are applied in the Monte Carlo simulations.

The density profiles of ⁹⁶Ru and ⁹⁶Zr nuclei are required in the Glauber simulations. In this study, the parameters for the Woods-Saxon density distributions of ⁹⁶Ru and ⁹⁶Zr are taken from Ref. [22, 45] (see Tab. 1), which are extracted from the isobar densities obtained by DFT calculations with the symmetry energy $E_{\text{sym}}\left({\rho }_c\right)=26.65$ MeV and its density slope parameter $L({\rho }_c)=47.3$ MeV at subsaturation cross density ${\rho }_{\text{c}}\sim 0.11$ fm^-3 [20]. Nuclear deformation is not included in the DFT calculations due to its strong model dependence. Following our previous work [22, 45], to introduce the nuclear deformation, the WS parameters R and a for the given β₂ (or β₃) are calculated to match the volume and root-mean-square (RMS) of the corresponding nucleus calculated by DFT, keeping the normalization factor ρ₀ fixed.

Centrality is determined by the multiplicity distributions in the Monte Carlo Glauber model, similar to what was done in heavy ion experiments [1]. For the optical Glauber model and optical Trento model in spherical nucleus-nucleus collisions, the centrality is related to the differential cross-section in impact parameter [39]. In this study, we use the 16 equal-size centrality bins from 0-80% centrality range. The results obtained with the fixed impact parameter corresponding to the centre of each centrality bin are used to mimic the observables in such a centrality bin. The impact parameters for ⁹⁶Ru and ⁹⁶Zr are listed in Tab. 2. The impact parameters are larger in ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions than in ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ collisions for a given centrality bin because the total cross section for ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions (${\sigma }_X^{\text{Z}rZr}=4.551$ barn) is larger than that for ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ collisions (${\sigma }_X^{\text{RuRu}}=4.369$ barn) based on the nuclear density parameters listed in Tab. 1. For simplicity, the impact parameter list is also used to determine the centrality for deformed nucleus-nucleus collisions in Optical Glauber simulations. We note that determining centrality in the deformed case is in principle more complicated, but this does not affect our conclusion (see discussion in the next section).

In order to achieve high precision of calculations in our study, the GPU parallel computing technology [46] is used in our Glauber model simulations.

Results and discussions

The $R_{\langle N_{\text{ch}}\rangle }$ obtained from the Optical Glauber model with x=0 and x=0.15, as well as those from the Optical Trento model with p=0 and p=-1, are shown in Fig 1. The case of the Optical Glauber simulation with x=0 is nothing but the Optical Trento simulation with p=1, and the multiplicity, in this case, is only proportional to the number of participants. The $R_{\langle N_{\text{ch}}\rangle }$ distribution depends on the model parameters x and p. The magnitudes increase with x and p over the whole centrality range. For the most central collisions, the $R_{\langle N_{\text{ch}}\rangle }$ increases with x in Optical Glauber simulations, but is independent of p in Optical Trento simulations. This is shown in the top 5% centrality range in Fig. 1(b). The magnitudes of $R_{\langle N_{\text{ch}}\rangle }$ converge to the same value from Optical Trento simulations with different values of p, but the trends in centrality differ.

Fig. 1Full-screen View

(Color online). The centrality dependence of $R_{\langle N_{\text{ch}}\rangle }$ from Optical Glauber (Opt.Glb) simulations with x=0.15 and x=0, and from Optical Trento (Opt.Trento) simulations with p=0 and p=-1 for (a) 0-80% centrality range and (b) top 5% centrality

The parameters x and p are usually constrained by the individual multiplicity distributions (in this case, of isobar collisions or their mean). The $R_{\langle N_{\text{ch}}\rangle }$ in the most central collisions has been proposed to probe the neutron skin and symmetry energy. However, the particle production mechanisms will introduce uncertainties as shown in previous study [5]. The results shown in Fig. 1 imply that the model uncertainty of $R_{\langle N_{\text{ch}}\rangle }$ is negligible if the particle production follows the Trento model. They also suggest that the precision measurements of the mean multiplicity distribution of ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions and their difference may offer insights on the particle production mechanisms in relativistic heavy ion collisions.

To study effects from event-by-event fluctuations, we present the $R_{\langle N_{\text{ch}}\rangle }$ distributions from Monte Carlo Glauber simulations with x=0.15 and x=0 in Fig. 2. The apparent large fluctuation effect (i.e. difference in $R_{\langle N_{\text{ch}}\rangle }$ between the Optical and the Monte Carlo Glauber) is because of the different total cross sections predicted by these two models [39]. For instance, the impact parameter cut would be smaller for a given centrality when a smaller total cross-section is used. Based on the nuclear density profiles in our study, the total cross sections for ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions are ${\sigma }_X^{\text{MC,RuRu}}=4.206$ barn and ${\sigma }_X^{\text{MC,ZrZr}}=4.348$ barn, which are 4.46% and 3.73% smaller than those in the Optical cases, respectively. To remove this trivial difference, we calculate the optical model using the total cross sections from the Monte Carlo Glauber. The results are shown in Fig. 2; they become similar now. Our results indicate that both the data of centrality dependent $R_{\langle N_{\text{ch}}\rangle }$ and the model prediction of total cross section are important for the model to data comparison in relativistic isobar collisions [1].

Fig. 2Full-screen View

(Color online). The centrality dependence of $R_{\langle N_{\text{ch}}\rangle }$ from Optical Glauber (Opt.Glb, solid curves) and Monte Carlo Glauber (MC.Glb, symbols) simulations with x=0.15 and x=0. The dashed curves are the results from Optical Glauber calculations with the total cross section from Mont Carlo Glauber simulations

The ${}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}$ nuclei are deformed [47, 48], as also suggested by recent studies [1, 6]. The quadrupole deformation parameter for ⁹⁶Ru is ${\beta }_{\mathrm{2,}Ru}=0.16$ and is negligible for ⁹⁶Zr. The octupole deformation parameter is negligible for ⁹⁶Ru and is ${\beta }_{\mathrm{3,}Zr}=0.20$ for ⁹⁶Zr. Figure 3 shows the $R_{\langle N_{\text{ch}}\rangle }$ from deformed isobar collisions, where the volume and RMS radius of the colliding nuclei are constrained to match those of spherical nuclei [22, 45]. For both the Optical and Monte Carlo calculations, the results from deformed isobar collisions are similar to those from spherical ones. Here we have used the same impact parameters listed in Tab 2 for deformed isobar collisions. As mentioned above, the centrality definition in the deformed case may be biased by multiplicity variations at a given impact parameter. However, the centrality is well defined in the Monte Carlo Glauber simulation, and the consistency in the corresponding results between spherical and deformed cases indicates that our centrality treatment in the Optical Glauber simulation is appropriate. The $R_{\langle N_{\text{ch}}\rangle }$ is considerably larger in the Optical case because of the aforementioned larger total cross section. The differences are ~10% between with and without nuclear deformation. Based on previous study [5], this would result in an uncertainty on the extracted slope parameter of symmetry energy of about 5 MeV. Our results indicate that, with appropriate magnitudes of deformation parameters, the effect of nuclear deformations on the uncertainties of $R_{\langle N_{\text{ch}}\rangle }$ distributions are small.

Fig. 3Full-screen View

(Color online). The centrality dependence of $R_{\langle N_{\text{ch}}\rangle }$ with spherical (Sph) and deformed (Def) colliding nuclei, calculated with (a) Optical Glauber (Opt.Glb) (b) Monte Carlo (MC.Glb) Glauber models

This is not the case for the anisotropic flow observables, where previous studies suggest that both nuclear deformation and event-by-event fluctuations are important [39, 49, 50]. Apart from nuclear deformation, most of the anisotropic flow produced in the most central heavy-ion collisions is due to the position fluctuations of nucleons in the colliding nuclei. In this paper, we study the effect of nuclear deformations and initial fluctuations on geometry eccentricities $R_{{ϵ}_n}$ (n=2,3). The results are shown in Fig. 4. Note that, for the Optical Glauber simulations, ϵ₃ vanishes for the spherical case due to symmetry, and for the deformed case the $R_{{ϵ}_3}<-1$ (because ${ϵ}_3^{\text{ZrZr}}\gg {ϵ}_3^{\text{RuRu}}$) is outside the frame of Fig. 4(b). Comparing the results between deformed and spherical nuclei, the general features are similar for both. This is consistent with the results shown in previous studies [6, 8], where the large quadrupole deformation β₂ in Ru increases the $R_{{ϵ}_2}$ in most central collisions while the large octupole deformation β₃ suppresses the values in mid-central collisions.

Fig. 4Full-screen View

(Color online). The centrality dependence of (a) $R_{{ϵ}_2}$ and (b) $R_{{ϵ}_3}$. The results are calculated with spherical (Sph) and deformed (Def) colliding nuclei using Optical Glauber (Opt.Glb) and Monte Carlo Glauber (MC.Glb) models. The results from Monte Carlo Glauber model with quark participant assumption (MC.qGlb) are presented as open symbols. From Optical Glauber simulations, ϵ₃ vanishes due to symmetry in spherical cases, and the $R_{{ϵ}_3}<-1$ in deformed cases are outside the frame of the plot

The magnitudes of $R_{{ϵ}_2}$ show significant differences between the Optical and Monte Carlo Glauber models. The $R_{{ϵ}_2}$ caused by the nuclear geometry is significantly diluted by the event-by-event fluctuations. The reason is rather straightforward: event-by-event fluctuations contribute to a large fraction of the eccentricities in Monte Carlo simulations, and these contributions are common between ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ collisions, suppressing the differences caused by nuclear deformation and/or collision geometry.

In the Monte Carlo Glauber simulation, the participating nucleons are usually treated as point-like particles as we did in the simulations above. However, nucleons have substructures [51-53]. More substructure constituents are expected to suppress the fluctuations. We use the quark participant assumption as it was done in Ref. [51], increasing the density by three times and decreasing the inelastic cross-section by a factor of 9. The results are also shown in Fig. 4. The Monte Carlo Glauber simulations under quark participant scenario significantly increase the $R_{{ϵ}_2}$ and decrease $R_{{ϵ}_3}$, bringing the values closer to the Optical limit. We note that despite of the cross-section differences, the $R_{{ϵ}_n}$ calculated in the Optical Glauber model are similar to those calculated in the Monte Carlo Glauber model with respect to reaction plane instead of participant plane [3], which is more related to the elliptic flow measurement with respect to zero-degree calorimeter (ZDC) instead of time projection chamber (TPC) in experiment [1].

There are more sources of initial fluctuations in heavy ion collisions, for example, the correlations between the quarks in a nucleon, as well as the correlations between the nucleons in the colliding nuclei. We find that the minimum distance between two nucleons in simulations, as one of the sources of nucleon-nucleon correlations, does not change our results. We note that the final-state effects, which are not the focus of the present study, are also important for the accurate extraction of the symmetry energy slope parameter. We postpone such investigation to a future work.

Summary

In this work, the effects of nuclear deformations and initial fluctuations on the ratio observables in relativistic isobar collisions are studied using the Optical Glauber and Monte Carlo Glauber models. The GPU parallel computing technology is used in our calculations to achieve high precision, especially for the deformed colliding nucleus. Our main findings are as follow.

• $R_{\langle N_{\text{ch}}\rangle }$ depends on particle production mechanisms and model parameters, while it converges to a single value in most central collisions in the Trento model (see Fig. 1).

• Initial fluctuations affect the total cross section and thus the $R_{\langle N_{\text{ch}}\rangle }$, however the effect is negligible for most central collisions (see Fig. 2).

• Because of the constraints of the total volumes and the measured RMS radii of the isobar nuclei, their deformations are found to have negligible effects on $R_{\langle N_{\text{ch}}\rangle }$ (see Fig. 3).

• The $R_{{ϵ}_2}$ and $R_{{ϵ}_3}$, driven by the nuclear diffusenesses and deformations, are suppressed by the initial finite number fluctuations, dependent on the degree of freedom of nucleons or constituent quarks (see Fig. 4).

Our results indicate unique (in)sensitivities of the ratio observables of $R_{\langle N_{\text{ch}}\rangle }$ and $R_{{ϵ}_n}$ to nuclear deformations and initial fluctuations in relativistic isobar collisions. These features, compared to data, may potentially probe the particle production mechanism and the physics underlying nuclear structure.