Methodology

Figure 1 illustrates the method for accessing the structure of ions colliding at relativistic energies. (A) Two nuclei are smashed in a high-energy collider (the large Lorentz contraction in the beam direction is not shown). (B) At the time of interaction, the nuclei are characterized by nontrivial geometries of nucleon configurations, including deformations and radial profiles. (C) The geometry of such configurations is reflected in the initial condition of the created QGP. The subsequent hydrodynamic expansion of this system, driven by pressure-gradient forces, converts the spatial asymmetries in the initial shape into the momentum asymmetries of emitted particles in the transverse plane. (D) Experimentally, transverse momentum (p_T) asymmetries can be revealed via a Fourier expansion of the particle distributions in azimuthal angle: $\frac{{\text{d}}^{2}N}{\text{d}{p}_{\text{T}}\text{d}\varphi }=\frac{\text{d}N}{\text{2πd}{p}_{\text{T}}}\left(1+2{\displaystyle \sum _{n=1}^{\infty }{v}_n\cos n\left(\varphi -{\text{Φ}}_n\right)}\right),$ (1) where the Fourier harmonics $V_n=v_n{e}^{in{\text{Φ}}_n}$ are coefficients of anisotropic flow. The most significant harmonics are V₂, elliptic flow, reflecting the elliptical asymmetry of the geometry of the QGP, and V₃, reflecting the triangular asymmetry [18-20]. We note that the total particle multiplicity, $N_{\text{ch}}={\displaystyle \int \text{d}{p}_{\text{T}}\text{d}N/\text{d}{p}_{\text{T}}}$, is proportional to the amount of energy deposited in a collision, which in turn is determined by the number of nucleons, N_part, participating in the interaction. The slope of the p_T spectra reflects the strength of the radial expansion, characterized by $\langle p_{\text{T}}\rangle \frac{1}{N_{\text{ch}}}\text{d}{p}_{\text{T}}{p}_{\text{T}}\text{d}N/\text{d}{p}_{\text{T}}$, which is inversely related to the transverse size of the overlap region [21]. Due to these relations, inherent to the hydrodynamic description, information about the structure of the colliding ions can be inferred from the detected final-state particles.

Fig. 1Full-screen View

(Color online) Schematic view of a relativistic heavy-ion collision, highlighting the role played by the collective properties of the colliding ions in shaping the geometry of the produced quark-gluon plasma (QGP). The Lorentz contraction of the two nuclei in the z-direction, by factor γ～100 at the BNL RHIC or over 1000 at the CERN LHC, is not shown. See text for a detailed description

The most direct way of observing the impact of the nuclear structure via this method is through comparing observables measured in collisions of species that are close in mass. Isobars, i.e., nuclides having the same mass number, are ideal candidates for such studies [16], as explicitly demonstrated by experimental data from ⁹⁶Zr+⁹⁶Zr and ⁹⁶Ru+⁹⁶Ru collisions, collected in 2018 at the BNL RHIC and released three years later [4]. Given two isobars, X and Y, and a given observable, $\mathcal{O}$ 𝒪, we ask the following question: $\boxed{\frac{{\mathcal{O}}_{X+X}}{{\mathcal{O}}_{Y+Y}}\stackrel{?}{=}1}$ (2) Model studies have established that any visible departure from unity in the ratio must originate from differences in the structure of the isobars. In the measurements of the STAR collaboration, structure influences are ubiquitously found. Ratios of more than ten observables taken between ⁹⁶Zr+⁹⁶Zr and ⁹⁶Ru+⁹⁶Ru have been measured, all displaying distinct and centrality-dependent deviations from unity, as reported in Fig. 2 [22]. Such rich and versatile information can provide a new type of constraint on the structure of these isobars, as also predicted by early model investigations, which we discuss below.

Fig. 2Full-screen View

(Color online) Ratios of observables taken between ⁹⁶Ru+⁹⁶Ru and ⁹⁶Zr+⁹⁶Zr collisions as a function of N_ch, as measured by the STAR Collaboration (Preliminary results). A total of ten ratios are shown

Nuclear structure input

The hydrodynamic model of heavy-ion collisions successfully reproduces a vast set of experimental measurements at the BNL RHIC and the CERN LHC [23]. The input to hydrodynamic simulations is the event-by-event distribution of nucleons in the colliding ions. Motivated by low-energy nuclear physics, a Woods-Saxon profile with a nuclear surface expanded in spherical harmonics is routinely employed, $p\left(r,\theta ,\varphi \right)\propto \frac{1}{1+e^{\left[r-R_0\left(1+{\beta }_2\left(\cos \gamma Y_2^0\left(\theta ,\varphi \right)+\sin \gamma Y_2^2\left(\theta ,\varphi \right)\right)+{\beta }_3{Y}_3^0\left(\theta ,\varphi \right)\right)\right]/a_0}},$ (3) where R₀ is the half-density radius, a₀ is the surface diffuseness, β₂ is the magnitude of the quadrupole deformation, γ determines the relative length of the three axes of the ellipsoid, and β₃ is the (axial) octupole deformation parameter, where β₃≠0 implies a breaking of parity symmetry in the intrinsic nuclear shape. Hydrodynamic simulations show that any changes in these parameters leave characteristic and detectable impacts on experimental observables such as those shown in Fig. 2 [24, 25].

Alternatively, hydrodynamic simulations can take pre-sampled nucleon configurations from ab initio calculations as input (See, e.g., Refs. [26-29] for such applications in ¹⁶O collisions). Here, diffuseness and deformations emerge directly from many-nucleon correlations in the sampled wave functions. Given the expected rapid progress in the reach and quality of ab initio calculations over the next few years [17, 30], this alternative approach should become broadly adopted in the modeling of heavy-ion collisions in future. Full exploitation of such predctions of state-of-the-art nuclear theory will demonstrate further the scientific relevance of the connection between high-energy observations and low-energy theories.

Signatures of intrinsic nuclear shapes

A crucial observable in high-energy heavy-ion collisions is the rms flow coefficient, $v_n=\sqrt{\langle {\left|V_n\right|}^2\rangle }$. Numerical and semi-analytical studies show that, for collisions at a given multiplicity (or centrality), vn is enhanced by the presence of nuclear deformations in the colliding ions, following[16, 31-33], $v_n^2\approx b_0+b_1{\beta }_n^2,$ (4) where b₀ and b₁ are positive coefficients that depend on centrality. The enhancement predicted by Eq. (4) would show up, in particular, when comparing collisions of deformed nuclei to collisions of spherical nuclei. A powerful way to do so is to compare deformed ion collisions with collisions of nearly spherical ²⁰⁸Pb ions. The left panel of Fig. 3 reveals an enhanced v₂ in ¹²⁹Xe+¹²⁹Xe collisions compared to ²⁰⁸Pb+²⁰⁸Pb collisions [34], as observed by the ALICE collaboration [35]. A state-of-the-art calculation [36] confirms the origin of this effect due to the large β₂ of ¹²⁹Xe.

Fig. 3Full-screen View

(Color online) Modification of multi-particle correlation observables in Xe+Xe collisions compared to the baseline with spherical nuclei provided by Pb+Pb collisions. Left: elliptic flow, v₂ [35]. Right: correlation between elliptic flow and the average transverse momentum, ρ₂ [45]

Concerning the triaxiality, γ in Eq. (1), revealing its presence requires the use of three-particle correlations. The most sensitive observable is the correlation of the shape of the QGP with its size [37], measurable experimentally via a correlation between $v_n^2$ and the fluctuation of transverse momentum, $\delta p_{\text{T}}=p_{\text{T}}-\langle p_{\text{T}}\rangle$, at a given multiplicity. This quantity is conveniently formulated as a Pearson coefficient [38], ${\rho }_n=\frac{\langle v_n^2\delta p_{\text{T}}\rangle }{\sqrt{\left(\langle v_n^4\rangle -{\langle v_n^2\rangle }^2\right)}\sqrt{\delta p_{\text{T}}\delta p_{\text{T}}}}$. For quadrupole deformation, theoretical work shows the following leading-order dependence [39] ${\rho }_2\approx {b^{\prime }}_0-{b^{\prime }}_1{\beta }_2^3\cos \left(3\gamma \right),$ (5) where ${b^{\prime }}_0$ and ${b^{\prime }}_1$ are positive coefficients. In the presence of large β₂, moving from oblate (γ=60°) to prolate (γ=0) shapes decreases ρ₂ in a substantial way. A recent measurement at RHIC shows precisely ρ_2<0 in central U+U collisions [40], which is explained naturally by the large prolate deformation of ²³⁸U [41], β₂～0.28, γ=0. The nucleus ¹²⁹Xe is particularly interesting for such a study, as its shape is considerably deformed and also triaxial, β₂=0.2 and γ≈30° [42-44]. In the right panel of Fig. 3, model calculations assuming oblate, triaxial, and prolate ¹²⁹Xe shape show a strong modification of ρ₂ in ¹²⁹Xe+¹²⁹Xe collisions with respect to the ²⁰⁸Pb+²⁰⁸Pb collisions [43]. Measurements from the ATLAS collaboration indeed confirm the triaxial scenario [45]. One important point is that the combined use of and ρ₂ can simultaneously constrain β₂ and γ.

In the octupole sector, much less is known from low-energy physics [46]. Direct evidence of octupole deformation in excitation bands of atomic nuclei is scarce, because octupole deformation rarely manifests as a mean-field effect (static deformations) [47, 48], as in a simple rotor model. However, dynamical octupole correlations at the beyond-mean-field level are present in essentially all nuclei [49], and should leave their fingerprint in the nucleon configurations from ab initio calculations. High-energy nuclear collisions, probing configurations of nucleons on an event-by-event basis, give access to all such non-static deformations in the ground states in the same way as the static ones.

One of the breakthrough outcomes of the isobar collision campaign at RHIC is reported in Fig. 4, also shown in the mid-bottom panel of Fig. 2. The ratio of vn taken between Ru+Ru and Zr+Zr collisions shows significant departures from unity. The data implies that ⁹⁶Ru has a larger β₂ than ⁹⁶Zr, as expected from low-energy experiments. A similar departure for n=3, showing an enhanced v₃ in Zr+Zr collisions, can only be ascribed to ⁹⁶Zr having a sizable β₃ [50], which is not predicted by mean-field energy density functional calculations [48, 51]. The results of the STAR collaboration demonstrate that heavy-ion collisions offer a clean access route to multi-nucleon correlations that are both difficult to quantify from traditional low-energy experiments and hard to predict from phenomenological models.

Fig. 4Full-screen View

(Color online) Preliminary ratios of flow coefficients, vn, taken between Ru+Ru and Zr+Zr collisions. The suppression of the v₃ ratio at large multiplicity is due to an enhancement of v₃ in Zr+Zr collisions

Radial profiles and relation to neutron distributions

The nuclear radial profile, determined by the R₀ and a₀ parameters in Eq. (2), influences the area and the density of the overlap region. In general, a smaller a₀ or R₀ for a fixed mass number leads to a sharper edge in the overlap geometry, leading to a more compact QGP, larger pressure gradients, and hence larger 〈p_T〉 and vn. The impact is more significant in off-central collisions where the overlap region is smaller, and sensitivity to a variation in R₀ and a₀ is larger. Indeed, model studies show that the probability distributions of N_part, and hence the distribution of N_ch, p(N_ch), as well as 〈p_T〉 and v₂, are largely impacted by variations in a₀ and R₀ [52, 53, 24].

Due to model-dependent systematics, constraining the radial nuclear profile in a single collision system is difficult. Such limitation is largely overcome by comparing experimental observables between systems close in size, such as isobars. Assuming the differences of radial parameters are small, deviation of isobar ratios from unity can be approximated by (taking ⁹⁶Ru and ⁹⁶Zr as an example) $\frac{{\mathcal{O}}_{\text{Ru}}}{{\mathcal{O}}_{\text{Zr}}}\approx 1+c_0\left(R_{0,\text{Ru}}-R_{0,\text{Zr}}\right)+c_1\left(a_{0,\text{Ru}}-a_{0,\text{Zr}}\right),\mathcal{O}\equiv p\left(N_{\text{ch}}\right),v_2,\text{ or }\langle p_{\text{T}}\rangle ,$ (6) where the coefficients c₀ and c₁ depend only on the mass number at a given centrality or multiplicity and are insensitive to the final state effects [54]. These simple equations describe well the isobar ratios, as verified in recent transport model simulations [24].

Energy density functional calculations suggest that ⁹⁶Zr has a larger diffuseness but a smaller radius than ⁹⁶Ru, i.e. $\text{Δ}{a}_0\equiv a_{0,\text{Ru}}-a_{0,\text{Zr}}<0$ and $\text{Δ}{R}_0\equiv R_{0,\text{Ru}}-R_{0,\text{Zr}}>0$ [56, 57]. As shown in Fig. 5, a transport calculation implementing such differences can quantitatively describe the measured ratios of v₂｛4｝ (the fourth-order cumulant of elliptic flow fluctuations) and of the distribution of charged particle number, p(N_ch). The cumulant v₂｛4｝ measures the flow originating from the intrinsic ellipticity acquired by the QGP due to the finite impact parameter of the collisions [58]. Figure 5 shows that such intrinsic ellipticity is impacted only by Δa₀ [59], but is insensitive to ΔR₀ and nuclear deformations. On the other hand, p(N_ch)-ratio is sensitive to both Δa₀ and ΔR₀ [56, 57]. A dedicated study finds that the ratio of p_T is also sensitive to both Δa₀ and ΔR₀ [53]. Therefore, the measurement of isobar ratios provides several independent determinations on the differences ΔR₀ and Δa₀, which can be confronted against the predictions of low-energy nuclear structure models.

Fig. 5Full-screen View

(Color online) Ratios of observables taken between Ru+Ru and Zr+Zr collisions. The inset in the left panel shows how the neutron excess of ⁹⁶Zr compared to ⁹⁶Ru yields a more diffuse nuclear surface, i.e., a larger a₀ and a slightly smaller R₀ in Eq. (1). Left: the impact of the larger a₀ of ⁹⁶Zr manifests predominantly in the fourth-order cumulant of v₂, v₂｛4｝, which originates from the intrinsic ellipticity of the QGP due to the finite collision impact parameter. Right: both a₀ and R₀ differences contribute to a broad hump in the ratio of p(N_ch). Calculations are obtained within the A Multi Phase Transport (AMPT) model [55]

The knowledge of nucleon distribution, in combination with the well-known proton distribution parameters a_p and R_0p from low-energy experiments, allows one to probe the difference between the rms radius of neutrons and protons in heavy nuclei, Δr_np= R_n-R_p, known as the neutron skin. The value of Δr_np is directly related to the slope of the symmetry energy, dubbed L, appearing in the equation of state (EOS) of nuclear matter [60]. Determinations of L are intensively pursued at low energy because this parameter plays a crucial role in the stability properties of neutron stars [61, 62]. Isobar ratios in high-energy collisions are expected to probe only the difference in the neutron skin, Δ(Δr_np)=Δr_np,Ru-Δr_np,Zr. Assuming Woods-Saxon distributions for protons and nucleons, Δr_np receives a contribution from both half-radius and surface diffuseness [24]: $\text{Δ}\left(\text{Δ}{r}_{\text{np}}\right)\propto \left(R_0\text{Δ}{R}_0-R_{0\text{p}}\text{Δ}{R}_{0\text{p}}\right)+7/3{\pi }^2\left(a\text{Δ}a-a_{\text{p}}\text{Δ}{a}_{\text{p}}\right)$. Therefore, collisions of isobars or, in general, of species of similar mass numbers allow one to access detailed information about radial profiles and neutron skins of nuclei systematically.

Stress-testing small system collectivity with ²⁰Ne

The neon-20 nucleus presents the most extreme ground state of all stable nuclides with A>10. It is a strongly-deformed object made of five α-clusters in a reflection-asymmetric α+¹⁶O molecular configuration [63-65]. In terms of the common quadrupole deformation coefficient, the ground state has β₂≈0.7, the highest of all stable ground states. The deformation of this nucleus is so large that its impacts can easily survive the large event-by-event fluctuations associated with sampling a small number of nucleons ($\propto 1/\sqrt{A}$). The extreme geometry of ²⁰Ne enables us to perform nontrivial tests of the initial-state modeling and the hydrodynamic response in small systems. In particular, one can compare ²⁰Ne+²⁰Ne collisions with collisions of nearly-spherical ¹⁶O nuclei (already recorded at RHIC, and are also planned for 2024 at the LHC). The ratios of observables between the two systems will be largely independent of final state transport properties, and hence directly access the variation in the initial condition caused by nuclear structure differences. Having data from ²⁰Ne+²⁰Ne collisions will maximize the scientific output of the ¹⁶O+¹⁶O run (and vice versa). On the side of nuclear structure theory, ab initio approaches have recently been pushed to describe light systems up to A～40[17], including ²⁰Ne [66]. Strong deformations in these approaches emerge from genuine n-body (up to A-body) correlations in the wavefunction generated by inter-nucleon interactions linked to QCD via an effective field theory. Precise measurement of multi-particle correlations in ²⁰Ne+²⁰Ne collisions will provide novel tests of the effectiveness of such ab initio calculations in capturing collective effects in strongly-correlated nuclei.

As a bonus, while collecting ²⁰Ne+²⁰Ne collisions in collider mode at 7 TeV, one can have the same collisions in fixed-target mode at around 0.07 TeV by injecting a ²⁰Ne gas in the SMOG system of the LHCb experiment [67]. This would enable a study of the $\sqrt{s_{\text{NN}}}$ dependence of the initial condition, longitudinal dynamics and geometry of small systems. Another, potentially superior way of imaging the structure of light nuclei is to collide them with heavy spherical nuclei, such as in ¹⁶O+²⁰⁸Pb or ²⁰Ne+²⁰⁸Pb [68, 69]. The shape of overlap region at small impact parameter directly captures the nucleon distribution in light nuclei. Ratios of observables between ¹⁶O+²⁰⁸Pb and ²⁰Ne+²⁰⁸Pb collisions will reveal the shape differences between ¹⁶O and ²⁰Ne. The main advantage over symmetric ¹⁶O+¹⁶O and ²⁰Ne+²⁰Ne collisions is that asymmetric “isobar”+Pb collisions produce much more particles and will have a better centrality resolution [67]. This idea may be extended to even smaller systems such as ⁸Be+²⁰⁸Pb and ¹²C+²⁰⁸Pb [68, 69]. Collisions of asymmetric systems is feasible at the BNL RHIC (as demonstrated by the p+Au, d+Au and ³He+Au runs [70]) but not at the CERN LHC in collider mode. However, asymmetric collisions can be performed in fixed-target mode by injecting oxygen-16 and neon-20 ions in the SMOG system of the LHCb detector [71] (a small sample of Ne+Pb data was collected in 2013).

Shape evolution along the Samarium isotopic chain

Certain isotopic chains in the nuclear chart exhibit strong variations in nuclear shapes. While this occurs mainly away from the stability line, the chain of eight stable samarium isotopes (Sm Z=62) features a transition from nearly-spherical to strongly-deformed nuclei with increasing neutron number, e.g. from ¹⁴⁴Sm with β₂≈0.09, to ¹⁵⁴Sm with β₂≈0.34, with a change in mass number of only about 7%. Since the hydrodynamic response is expected to be essentially constant over the isotopic chain, these systems offer a strong lever-arm to probe in detail how the initial condition of QGP responds to varying nuclear shapes, e.g. by predicting the coefficients b₀ and b₁ in Eq. (4) using two Sm isotopes and then make predictions of β₂ for other isotopes [33]. The β₂ differences among isotopes can be extracted from ratios of flow observables, as done for the BNL RHIC isobar run. The extracted differences from heavy-ion collisions can be compared with nuclear structure knowledge, to study whether shapes evolve similarly when adding neutrons one-by-one in low-energy experiments and high-energy collisions. We stress that these nuclei have been subject of much investigation at low energy, where their properties are nicely consistent across experiments and theoretical frameworks.

It is worth noting, then, that scanning the Sm isotopic chain in high-energy collisions would provide new experimental insight onto the octupole deformations of such nuclei. As demonstrated by the isobar ratios, nontrivial results are expected. Clear observation of octupole, and potentially hexadecapol deformations for such nuclei would showcase the discovery potential of high-energy nuclear collisions as a tool to observe the manifestations of many-body correlations of nucleons in the ground state of nuclei, in a way that is fully complementary to low-energy structure experiments. In turn, this will provide new experimental constraints to test future ab initio calculations of such large and deformed systems.

The neutron skin of ⁴⁸Ca and ²⁰⁸Pb in high-energy collisions

In low-energy experiments, the neutron skins of ⁴⁸Ca and ²⁰⁸Pb, two doubly-magic nuclei with a considerable neutron excess, have been the subject of much work. Dedicated experiments at Jefferson Lab have been devoted to measuring the neutron skin of these species [72, 73]. The measured value for ²⁰⁸Pb is Δr_np=0.28 ± 0.07 fm, which is systematically larger than predictions from energy density functional theories. The properties of neutron stars (e.g., the tidal deformability) resulting from such a constraint on the EOS turn out to be slightly at variance with those inferred from pulsar and gravitational wave observations, which has sparked intense debate in the community [74, 75]. The neutron skin of ⁴⁸Ca is instead more in line with the theoretical expectations. We aim to provide new constraints on the neutron skins of ⁴⁸Ca and ²⁰⁸Pb by utilizing high-energy collisions.

Providing a robust estimate of the neutron skin of ⁴⁸Ca in high-energy nuclear collisions is rather straightforward. The isotopic chain of calcium has two doubly-magic nuclei, ⁴⁸Ca and ⁴⁰Ca. The latter has the same number of protons and neutrons, and its neutron skin is much smaller than that of ⁴⁸Ca. However, experiments reveal that ⁴⁸Ca and ⁴⁰Ca have essentially the same charge radius with a difference less than 0.001 fm [76, 77], such that neutrons alone determine the differences in size between these two isotopes. As discussed in Sect. 2.4, heavy-ion collisions allow one to experimentally access differences in the neutron skins between nuclei of similar mass. Therefore, if Δr_np(⁴⁸Ca)≫Δr_np(⁴⁰Ca)≈0, collisions of such nuclei could isolate $\text{Δ}{r}_{\text{np}}\left({}^{48}\text{Ca}\right)-\text{Δ}{r}_{\text{np}}\left({}^{40}\text{Ca}\right)\simeq \text{Δ}{r}_{\text{np}}\left({}^{40}\text{Ca}\right).$ (7) We estimate that this quantity can be accessed with an uncertainty of about 0.02 fm. Any significant deviations from the expectations of low-energy theories or experiments should be ascribed to the modification of the partonic structure of nucleons in nuclear environment at high energy.

Following this idea, the constraints on neutron skin of ²⁰⁸Pb could be obtained by comparing data from ²⁰⁸Pb+²⁰⁸Pb with data from ¹⁹⁷Au+¹⁹⁷Au, as the two species are nearly isobars. Therefore, having such collisions at the same beam energy would allow us to determine the difference Δr_np.Pb - Δr_np.Au from observables such as v₂｛4｝. This information could be combined with an additional estimate of the neutron skin from a method recently developed by the STAR collaboration [78], also at high energy. This method employs the production of ρ⁰ mesons in photo-nuclear processes in ultra-peripheral collisions using the newly developed spin interference enabled nuclear tomography. The cross section for ρ⁰ production in dipole-nucleus scattering contains a coherent component determined by the gluon distribution of the target nucleus. Fits of the coherent diffractive |t| distribution within a Woods-Saxon geometry model in ¹⁹⁷Au+¹⁹⁷Au collisions lead to Δ_np(¹⁹⁷Au)=0.17±0.03(stat.)±0.08(sys) fm. This method could be readily applied to other species such as ²⁰⁸Pb via ²⁰⁸Pb+²⁰⁸Pb collisions. It would measure the neutron skin of ²⁰⁸Pb with an uncertainty that is similar to or even better than that obtained by the PREX-II experiment. We emphasize that the systematic errors are largely correlated in this technique. The experiment should be able to demonstrate whether the extracted neutron skin difference between ²⁰⁸Pb and ¹⁹⁷Au is compatible with low energy models and measurements (including PREX-II for Pb [72]). We note that a short Pb+Pb collision run at RHIC would be sufficient for this purpose. This is a cost-effective experiment with significant impacts on the nuclear physics community as a whole.

Furthermore, it is worth noting that at the energy reached at the LHC, electro-weak (EW) bosons are abundantly produced in nucleus-nucleus collisions via $u\overline{d}\to W^{+}$, $d\overline{u}\to W^{+}$, and $q\overline{q}\to Z$ processes. These real EW bosons probe the weak charge distributions in the heavy-ion initial state, i.e., the sum of weak charge distributions from the two colliding nuclei in the overlap region. Therefore, isobar ratios of W and Z boson yields as a function of centrality may provide direct access to the radial distribution of valence and sea quarks, offering an access route to charge radius, mass radius and thus the neutron skin, and compare with that extracted from the PREX-II analysis of the neutral weak form factor of ²⁰⁸Pb [79] associated with virtual Z boson. Based on 1 nb^-1 5.02 TeV Pb+Pb data (from one month of typical heavy-ion running), we expect each LHC experiment to deliver about 600k W bosons and 20k Z bosons reconstructed in the lepton decay channel [80, 81]. If isobar or isobar-like collisions, such as ⁴⁰Ca+⁴⁰Ca vs. ⁴⁸Ca+⁴⁸Ca, or ²⁰⁸Pb+²⁰⁸Pb vs. ¹⁹⁷Au+¹⁹⁷Au become available at LHC, it will be possible to compute ratios of W and Z boson yields as a function of centrality with a statistical uncertainty of order 1% ($\sqrt{1/10000}$). This could be achieved in particular at the high-luminosity LHC, in Run5 or beyond (<2035) [82, 83].

Initial conditions of heavy-ion collisions

The success of the hydrodynamic framework of heavy-ion collisions enables us today to perform quantitative extractions of the transport properties of the QGP via multi-system Bayesian analyses [11-15]. A major limitation of such extractions is the lack of precise control on the initial condition of the QGP prior to the hydrodynamic expansion. Insights about the energy deposition from two collided nuclei come from the color glass condensate (CGC) effective theory of high-energy QCD [84]. There, for a given boosted nuclear profile described by the thickness function $T={\displaystyle \int \rho \left(x,y,z\right)\text{d}z}$, the average energy density deposited in the transverse plane in the collision of, say, nuclei A and B, at the instant immediately after the collision occurs is of the form [85] $\langle T^{00}\rangle \left[{\text{GeV/fm}}^3\right]\propto T_A{T}_B.$ (8) Bayesian analyses of heavy-ion data, while constraining transport properties of the QGP, attempt as well to constrain the initial conditions of the collisions. The prediction of the CGC in Eq. (8) can be tested via a generic parameterized Ansatz for the energy density, such as the TENTo Ansatz for the energy density per unit rapidity, dE/dy[GeV/fm²], namely $\text{d}E/\text{d}y\propto {\left({\widehat{T}}_A^p+{\widehat{T}}_B^p\right)}^{1/p}$ [86], or its generalized version, $\text{d}E/\text{d}y\propto {\left({\widehat{T}}_A^p+{\widehat{T}}_B^p\right)}^{q/p}$ [23, 87], where $\widehat{T}$ represents the thickness function constructed solely from the participant nucleons within the colliding ions. The values of p and q and other model parameters such as nucleon width w and inter-nucleon minimum distance d_mincan then be inferred from the analysis of high-energy collision data.

However, information about the content of the colliding nuclei, which can not be predicted based on the CGC alone, yields a significant uncertainty in our understanding of the energy deposition itself and, in turn, of the QGP transport parameters resulting from fits to data. One example is provided by v₃ in Zr+Zr collisions. If one attempted to reproduce the measured v₃ in hydrodynamic calculations without implementing any β₃ parameter for such a nucleus, one would correct a 10% enhancement of such an observable in central collisions by biasing the extraction of other QGP transport or initial-state properties dramatically. It is the knowledge of the presence of a large octupole deformation from the isobar ratio v_3,Zr/v_3,Ru, that enable us to avoid biasing the extracted QGP features. Bayesian approaches have not yet systematically explored the impact of nuclear shape and radial distributions. Nuclear structure knowledge should be used systematically as a new lever arm to probe the initial condition of collisions of species that are close in mass, and thus obtain better determinations of the QGP transport coefficients.

As discussed in Sect. 2.1, the deviation of isobar ratios from unity probes directly the structural differences between the two species, and the way the initial condition is shaped by two colliding ions. Numerical work shows that the ratios of many observables can be expressed in terms of the differences of Woods-Saxon parameters as a generalization of Eq. (6) [24], $\begin{array}{c}\frac{{\mathcal{O}}_{\text{Ru}}}{{\mathcal{O}}_{\text{Zr}}}\approx 1+c_0\left(R_{0,\text{Ru}}-R_{0,\text{Zr}}\right)\\ +c_1\left(a_{0,\text{Ru}}-a_{0,\text{Zr}}\right)+c_2\left({\beta }_{2,\text{Ru}}^2-{\beta }_{2,\text{Zr}}^2\right)+c_3\left({\beta }_{3,\text{Ru}}^2-{\beta }_{3,\text{Zr}}^2\right),\end{array}$ (9) with $\mathcal{O}\equiv p\left(N_{\text{ch}}\right)$, v₂ or 〈p_T〉. Crucially, these ratios are insensitive to variations of QGP transport properties [54]. Therefore, the left-hand side of Eq. (9) captures the variations of initial conditions in the isobar systems, which are related to the structure parameters on the right-hand side. The coefficients cn reflect how the initial condition changes when the nuclear structure is varied between the isobars. Via isobar collisions, thus, one can conveniently separate the role played by low-energy nuclear structure input (R₀, a₀, ${\beta }_2^2$, ${\beta }_3^2$) and role played by constraints from the knowledge of high-energy heavy-ion collisions (ci coefficients) for observables (isobar ratios) that depend solely on the initial condition of the QGP, as illustrated in Fig. 6.

Fig. 6Full-screen View

(Color online) Impact of isobar-like collisions on the goal of the heavy-ion program. Better control on the initial condition can be achieved by exploiting the constraints from both the ratios of final-state observables and the nuclear structure knowledge

Collisions from pp, p+A and A+A have been collected at the BNL RHIC and the CERN LHC. For A+A collisions, we have ²³⁸U+²³⁸U and ¹⁹⁷Au+¹⁹⁷Au collisions at the BNL RHIC, and ¹²⁹Xe+¹²⁹Xe and ²⁰⁸Pb+²⁰⁸Pb collisions at the CERN LHC. However, none of these pairs are close enough in their mass number¹, which means that the final-state effects do not completely cancel in the ratios [34]. The model dependencies of these residual effects are then significant enough to preclude a precise extraction of the initial condition. Colliding isobars, or in general species close in mass, such as ¹⁹⁷Au and ²⁰⁸Pb, represents an ideal way to constrain the initial condition across the nuclide chart.

Concerning the possibility of having Pb+Pb reference data at the BNL RHIC, there are two arguments in the context of hot QCD studies to motivate such an effort, in addition to the neutron skin case pointed out in Sect. 2.3. 1) Being doubly magic, ²⁰⁸Pb is essentially spherical. In contrast, ¹⁹⁷Au has a modest oblate deformation. For the high-precision studies of Au+Au collisions expected from the upcoming sPHENIX program, it would be important, then, to have Pb+Pb collisions as a tool to calibrate the initial condition of Au+Au collisions and ensure that the expectations of the low-energy nuclear theory are compatible with the observations at high energy. 2) Having Pb+Pb systems would also provide a bridge to compare the outcome of Pb+Pb collisions at the BNL RHIC to that of Pb+Pb collisions at the CERN LHC, to study the beam energy dependence of observables. For both these goals, a short Pb+Pb run at the BNL RHIC would be sufficient.

Last but not least, isobar or isobar-like collisions may serve as novel probes of the hard sector, via the analysis of observables such as the production of leading hadron, jets, photons, and heavy flavors. It has already been shown that collective flow of D-mesons is sensitive to the deformation of the nucleus [88]. Additionally, by constructing ratios of selected observables at a given centrality, or multiplicity, final state effects such as jet quenching are expected to cancel. Deviations from unity in the constructed ratios will provide access to flavor-dependent Nuclear Parton Distribution Function (nPDF), tailored for each underlying hard-scattering process. Interestingly, the precision determination of the impact parameter from bulk particles in coincidence with hard processes means that we can use isobar ratios to detect differences in the transverse spatial distribution of partons at given longitudinal momentum fraction between the two isobars.

One such example is already discussed in Sect. 3.3 in the context of W and Z bosons for neutron skin measurements. For some of the hard probes, such as high-p_T charged hadrons or inclusive jets, production will be so abundant that even a short run would permit one to determine the isobar ratio with a statistical precision of 1% or better as a function of centrality, with large cancellation of the systematical uncertainties. One could also study how the isobar ratio evolves with rapidity to detect potential modifications to the nuclear structure inputs due to nPDF or gluon saturation. Isobar ratios of more differential measurements, such as dijet or photon-jet measurements, could probe in more detail the correlation between the final-state medium effects, such as quenching, and the geometry of the hard-scattering processes, such as the path length. For this purpose, collisions of different species should be taken at the same $\sqrt{s_{\text{NN}}}$, with similar pileup and detector conditions. Model studies are forthcoming to put this argument on a more quantitative ground.

Impact on future experiments: EIC and CBM FAIR

Collisions of isobars may provide valuable input to the physics of the planned EIC. One important goal of the EIC program is to understand the partonic structure of nuclei at very high energy [89, 90]. At small longitudinal momentum fraction, x, the density of gluons may saturate and form the so-called color glass condensate (CGC). EIC will probe gluon saturation using a range of scattering processes in electron-nucleus collisions. In heavy-ion collisions, the modification of parton distributions in nuclei (nPDF) impacts the initial conditions of the QGP, which in turn are imaged via the isobar ratios of bulk and high-p_T observables. In this way, one can gain access to the transverse spatial distribution of partons. Exploiting isobar ratios as a function of rapidity, and in particular as a function of $\sqrt{s_{\text{NN}}}$, may provide a unique probe of x-dependence of nPDF and gluon saturation. Collisions of the same isobar pair, for example Ru+Ru and Zr+Zr, at different energies could be realized at CBM FAIR [91] at $\sqrt{s_{\text{NN}}}$≤4.9 GeV, BNL RHIC at $\sqrt{s_{\text{NN}}}$=200 GeV, and the CERN LHC at $\sqrt{s_{\text{NN}}}$>5 TeV. Any differences between RHIC and LHC in the isobar ratios for soft and hard probes could be used to infer the energy dependence of initial conditions and in turn that of the partonic structure within nuclei. This study will complement the EIC program by providing additional information on the spatial structure of dense gluonic matter. In turn, this information will provide valuable input for the CBM experiment at FAIR in the study of the QCD phase diagram at low temperature and high baryon density [91], in particular, to inform theoretical models such as SMASH [92], AMPT [93], and hydrodynamics [94-98] that aim to describe the dynamical and transport properties of nuclear matter in such conditions.