Determination of the neutron skin of 208 Pb from ultrarelativistic nuclear collisions

Emergent bulk properties of matter governed by the strong nuclear force give rise to physical phenomena across vastly diﬀerent scales, ranging from the shape of atomic nuclei to the masses and radii of neutron stars. They can be accessed on Earth by measuring the spatial extent of the outer skin made of neutrons that characterises the surface of heavy nuclei. The isotope 208 Pb, owing to its simple structure and neutron excess, has been in this context the target of many dedicated eﬀorts. Here, we determine the neutron skin from measurements of particle distributions and their collective ﬂow in 208 Pb+ 208 Pb collisions at ultrarelativistic energy performed at the Large Hadron Collider, which are sensitive to the overall size of the colliding 208 Pb ions. By means of state-of-the-art global analysis tools within the hydrodynamic model of heavy-ion collisions, we infer a neutron skin ∆ r np = 0 . 217 ± 0 . 058 fm, consistent with nuclear theory predictions, and competitive in accuracy with a recent determination from parity-violating asymmetries in polarised electron scattering. We establish thus a new experimental method to systematically measure neutron distributions in the ground state of atomic nuclei.

Understanding the distribution of neutrons within heavy atomic nuclei has profound implications for our knowledge of the neutron-rich matter that shapes exotic astrophysical objects such as neutron stars.The neutron skin that forms on the surface of heavy nuclei, whereby neutrons are located more diffusely and more on the outside [1,2], represents in particular a sensitive probe of the equation of state (EOS) of neutron matter, whose pressure determines the spatial extent of the neutron distributions.Indeed, nuclear models predict a strong correlation between the neutron skins of heavy nuclei and the masses and radii of neutron stars [3,4].
While proton distributions in nuclei can be determined in a model-independent way from electron scattering experiments [5], accessing neutron distributions poses a far greater challenge.As a consequence, we have only limited experimental constraints on the neutron skin of nuclei, ∆r np , defined as the difference in root mean square (rms) radii between protons and neutrons.The doublymagic nucleus 208 Pb (Z = 82, N = 126) has both protons and neutrons filling up their respective shells and represents an optimal study subject in this context.A recent, precise deduction of the neutron skin of 208 Pb has been achieved by the PREX collaboration [6] from the measurement of parity-violating asymmetries in polarised electron scattering.On the side of theory, the first calculation of 208 Pb and its neutron skin in the context of ab initio nuclear theory was also recently performed [7].These results, along with information coming from pulsar and gravitational wave observations, portray a picture of nuclear matter that hints at potential tensions [8][9][10].
In this article, we determine the neutron skin of 208 Pb from a new type of probe.We use data collected in 208 Pb+ 208 Pb collisions performed at ultrarelativistic energy at the CERN Large Hadron Collider (LHC).These collisions produce short-lived quark-gluon plasma [11][12][13] (QGP), the hot phase of quantum chromodynamics (QCD), which behaves like a near-ideal relativistic fluid [14,15] before fragmenting into observable particles.In high-energy scattering, interactions are mediated by gluons, such that the combined distribution of protons and neutrons (altogether called nucleons) within the colliding 208 Pb ions determines the shape and the size of the created QGP.Employing the latest advances in simulation and Bayesian inference tools within the hydrodynamic framework of heavy-ion collisions we reconstruct the geometry of the QGP by using the detected particle distributions.In conjunction with the precise knowledge of the proton density this enables us to place a tight constraint on the neutron skin of 208 Pb.
The neutron skin and the quark-gluon plasma -Our understanding of the QGP formed in 208 Pb+ 208 Pb collisions is highly developed thanks to the wealth of experimental data collected in the past decade by all LHC experiments, and in particular by the ALICE experiment dedicated to nuclear physics [16].Following Fig. 1, in an ultrarelativistic heavy-ion collision in the lab frame (Fig. 1a), interactions of gluons deposit energy density in the area of overlap in the so-called transverse plane, perpendicular to the beam direction (Fig. 1b).The deposition of energy density depends on the collision's impact parameter b, on the structure of the colliding nuclei and on the dynamics of the interaction itself.
Phenomenological studies have established a picture where the colliding ions are treated, in each collision (or event), as a superposition of nucleons that participate in the interaction.Both boosted nuclei are thus associated with a profile of matter in the transverse plane, T L,R (x ⊥ ), given as the sum of their participant nucleon profiles, typically taken as Gaussians with a width denoted by w.The interaction process and the subsequent energy depositions are then parameterised following some flexible prescription which can be fine-tuned directly from experimental data.Here we use a T R ENTo-type Ansatz for the energy density of the QGP [17,18], where L, R denote the two colliding ions, while p and q are model parameters.As the positions of the participant nucleons shaping the functions T L,R are sampled in each collision from the neutron and proton densities in the ground state of the scattering ions, the energy density e(x ⊥ ) is sensitive to their spatial distribution.This can be seen by eye in the density plot of Fig. 1b, representing an average energy density over many collisions.The scenario where the colliding 208 Pb nuclei have a narrower neutron skin leads to a QGP with a sharper profile over the plane and a higher density peak.
Starting from the initial condition discussed in Fig. 1b, the QGP then evolves as a relativistic viscous fluid (with transport properties, such as shear and bulk viscosities, that are also model parameters).For a single event, snapshots of the hydrodynamic expansion obtained using our hydrodynamic code are depicted in Fig. 1c.Cooling of the QGP lasts until the confinement crossover is reached, after which at a fixed switching temperature the fluid is converted into a gas of QCD resonance states that can further re-scatter or decay to stable particles.Out of this process, experiments can only detect final event-by-event stable particle spectra, typically denoted by: where p T is the transverse momentum, η is the particle pseudorapidity (η ≡ − ln tan(θ/2) with θ the polar angle in the (x ⊥ , z) plane of Fig. 1a), and the subscript ch indicates that only charged particles are included.We have conveniently decoupled the spectrum into a distribution of transverse momenta, p T ≡ |p T |, which quantifies the explosiveness of the QGP expansion, and an azimuthal component developed in Fourier modes, where v n are the so-called anisotropic flow coefficients that quantify the anisotropy of the particle emission.
Experimentally the first step is to sort the collisions in centrality classes based on the number of particles that they produce, where 0% centrality corresponds to events with the highest number of particles at almost zero impact parameter.As a function of centrality one can then measure among others the distributions of p T and v n coefficients for different particle species (pions, kaons, protons and more).This generates a wealth of experimental information from which the hydrodynamic model parameters (here, we have 26 in total) can be inferred.The central idea of this manuscript is that of promoting the neutron skin of 208 Pb to a model parameter that we constrain from LHC data.2. Signature of the neutron skin on bulk particle production in ultrarelativisitic 208 Pb+ 208 Pb collisions.Varying only the neutron skin size at our optimal parameter settings we show the charged particle multiplicity (left), the mean transverse momentum (middle) and the elliptic flow as measured by v2{k} (right) with a comparison to ALICE data [19,20].A larger neutron skin leads to more collisions, but per collision the multiplicity is lower at larger centralities.The larger size of the QGP leads to a reduced transverse momentum on average.Smearing of the elliptical shape leads to reduced elliptic flow as measured by v2{2} and v2{4}.Theoretical error bars are statistical only, experimental uncertainties include systematics as well.
The neutron skin is introduced by considering variations in the neutron diffuseness, a n , in the two-parameter Fermi distributions that model the point-neutron and point-proton densities in the colliding 208 Pb nuclei: We take a p = 0.448 fm, R p = 6.680 fm (corresponding to an rms proton radius r p = 5.436 fm), and R n = 6.690 fm, which is motivated by the experimental result that the neutron skin is caused by a more diffuse profile rather than a larger half-width radius [1,2].Before proceeding with a full Bayesian analysis we simulate the QGP formation and evolution for three different values of ∆r np while keeping all other model parameters fixed.First, a larger neutron skin leads to a larger total hadronic cross section, σ tot (see Fig. 1b for an increase from 7.75 to 8.67 b), because it increases the overall number of events occurring at higher impact parameters.
We follow now Fig. 2, showing experimental and model results for quantities that characterise the bulk of particle production from the measured spectra.The larger σ tot for the larger neutron skin induces larger impact parameters at the same centrality.As a consequence, fewer particles are produced for larger values of ∆r np , as clearly visible in the total multiplicity in Fig. 2 (left panel).A second effect of a larger skin, highlighted in Fig. 1b, is that it leads to more diffuse QGP droplets, which leads to weaker pressure gradients and a slower hydrodynamic expansion.This translates into a lower average momentum of the detected particles, as seen in the middle panel of Fig. 2. In addition Fig. 1 shows that a larger neutron skin reduces the ellipticity of the QGP.This leads to a reduction of the elliptic flow, measured in experiment as a two-particle azimuthal correlation (v 2 {2}, the rms value  1)) is highly anti-correlated with ∆rnp, as both have a strong effect on the centrality dependence of observables (see also Fig. 2).
of the distribution of v 2 ) or as a four-particle correlation (v 2 {4}).Indeed Fig. 2 (right) shows the expected reduction and moreover we find that a larger neutron skin enhances the difference between v 2 {2} and v 2 {4}, which corresponds to larger elliptic flow fluctuations.
Bayesian inference of the 208 Pb neutron skin -Due to the interplay and cross-correlations among parameters and observables, constraining the model from experiment requires advanced Bayesian analysis tools as pioneered in earlier works [15,21].Our analysis makes use of 653 data points in 208 Pb+ 208 Pb collisions and a single data point (the total cross section) of protonnucleus (p+ 208 Pb) collisions.We use 3000 design points for the Gaussian Processes to emulate our collisions as a function of the 26-dimensional parameter space.See the Supplemental Material (SM) for a specification of all data, parameters and their inferred distributions.
The posterior distributions are displayed in Fig. 3 for a subset of parameters that correlate highly with ∆r np .These are the parameters appearing in the energy deposition formula, Eqn.(1), namely, the energy deposition parameters p and q, as well as the nucleon size, w.In fact, the p parameter and ∆r np are the most negatively correlated across our entire parameter space.This is not surprising, as both parameters strongly influence the centrality dependence of observables, whereby a larger neutron skin in particular affects off-central collisions by increasing the total cross section.
Here we briefly revisit Fig. 2, where the middle curve represents our most likely value estimate.In the SM we present the full posterior distributions of our set of 653 data points.There, it can also be seen that the reason for the mismatch between the computed ⟨p T ⟩ and experimental data in Fig. 2 lies in a slight overestimate of yield of protons, which comes with a larger p T .There is also a significant posterior uncertainty in the anisotropic flow, which is dominated by the emulator uncertainty.
In Fig. 4 we put our new result in context of other state-of-the-art determinations of the skin of 208 Pb.From the posterior distribution we obtain ∆r np = 0.217± 0.058 fm, corresponding to a point-like rms neutron radius r n = 5.653±0.058fm.Our result is compatible with both the ab initio determination [7] and the PREX result [6], which is competitive in precision.With regards to the EOS of neutron matter, from the relation between ∆r np and the slope parameter, L, of the symmetry energy around the nuclear saturation density [22], we obtain L = 79 ± 39 MeV, representing the first collider-based constraint on this parameter from high-energy data.
We comment now on the robustness of this result.The total 208 Pb+ 208 Pb and p+ 208 Pb cross sections [23,24] pose important constraints on the neutron skin.Indeed, excluding these two measurements we obtain ∆r np = 0.31 ± 0.10 fm, whereas using exclusively these two data points results in ∆r np = 0.03±0.12fm.Our result comes hence from constraints due to a combination of observables, where the cross section prefers a smaller neutron skin, while other observables prefer a larger value (this is similar for w [25]).For the first time in Bayesian analyses we include the ρ 2 observable [26,27], a sensitive probe of the initial conditions [25,[28][29][30][31] [6] and the estimate of ab initio nuclear theory (with an error bar corresponding to a 68% credibility interval) [7].
sures the correlation between v 2 {2} and ⟨p T ⟩.Without this observable, the analysis yields a consistent result, ∆r np = 0.243±0.059fm.Also, as introduced in Ref. [25], we weight the targeted observables according to a prescription that models unknown theoretical uncertainty with respect to p T -differential observables in particular.Turning this weighting off, we find a consistent albeit slightly smaller neutron skin, ∆r np = 0.160 ± 0.057 fm.Further indication of the robustness of our finding comes from the fact that targeting a subset of p Tintegrated-only observables, corresponding to 233 AL-ICE data points, we obtain ∆r np = 0.216 ± 0.057 fm.This suggests that the extraction of ∆r np is likely insensitive to theoretical uncertainties in the particlisation of the QGP at the switching temperature [32].Lastly, we note that our T R ENTo Ansatz of Eqn. ( 1) is very versatile, and may lead to a relatively conservative estimate of the uncertainty on ∆r np .Implementing in the future a model of initial conditions motivated by first-principles arguments and with fewer parameters [33], may lead to stronger constraints than presented here.
Future skin determinations at the LHC -Albeit computationally expensive, it would be interesting to vary the full neutron profile, by changing R n , or via a multi-parameter function as in Ref. [2].The SM presents an exploratory study of this kind.Increasing the neutron half-width radius does not affect the average multiplicity and the elliptic flow, but leads to a decrease in the mean transverse momentum and a higher hadronic cross section.We speculate, thus, that a more complete analysis could eventually lead to a slightly smaller neutron skin.
We expect our analysis to trigger a program of complementary measurements of skin effects at the LHC.A method pioneered by the STAR collaboration utilises the photo-production of vector mesons in ultra-peripheral nucleus-nucleus collisions to infer the average gluon density in the colliding nuclei, and hence the neutron skins [34].The extracted skin of 197 Au is in good agreement with nuclear theory predictions [35].Therefore, the same method could be exploited at the LHC to perform an independent extraction of the skin of 208 Pb.
In addition, the global analysis presented here uses so-called soft observables that depend on particles with transverse momentum of order of the QCD deconfinement temperature, around 150 MeV.With highluminosity LHC runs it may be possible to constrain the neutron skin as well via hard observables, such as high transverse momentum electroweak bosons [36,37].The charge of the produced electroweak bosons can serve as a direct probe of the number of neutron-neutron interactions.By selecting collisions at relatively large impact parameter, it is then possible to determine the dominance of neutrons at the outer edges of the 208 Pb nucleus.
It is likely that the nucleus 48 Ca and other ions will be collided at the high-luminosity LHC in the next decade [38].This will enable an extended analysis that in particular can be compared with the dedicated CREX measurement of the neutron skin of 48 Ca [39].Comparing many different collision systems will furthermore permit us to study ratios of observables that cancel most of the systematic theoretical uncertainties [40][41][42][43], leading to improved determinations of ∆r np across the nuclear chart.
The Trajectum framework consists of an initial stage, a hydrodynamic stage and finally a freeze-out to hadrons that are then evolved by the SMASH code [44][45][46].Here we give a brief summary of all parameters involved, how they are constrained, and by which experimental data.All parameters are displayed in boldface.Full details can be found in Ref. [18].
For the initial stage the nucleons are distributed within the colliding nuclei according to Eqn. 2 with parameter a n , whereby each sampled nucleon pair has a minimum distance of d min .In addition, in Eqn. 2, we consider the half-width radius expanded in spherical harmonics up to the axial quadrupole deformation, , where β 2 quantifies the magnitude of the ellipsoidal deformation.Given the structure of 208 Pb [47], we let β 2 fluctuate around zero with a standard deviation ⟨β 2 2 ⟩ − ⟨β 2 ⟩ 2 .The nucleons that collide are determined by the measured σ NN cross section as in [17,25] and are called participants.Each participant then consists of n c constituents, each associated with a transverse Gaussian profile with width v = 0.2 fm + χ struct (w − 0.2 fm).The center coordinates of the nucleon constituents are distributed according to a Gaussian distribution.This leads to an average Gaussian nucleon profile with a width w.Superimposition of the nucleon constituent profiles then leads to the thickness functions T in Eqn. 3. The normalisation of each Gaussian profile fluctuates and is equal to N γ/n c , where γ is sampled from a gamma distribution with mean 1 and width σ fluct √ n c .The final energy density is then given by Eqn. 3 (main text) with parameters p and q.
For a time τ hyd this energy density then either evolves via free streaming or according to a holographic prescription that starts assuming viscous hydrodynamics (see [18]).The way this is chosen is with a continuous parameter r hyd between 0 and 1 which interpolates between free streaming (r hyd = 0) and the holographic prescription (r hyd = 1).
Next the stress-energy tensor is evolved according to second order hydrodynamics with 9 parameters that determine the transport coefficients, namely η/s, (η/s) slope , (η/s) δslope , (η/s) 0.8 GeV , (ζ/s) max , (ζ/s) m×w , (ζ/s) T0 , τπsT η and τππ τπ .Of these the first four determine the shear viscosity over the entropy density, η/s, as a function of the temperature, by fixing the average value, the slope between temperatures of 150 and 300 MeV, the difference in slope at higher temperatures and its value at 800 MeV and higher.The next three determine the maximum of the bulk viscosity over the entropy density, ζ/s, the width of its peak times the maximum, and the location of the maximum, respectively.The last two parameters are dimensionless combinations of first and second order transport coefficients, whereby in particular the first determines the relaxation rate towards first order viscous hydrodynamics.In the hydrodynamic phase the equation of state can vary with the parameter a EOS , which varies a n in the EOS prescription of [15] (not to be confused with the parameter a n discussed in the main text) while keeping the degrees of freedom at high temperatures fixed.
At the switching temperature T switch the fluid is frozen out to hadrons according to the Cooper-Frye procedure [48] with the Pratt-Torrieri-Bernhard (PTB) prescription for viscous corrections [15,49].Lastly, the hadrons are evolved with the SMASH code, whereby we scale all interaction cross sections by a factor f SMASH .
Finally, in the analysis the number of events can vary with a normalising parameter cent norm , which is also the only parameter that has a non-trivial prior likelihood distribution, namely a Gaussian of width 1%.
We then compare the output of the code with experimental data.In this work we use the hadronic cross section measurements for 208 Pb+ 208 Pb and p+ 208 Pb collisions [23,24], the charged particle yields at 2.76 [50] and 5.02 TeV [51], the identified particle yields dN ch /dy and mean transverse momentum ⟨p T ⟩ for pions, kaons and protons at 2.76 TeV [52], as well as unidentified transverse energy E T [53] and fluctuations of mean p T [54] at 2.76 TeV.We also include the integrated anisotropic flow coefficients v 2 {2}, v 2 {4}, v 3 {2} and v 4 {2} at both 2.76 and 5.02 TeV [55] and p T -differential observables with bin boundaries at (0.25, 0.5, 0.75, 1.0, 1.4, 1.8, 2.2, 3.0) GeV.In particular, this includes spectra for pions, kaons and protons at 2.76 TeV [52], as well as v 2 {2}(p T ) for pions, kaons and protons, and v 3 {2}(p T ) for pions (these data are only available for p T > 0.5 GeV) [56].Finally we also include "statistically difficult" observables being the ρ 2 (v 2 {2}, ⟨p T ⟩), N SC(2, 3) and N SC(2, 4) correlators, where N SC(i, j) measures the normalised correlations between v 2 i and v 2 j [27,57].The posterior for our parameters is then given by Bayes' formula with P(x) the (flat) prior probability density and where with y(x) the predicted data for parameters x, y exp the n experimental data points and Σ(x) is the sum of the experimental and theoretical covariance matrices.The covariance matrices are constructed as in [15].The standard procedure is then to run the model at a number of design points in the parameter space as determined by a latin hypercube and use those design points to construct an emulator for y(x) to evaluate Eqn.(3) using the parallel tempered emcee code [58,59] work we used 3000 design points whereby each design point has 1M initial stages of which we evolve 18k using hydrodynamics and finally simulate about 100k SMASH events to get high statistics even for ultracentral collisions.Every 1 in 15 design points uses 10 times more statistics to allow for "statistically difficult" observables that are then emulated using 200 design points.The final posterior distributions and their correlations for the standard settings of our fitting procedure are displayed in Fig. 5.We refer to the main text of the paper for details on variations in these settings.
In Fig. 6 we compare the resulting posterior likelihood distributions for the proton and neutron densities to measurements, as presented in [2].Unlike the analysis of [2], we vary only a n , so that consequently our neutron density is much more restricted and the proton density is fixed.Nevertheless, in the region where our variation is relevant (i.e., where the diffuseness of the neutron skin matters) the results of the two methods agree remarkably well both in values and in uncertainties.
It is an important question what could happen when varying the neutron profile with more parameters, such as done in [2].An indication can be found in Fig. 7, where we vary R n in the Woods-Saxon profile such that the total ∆r np variation matches Fig. 2 in the main text of the paper.It can be seen that the average multiplicity and the elliptic flow remain unaffected, while a smaller R n results in a larger mean transverse momentum.This is in agreement with the impact of this parameter on the number of participant nucleons and initial energy density eccentricity, which we verified to be also unaffected (note, however, that a smaller nucleus leads to a higher number of binary collisions).It is important that increasing R n in this figure increases as well the total hadronic cross section from 7.94 to 8.40 b, so that even if the multiplicity and the anisotropic flow remain unaffected, a potential Bayesian analysis would penalize the large R n by the measurement of 7.67 b [25,60].Based on this figure a somewhat smaller ∆r np would be possible, but perhaps with a larger bulk viscosity to compensate for the higher mean transverse momentum.A full analysis is however beyond the scope of this work.FIG. 6. Neutron density measured in polarised proton scattering from Zenihiro et al [2] compared to that extracted from LHC data.Since we use a single parameter variation (an in Eqn. 2) our profile does not vary in the central region.Nevertheless, around the skin at r ≈ 7 fm our likelihood band is consistent with [2].The accurately known proton densities (green) are also in good agreement.
the considered observables.Particularly important in relation to the main text is that the mean transverse momentum for identified pions, kaons and protons have an excellent description (row 4, column 3) as opposed to the charged hadron mean transverse momentum in the main text.The reason here is that our model overestimates the number of protons (row 5, column 1) by about 10%.Given the high proton p T , this small excess leads to the overestimate of experimental data observed in Fig. 2 of the main text.This issue is mainly due to effects in the hadronic phase, such as proton-antiproton annihilation and regeneration (see e.g.[61]), whose detailed inclusion is beyond the scope of Trajectum.FIG. 8. We show all the observables used in our analysis from measurements by the ALICE collaboration [50][51][52][53][54][55][56] and the ATLAS collaboration [27,57], with the exception of the total hadronic cross sections (see [25]).Theoretical uncertainties come from the posterior distribution from Fig. 5, which includes an emulator uncertainty that is dominant for most flow observables.

FIG. 1 .
FIG. 1. Neutron skin and collective flow in relativistic nuclear collisions.a: Two ions collide with impact parameter b = 8 fm.Both ions are Lorentz-contracted by a factor γ ≈ 2500, and the relevant dynamics hence effectively takes place in the transverse plane, x ⊥ = (x, y).b:The collision deposits energy in the interaction region depending on the extent of the neutron skin of the 208 Pb nuclei.We consider ∆rnp = 0.086 fm (top) and ∆rnp = 0.384 (bottom).The neutron skin is varied by keeping the half-width neutron radius, Rn, constant while changing the neutron diffuseness, as displayed by the dotted lines (see also Eqn. (2) below).A larger neutron skin leads to a considerably larger total hadronic cross section, σtot, and the resulting QGP is in addition more diffuse and less elliptical.c: We show a single QGP evolving hydrodynamically and being converted into particles (marked in the figure with their respective symbols) as it cools, while expanding both in z and in the transverse plane.The observation of millions such events leads to characteristic azimuthal anisotropies in the momentum distribution of the produced particles, the most important of which is quantified by the rms value of its second Fourier component, the elliptic flow v2{2}, which reflects the ellipticity of the QGP.

FIG. 3 .
FIG.3.Inferred neutron skin and energy-deposition parameters.We show the posterior distribution of the neutron skin ∆rnp, the nucleon width w and the energy deposition parameters p and q, together with their expectation values (see top) and correlations.Uncertainties correspond to the standard deviations of the posterior distributions.Especially the p parameter (see Eqn. (1)) is highly anti-correlated with ∆rnp, as both have a strong effect on the centrality dependence of observables (see also Fig.2).

FIG. 5 .
FIG. 5. Correlation matrix among all 26 model parameters.Detailed information about the prior ranges is provided in Tab.I.
which mea-FIG. 4. State-of-the-art determinations of the neutron skin of 208 Pb.We show the final likelihood distribution of the neutron skin as determined from the LHC data as compared to the values obtained experimentally by the PREX collaboration (including both experimental and theoretical uncertainties in the extraction) . In this Δr np [fm] .7.We vary Rn such that the total ∆rnp variation matches Fig.2in the main text, keeping all other parameters fixed.The size itself has a visible effect on the mean transverse momentum, but little effect on either multiplicity or anisotropic flow. FIG