Methods for systematic study of nuclear structure in high-energy collisions

TL;DR

This paper introduces an efficient Monte Carlo-based method to systematically study nuclear structure in high-energy collisions, significantly reducing computational demands and enabling detailed analysis of nuclear properties.

Contribution

It proposes a novel nucleon position shifting technique in Monte Carlo simulations to facilitate systematic nuclear structure studies with fewer simulated events.

Findings

Method reduces the number of required simulated events by orders of magnitude.

Enables detailed sensitivity analysis of nuclear properties in high-energy collisions.

Facilitates systematic exploration of nuclear structure using high-precision collision data.

Abstract

There is increasing interest in using high-energy collisions to probe the structure of nuclei, in particular with the high-precision data made possible by collisions performed with pairs of isobaric species. A systematic study requires a variation of parameters representing nuclear properties such as radius, skin thickness, angular deformation, and short-range correlations, to determine the sensitivity of the various observables on each of these properties. In this work we propose a method for efficiently carrying out such study, based on the shifting of positions of nucleons in Monte-Carlo samples. We show that by using this method, statistical demands can be dramatically reduced -- potentially reducing the required number of simulated events by orders of magnitude -- paving the way for systematic study of nuclear structure in high-energy collisions,

Tables

Table 1. Table 1: Sample parameter set relevant for studying the isobar systems of 96 Ru and 96 Zr of Ref.  STAR:2021mii . The resulting densities are illustrated in Fig.  3 .

System	$R$ (fm)	$a$ (fm)	${\beta }_2$	$\gamma$ (∘)	${\beta }_{3,0}$
96Ru	5.09	0.46	0.16	30	0
96Zr	5.02	0.52	0.06	0	0.20

Table 2. Table 2: Reduction in the required number of events (“Improv. Factor”) for computing the ratio of RMS eccentricity ε n ​ { 2 } subscript 𝜀 𝑛 2 \varepsilon_{n}\{2\} using the method of shifting nucleons compared to independent nuclei, in a simple benchmark test of b = 0 𝑏 0 b=0 collisions in a participant Glauber model. Nuclei were prepared from an uncorrelated, spherical distribution and modified so that a single parameter becomes non-zero. Also shown is the corresponding change in ε n ​ { 2 } subscript 𝜀 𝑛 2 \varepsilon_{n}\{2\} and the average nucleon shift when changing each parameter from a value of 0. For the β 3 subscript 𝛽 3 \beta_{3} change, column 3 represents the resulting change in ε 3 ​ { 2 } subscript 𝜀 3 2 \varepsilon_{3}\{2\} , otherwise it is ε 2 ​ { 2 } subscript 𝜀 2 2 \varepsilon_{2}\{2\} . The improvement factor does not change significantly with harmonic, so only column 3 is harmonic-dependent. For the case of correlated nucleons, the strength parameter is set to C str = − 1 subscript 𝐶 str 1 C_{\rm str}=-1 , corresponding to almost complete exclusion within the correlation length.

	Param.	${\epsilon }_n\left\{2\right\}$	Improv.	Avg.
Par.	Change	Change	Factor	Shift
${\left({\beta }_2\right)}^2$	${\left(0.005\right)}^2$	0.02%	170	0.008 fm
${\left({\beta }_2\right)}^2$	${\left(0.01\right)}^2$	0.10%	100	0.02 fm
${\left({\beta }_2\right)}^2$	${\left(0.02\right)}^2$	0.39%	42	0.03 fm
${\left({\beta }_2\right)}^2$	${\left(0.05\right)}^2$	2.3%	12	0.08 fm
${\left({\beta }_2\right)}^2$	${\left(0.1\right)}^2$	8.8%	4.7	0.17 fm
${\left({\beta }_2\right)}^2$	${\left(0.2\right)}^2$	31%	2.1	0.33 fm
${\left({\beta }_3\right)}^2$	${\left(0.01\right)}^2$	0.05%	79	0.01 fm
${\left({\beta }_3\right)}^2$	${\left(0.05\right)}^2$	1.6%	13	0.06 fm
${\left({\beta }_3\right)}^2$	${\left(0.1\right)}^2$	6.3%	5.0	0.12 fm
${\left({\beta }_3\right)}^2$	${\left(0.2\right)}^2$	23%	2.2	0.25 fm
${\left(C_{\mathrm{len}}\right)}^3$	(0.2 fm$){}^3$	0.13%	2900	0.002 fm
${\left(C_{\mathrm{len}}\right)}^3$	$\times 2$	0.27%	1100	0.005 fm
${\left(C_{\mathrm{len}}\right)}^3$	$\times 4$	0.53%	350	0.009 fm
${\left(C_{\mathrm{len}}\right)}^3$	(0.4 fm$){}^3$	1.1%	180	0.017 fm
${\left(C_{\mathrm{len}}\right)}^3$	$\times 2$	2.0%	98	0.032 fm
${\left(C_{\mathrm{len}}\right)}^3$	$\times 4$	3.8%	54	0.059 fm
${\left(C_{\mathrm{len}}\right)}^3$	(0.8 fm$){}^3$	7.3%	25	0.11 fm
${\left(C_{\mathrm{len}}\right)}^3$	$\times 2$	14%	13	0.19 fm

Equations

$$\rho (r)$$

$$N(R,a)$$

$$P_s(\mathbf{x})=P_s(r)$$

$$P_g(\mathbf{x})=P_g(r)$$

$${\rho }_c(\mathbf{x})$$

$$\displaystyle=\frac{3}{8\pi^{3/2}R_{s}^{3}}\biggl{[}\frac{\sqrt{2}w}{r}\left(e^{-\frac{(r+R_{s})^{2}}{2w^{2}}}-e^{-\frac{(r-R_{s})^{2}}{2w^{2}}}\right)$$

$$\displaystyle\ \ +\sqrt{\pi}\left\{{\rm erf}\left(\frac{r+R_{s}}{\sqrt{2}w}\right)-{\rm erf}\left(\frac{r-R_{s}}{\sqrt{2}w}\right)\right\}\biggr{]}$$

$$D_{\mathrm{KL}}({\rho }_1||{\rho }_2)$$

$$R_s(R,a)$$

$$w(R,a)$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0.$$

$$d\mathbf{x}=\mathbf{v}dt.$$

$$\mathbf{v}(r,\theta ,\varphi )$$

$$\frac{\partial \rho }{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\rho v_r\right),$$

$$v_r(r,t)$$

$$\rho (r,t)$$

$$r_s$$

$$R_s$$

$$\frac{\partial \rho }{\partial t}$$

$$v_r(r,t)$$

$$R\to R\left(1+\sum _{\ell ,m}{\beta }_{\ell ,m}{Y}_{\ell ,m}(\theta ,\varphi )\right),$$

$$\tilde{\rho }(r,\theta ,\varphi )$$

$$\rho (\mathbf{x},t)$$

$$\mathbf{v}=\nabla \Phi .$$

$$\frac{\partial \rho }{\partial t}+\nabla \rho \cdot \nabla \Phi +\rho \Delta \Phi =0.$$

$$\Phi (r,\theta ,\varphi )=\sum R{\beta }_{\ell ,m}{f}_{\ell ,m}(r)Y_{\ell ,m}(\theta ,\varphi )$$

$$f_{\ell ,m}^{\prime \prime }+f_{\ell ,m}^{\prime }\left(\frac{2}{r}+\frac{{\rho }^{\prime }}{\rho }\right)-\frac{\ell (\ell +1)}{r^2}{f}_{\ell ,m}-\frac{{\rho }^{\prime }}{\rho }$$

$$\underset{r\to \infty }{\lim }\frac{{\rho }^{\prime }}{\rho }$$

$$\underset{r\to 0}{\lim }\frac{{\rho }^{\prime }}{\rho }$$

$$\underset{r\to 0}{\lim }{f}_{\ell ,m}=C{r}^{\ell }+D{r}^{1-\ell }$$

$$f_{\ell ,m}^{\prime \prime }+\frac{1}{a}{f}_{\ell ,m}^{\prime }-\frac{1}{a}$$

$$\underset{r\to \infty }{\lim }{f}_{\ell ,m}=r+C^{\prime }{e}^{\frac{r}{a}}+D^{\prime }.$$

$$f_{\ell ,m}(r\to 0)$$

Full text

1 Introduction

An emerging direction of research in recent years is to use data from high-energy nucleus-nucleus collisions in order to infer properties of colliding nuclei, such as their deformation Giacalone:2019pca ; Giacalone:2020awm ; Giacalone:2021udy ; Bally:2021qys ; Jia:2021qyu ; Giacalone:2021uhj ; Zhang:2021kxj ; Bally:2023dxi ; Samanta:2023tom . In particular, recent measurements in isobaric collisions of 96Ru and 96Zr at the Relativistic Heavy-Ion Collider STAR:2021mii have shown that small differences in the nuclear properties are potentially measurable Hammelmann:2019vwd ; Xu:2021vpn ; Nijs:2021kvn ; Zhao:2022uhl ; Jia:2022qgl ; Liu:2022kvz ; Nie:2022gbg .

For such studies, it is desirable to have the ability to continuously change nuclear structure parameters, and precisely quantify the small differences they induce on collision observables. Modern modeling of high-energy collisions involves Monte Carlo simulations Moreland:2018gsh ; JETSCAPE:2020mzn ; Nijs:2021kvn ; Parkkila:2021yha . If a new, independent set of Monte Carlo simulations is generated for each parameter value, a huge amount of statistics is required to resolve the small differences in final observables. For example, in Ref. Jia:2022qrq , hundreds of millions of simulation events were generated for each set of possible nuclear parameter values, in order to study the precise observable ratios measured in the RHIC isobar run.

We introduce a method which significantly reduces the required number of simulations, potentially by orders of magnitude. The most important source of event-by-event fluctuations in a high-energy collision is the random position of each nucleon at the time of collision PHOBOS:2006dbo ; Alver:2010gr ; Luzum:2011mm ; Luzum:2013yya . We show that changes in nuclear properties can be implemented by slightly shifting these positions.

Nucleon positions are typically sampled independently according to a Woods-Saxon distribution. We first introduce, in Sec. 2, a fast method for carrying out this sampling. In Sec. 3, we describe how nucleon positions must be shifted in order to implement a change in the one-body distribution, such as a global deformation. In Sec. 4, we describe how they must be shifted in order to implement short-range two-body correlations. The application of our method to heavy-ion collisions is discussed in Sec. 5. In Sec. 6 we discuss a complementary method to explore small regions in parameter space without any extra simulations by reweighting events. In Sec. 7, we show on specific examples that our methods are considerably more efficient than traditional methods.

2 Preparing the nucleus

The idea is to prepare a set of nuclear configurations — the position of each nucleon in a nucleus — each consisting of $A$ independent nucleons governed by a spherically-symmetric probability distribution $\rho$. This existing set can then be modified by changing the positions of these nucleons, such that the new set of nuclei respect a different 1-body density $\rho$, as well as a non-trivial correlation function $C(r)$.

The typical profile used to describe a nucleus is a Woods-Saxon

[TABLE]

which is parameterized by a radius $R$ and diffusiveness parameter $a$. The proportionality factor $N(R,a)$ is defined by noting that total probability must sum to unity: $\int d^{3}x\,\rho=1$. Explicitly,

[TABLE]

where Li3 is the polylogarithm of order 3.

This distribution can be independently sampled (e.g., using acceptance-rejection sampling ) in order to generate a discrete set of nuclei governed by distribution $\rho$.

Here we present an alternative method for sampling (and manipulating) a spherically-symmetric distribution. This is done by approximating a Woods-Saxon distribution by the convolution of a spherical step function and a 3D Gaussian. This has the advantage of giving closed-form analytic expressions, which gives a high degree of analytic control. For example, the relation between the point nucleon density and, e.g., the charge density can be written analytically. Additionally, nuclear configurations can be computed without rejection sampling of a Woods-Saxon distribution. Finally, parameters $R,a$ can be easily changed.

The idea is to treat the coordinate of each nucleon as a random variable that, rather than being governed by a Woods-Saxon, is a sum of two other random vectors governed by a spherical step function and spherical Gaussian, respectively.

[TABLE]

with $r=\lvert{\bf x}\rvert$. The position of a nucleon is given by the vector sum of positions drawn from each of these distributions. The sum of the two random vectors thus follows the convolution

[TABLE]

With the correct choice of radial and width parameters $(R_{s},w)$ a Woods-Saxon can be well approximated by this function $\rho_{c}(r)\simeq\rho(r)$.

A natural way to compare probability distributions is the Kullback–Leibler (KL) divergence Kullback:1951zyt

[TABLE]

Two identical distributions have KL divergence $D_{\rm KL}(\rho_{1}||\rho_{1})=0$, while increasingly different distributions cause it to increase.

We can choose the best parameters $(R_{s},w)$ to describe a given Woods-Saxon $(R,a)$ by minimizing the KL divergence $D_{\rm KL}(\rho,\rho_{c})$. A rough rule of thumb for the conversion when $a/R$ is small is

[TABLE]

In practice we obtain the correct values on-the-fly with numerical minimization.

A comparison to traditional Woods-Saxon distributions is shown in Fig. 1. The fit is worst for certain large values of skin thickness $a$ (0.2–0.3 fm), but works well for values of interest.

In all of the numerical examples shown in this work, we use this method for generating seed nucleon configurations, which can then be modified as desired by shifting nucleon positions.

3 Changing shape

Once a discrete set of nuclear configurations are prepared, they can be modified to follow the desired statistics.

First we consider a continuous change in the 1-body density — that is, a change that can be characterized by some continuous parameter $t$, $\rho=\rho({\bf x},t)$.

We will be transporting particles in order to affect the change in density that is brought about by a change in parameter $t$. Conservation of particles demands that this transformation obeys a continuity equation

[TABLE]

where the shift in position $d{\bf x}({\bf x},t)$ of particles with respect to parameter $t$ is described by vector field ${\bf v}({\bf x},t)$,

[TABLE]

3.1 Radial deformation

We start by describing how to modify the parameters $R,a$ of a spherical Woods-Saxon, in the case where it is desired to start with a given Woods-Saxon sampling instead of the alternate step+Gauss distribution (6).

In this case, there is rotational symmetry both before and after the transformation, and we can choose a purely radial shift for the nucleons.

[TABLE]

The continuity equation (10) becomes

[TABLE]

which we can simply integrate to obtain

[TABLE]

For the specific case of a Woods-Saxon distribution, we note that for fixed value of $a/R$, different distributions are related to each other by a simple scale transformation $r\to Cr$.

To modify $a/R$, we define

[TABLE]

so that

[TABLE]

3.2 Angular deformation

Starting with a set of spherical nuclei prepared using one of the above methods, we can modify the positions of nucleons to obtain nuclei described by a deformed distribution. A nucleus with angular deformation is typically parameterized by deforming a Woods-Saxon (1), replacing the radius parameter

[TABLE]

where $Y_{\ell,m}(\theta,\phi)$ are the real form of spherical harmonics, and each multipole component is characterized by a coefficient $\beta_{\ell,m}$. That is,

[TABLE]

We first note that we can relate the deformed distributions to the spherical Woods-Saxon (1) with a coordinate transformation

[TABLE]

We also note that the normalization of the probability distribution depends (weakly) on coefficients $\beta_{\ell,m}$ in addition to parameters $R$, $a$. However, our transformation will consist only of transporting particles, so probability is naturally conserved. We therefore omit the normalization factor, which will cancel in the final expressions.

We therefore define our $t$-parameterized distribution as

[TABLE]

so that $t=0$ represents a spherical distribution and $t=1$ represents the desired deformed distribution, with the two limits continuously connected by intermediate values.

We then determine a vector field ${\bf v}$ which satisfies the continuity equation (10) for the above expression of $\rho({\bf x},t)$. This single first-order differential equation for the vector field ${\bf v}$ is an underdetermined problem, with no unique solution. We posit a solution of the form of the gradient of a scalar field $\Phi$.111In the language of fluid dynamics, this corresponds to irrotational flow.:

[TABLE]

An advantage of this prescription is that symmetry arguments can be used to specify the dependence of $\Phi$ on the coordinates $(r,\theta,\phi)$. In the case of a radial deformation, for instance, the velocity field Eq. (12) can be written as a gradient, where $\Phi$ solely depends on the radial coordinate $r$.

Inserting Eq. (24) into the continuity equation (10), one obtains the following equation for $\Phi$:

[TABLE]

If the deformation is sufficiently small, we can linearize this equation and replace $\rho$ with the original spherical distribution (corresponding to the value at $t=0$) in the last two terms. The second term then reduces to $\rho^{\prime}(r)\partial\Phi/\partial r$.

In the case of an angular deformation, we decompose $\Phi$ into multipole components

[TABLE]

Inserting this equation into Eq. (25) and using Eq. (23), the equations for the different multipole component decouple, and they satisfy the following second-order equation in $r$

[TABLE]

Note that $f_{\ell,m}$ does not, in fact, depend on $m$.

Defining a unique solution requires two boundary conditions. For this, we consider the behavior at small and large $r$. Note first that

[TABLE]

since we expect $R\gg a$.

So in the limit of small $r$, $f_{\ell,m}$ is approximately a solution of Laplace’s equation

[TABLE]

At large $r$, we have

[TABLE]

with general solution

[TABLE]

We would like to choose boundary conditions which eliminate both the divergence at small $r$, that is, $D=0$ and the divergence at large $r$, that is, $C^{\prime}=0$. This can be achieved with conditions

[TABLE]

The general solution can not be found analytically, and must instead be solved numerically. In practice, this can be done by separately solving the full inhomogeneous and the homogeneous equations, integrating backward from some large $r_{\rm max}$, and taking the linear combination of the two solutions that satisfies the desired boundary condition at some small $r_{\rm min}$. That is, we solve

[TABLE]

and take as the final solution the linear combination

[TABLE]

The solution is very stable and typically does not have any significant dependence on $r_{\rm min}$ or $r_{\rm max}$.

The final shift is thus

[TABLE]

For small deformation, we can set $t=1$, shift all particles, and the new set of configurations will consist of independent nucleons with (approximately) the desired 1-body probability density. For illustration, in Fig. 2 we show example vector plots of the shift field $d{\bf x}$ for the case of $\beta_{2,0}$ only and $\beta_{3,0}$ only.

For larger deformations ($\sum\beta\gtrsim 0.2$), it is advantageous to break into multiple ($N$) steps of size $t/N$, evaluating the shift $d{\rm x}$ at the new position after each step. In principle this requires solving the differential equation for ${\bf v}$ at the new nonzero value of $t$ for each subsequent step. In this case, the equations for different multipoles $f_{\ell,m}$ become coupled to each other and it is necessary to simultaneously solve the coupled differential equations. In practice, good accuracy can be obtained by breaking the shift into several steps, but reusing the same ($t=0$) shift field to compute each step.

For illustration, in Fig. 3 we explicitly show the resulting density from the Monte Carlo implementation in the relevant case of the isobar pair 96Ru and 96Zr from Ref. STAR:2021mii using parameters shown in Table 1. Note that quadrupole deformations are typically parameterized with magnitude $\beta_{2}$ and angle $\gamma$, related to the spherical harmonic coefficients as

[TABLE]

4 Short-range correlations

Interactions between nucleons can induce correlations. Such correlations are encoded in $N$-body distributions beyond the 1-body density $\rho{(\bf x)}$. Let $\rho_{2}({\bf x}_{1},{\bf x}_{2})$ denote the distribution of nucleon pairs, from which we can define the correlation function $C({\bf x}_{1},{\bf x}_{2})$,

[TABLE]

In particular, short-range interactions between nucleons (as well as their Fermi-Dirac statistics) are expected to induce short-range correlations between nuclei. That is, you will find fewer (or more) pairs of nucleons within a short distance $r$ from each other than would be expected from independent particles governed by distribution $\rho$.

As was possible for the 1-body distribution, we can induce a correlation between particles by slightly shifting their positions. The simplest way is to start with a collection of uncorrelated nucleons (prepared, e.g., as described in the previous sections), and change the distance between pairs of nucleons such that the resulting distribution obeys the desired correlation function.

We choose a prescription such that pairs of uncorrelated particles with relative position ${\bf r=}r\hat{r}={\bf x}_{2}-{\bf x}_{1}$ are shifted radially from their center point so that their new separation is $\tilde{r}(r)$. The new separation is chosen as a monotonically-increasing function that preserves the number of pairs,

[TABLE]

For a chosen correlation function $C$ (which in principle could depend on absolute position, orientation, spin, etc.), one can invert this expression to solve for $\tilde{r}$.

We choose a symmetric shift, so that each particle in the pair moves half the required separation shift. The total shift of a given particle $i$ is the vector sum of the shifts implied by its pairing with all the other particles in the nucleus,

[TABLE]

Further, we note that the correlation cannot be an arbitrary function. In particular, the number of nucleons (and therefore the number of pairs) in each nucleus is fixed, and this places a restriction on the correlation function. In our formalism, this corresponds to the fact that $\rho_{2}$ and $\rho$ are probabilities, which therefore integrate to unity, and which implies the sum rule

[TABLE]

In order to implement an arbitrary short-range correlation $C_{\rm short}({\bf r})$ , we add a small, constant offset $C_{\infty}$ to ensure that this sum rule is satisfied (if it’s not already), that is:

[TABLE]

Inserting this equation into Eq. (48) and assuming that the range of the correlation is much smaller than the system size, one obtains:

[TABLE]

where $C_{\rm vol}$ is the volume integral of the short-range correlation:

[TABLE]

To illustrate the method, we choose a simple example of a step-function correlation function with a variable correlation length $C_{\rm len}$ and strength $C_{\rm str}$ (see Fig. 5): That is, we define:

[TABLE]

Eq. (51) then gives

[TABLE]

We note that for the short-range correlations of interest, the offset $C_{\infty}$ defined by Eq. (50) is quite small. For example, for $C_{\rm str}=-1$ and $C_{\rm len}=0.4$ fm, we obtain $C_{\infty}\simeq 3\times 10^{-3}$, which is seemingly negligible. Nevertheless, its inclusion ensures that the 1-body distribution $\rho(\vec{x})$ is more precisely maintained.

With this choice, the final pair separation ${\bf\tilde{r}}$ for a given initial separation ${\bf r}$ is

[TABLE]

While this method does not guarantee that the 1-body distribution will remain fixed, as long as one inputs a valid correlation function that respects the sum rule (48), both the 1-body and 2-body distributions are reproduced to a good approximation for correlations of short range. See, e.g., Figs 6 and 7.

Indeed, there are advantages to this method of implementing correlations compared to other common methods, even beyond the increase of computation efficiency that is the main concern of this work. The excellent simultaneous control of the average density and correlations is one example. A common simplified method for mimicking a short-range nucleon-nucleon correlation is to simply disallow nuclei with nucleon pairs less than some exclusion distance $d_{\rm min}$. Besides being restricted to a very specific correlation function, this method can bring other problems. Implementing this in the most straightforward way Alver:2008aq , for example, results in an unwanted modification of the 1-body distribution if $d_{\rm min}$ is not sufficiently small.

One can use clever methods to keep the nucleon density fixed while ensuring that no pairs have distance less than $d_{\min}$, as done in Ref. Moreland:2014oya . However this results in an uncontrolled 2-body distribution, which can develop complicated and unintended structures. Further, the method cannot be used for nuclear shapes without an axial symmetry (for example, for $\beta_{22}\neq 0$), since it involves randomly reassigning the azimuthal angle of offending nucleons.

The sampling method of Refs. Alvioli:2009ab ; Alvioli:2010yk , on the other hand, can in principle give more precise control over the 1-body and 2-body distributions. However, it remains difficult to efficiently and systematically study changes in the nuclear correlations, where the method described here has the clear advantage.

5 Application to high-energy nuclear collisions

Once the nuclear configurations are prepared, they can be used to simulate collision events. The details depend on the specific model used for the collision dynamics, but generally, in order to minimize statistical uncertainty, it is useful to reuse as many of the fluctuations as possible. Besides the nucleon-position fluctuation discussed in the previous sections, there are others that can arise.

Fluctuations that will be present in any simulation include the random impact parameter of each collision and the orientation of each nucleus that collides. One should ensure that when preparing collisions corresponding to different nuclear parameters, one uses the same pair of nuclei prepared as described above, with the same impact parameter, and the same Euler angles for each respective nucleus.

Note that this will not always be possible. For example, at a given impact parameter, one pair of nuclei may suffer a collision while a corresponding pair of deformed nuclei might not. To avoid bias, one must nevertheless include such unpaired events, despite the increase in statistical uncertainty. Similarly for models with a notion of participating versus spectator nucleons, there might be a different number of participants in corresponding collisions.

Despite this, the decrease in statistical uncertainty from using our methods can still be dramatic, as we show in Sec. 7.

The treatment of any other fluctuations will be specific to each model. Generally, one should reuse fluctuations as much as possible, using the same random number associated with each pair of nuclei (in the case of probabilistic cross section), individual nucleons, transverse position, etc., even if the relevant cumulative probability distribution is slightly different.

6 Probability Reweighting

Instead of performing separate simulations for each nucleus (i.e., each point in parameter space), it can be possible to study multiple parameter values from a single set of simulations.

The idea is the following. If we sample nuclei randomly from some $A$-body probability distribution $\rho_{A}({\bf x}_{1},{\bf x}_{2},\ldots{\bf x}_{A})$ we can reinterpret the resulting set of nucleon configurations as a weighted sampling from a different distribution $\rho^{\prime}_{A}$. The relative weight is simply $\rho^{\prime}/\rho$ evaluated at the sampled nucleon positions.

In the limit of infinite sampling, $\rho^{\prime}$ can be reconstructed exactly, as long as the distributions both have non-zero support in the same region. Thus, in principle one can explore many points in parameter space “for free”, by computing a set of simulations at one point in parameter space (i.e., from a single $A$-body density), and computing results for other parameter sets by simply reweighting the events when computing observables.

However, convergence can be very poor if the probability distributions are not similar — especially if $\rho$ has negligible density in some region where $\tilde{\rho}$ has non-negligible density. It can take many, many samples to fill out the distribution, and convergence can be extremely poor, wiping out any potential efficiency gains.

Here we make a preliminary study of this method to gauge its feasibility in the context of the nuclear structure studies of interest for this work. We use the $A$-body density

[TABLE]

with $\tilde{\rho}$ the (potentially) deformed Woods-Saxon (22) and $C$ a desired 2-particle correlation. We sample one distribution (an uncorrellated, spherical Woods-Saxon), and estimate observables relevant to a different distribution by reweighting events constructed from the sampled nuclei.

Details and numerical benchmarks are presented in Section 7.

7 Quantifying the efficiency gain

In order to quantify the efficacy of these methods, we implement a simple Glauber model and calculate the ratio of observables obtained from collisions of nuclei with two different nuclear properties.

For simplicity, we set impact parameter to 0, so that each nucleus (before sampling) is centered at the same transverse coordinate. After generating a set of nuclear configurations, we select two nuclei to collide and give each nucleus a random (3D) orientation.

A nucleon-nucleon collision is said to occur if the transverse distance is $r<2\sqrt{\sigma_{\rm NN}/\pi}$, where $\sigma_{NN}$ = 0.62 fm2 represents the inelastic nucleon-nucleon cross section.

Among the set of participating nucleons, we then calculate the participant eccentricities

[TABLE]

for $n=2$ and $n=3$, where $(r_{i},\phi_{i})$ are the polar coordinates of participant nucleon $i$ in the transverse plane.

With the eccentricity computed in each collision event, we compute the RMS eccentricity

[TABLE]

The ratio of $\varepsilon_{n}\{2\}$ in two different collision systems is a good approximation for the ratio of final measured observables $v_{n}\{2\}$, and so represents a quantity of high relevance.

To estimate the uncertainty in the eccentricity ratio, we use jackknife resampling.

We quantify the efficacy of our method by comparing the statistical uncertainty for the case of traditionally-prepared nuclei — that is, where all nuclei are sampled completely independently — to the case where we use the method of shifted nucleons to relate different nuclear structure. In the latter case, we use the same seed nuclei (before shift) for both nuclei in the colliding pair, and the same rotation angle for each corresponding nucleus.

For these tests we start with a baseline spherical nucleus of uncorrelated nucleons, with Woods-Saxon parameters $R=5.09$ fm and $a=0.46$ fm. We then compare this baseline case to the case where one parameter ($\beta_{2}$, $\beta_{3}$, $C_{\rm len}$) is non-zero. For correlations we study a (near) full exclusion, $C_{\rm str}=-1$, but note that preiliminary tests indicate that eccentricities only depend on the combination $C_{\rm str}C_{\rm len}^{3}$.

The deformed/correlated nuclei are either generated completely independently or by shifting nucleons. Observables are then computed for these cases, as well as the case where the baseline spherical, uncorrelated collisions are reweighted to compute eccentricity ratios corresponding to the respective deformed/correlated cases.

In Fig 8, we show an example of adding short range correlations ($C_{\rm length}$ = 0.2 fm, $C_{\rm strength}$ = -1, numerator of eccentricity ratio) to the case where nucleons are uncorrelated. The shift in eccentricity with the addition of correlations is $\sim$0.1%, which requires a very large number of collision events to resolve using the traditional method. In contrast, using the method of shifting nucleons, the statistical uncertainty is reduced by a factor $\sim$32. The required number of events to obtain a certain statistical precision is thus reduced by a factor $\sim$2900. That is, a calculation with only $\sim$35 events has smaller uncertainty than 100000 events sampled independently. Such a dramatic decrease in computation requirements has clear implications on the feasibility of systematic study.

However, this efficiency gain depends on the size of the transformation between the two nuclei. In particular, as the typical shift in nucleon position increases in distance, there is a higher likelihood of participant nucleons turning to spectators, or vice versa. This diminishes the correlation in statistical error that exists before and after the transformation.

In order to give an idea of the efficiency gain in various contexts, we list a number of examples in Table 2.

The change in $\varepsilon_{n}\{2\}$ generally tends to scale linearly with $\beta^{2}$ and $C_{\rm str}C_{\rm len}^{3}$. The improvement factor tends to scale like 1/$\sqrt{\Delta(\beta^{2})}$ and $1/\Delta(C_{\rm str}C_{\rm len}^{3})$. That is, while smaller changes create greater difficulty for traditional independent sampling, when using the better method uncertainties are actually decreased.

For a given average nucleon shift, the study of correlations sees more efficiency gains than the study of angular deformation. This is because in the case of angular deformations, nucleons in the bulk of the nucleus do not shift, while nucleons near the edge shift more than the average (see, e.g., Fig. 2). It is precisely these nucleons at the edges that have the largest effect on eccentricities (especially in the central collisions tested here), and so the loss of syncronization from participant-spectator conversion due to the nucleon shift causes a decrease in correlation of statistical uncertainty, and therefore a decrease in efficiency gain.

The method of reweighting improves even more quickly with decreasing changes in nuclear properties than the shifting-nucleons method. So it is more efficient than shifting nucleons for small changes ($\Delta\beta\lesssim$ 0.01), but loses efficacy quickly for larger changes, becoming worse than independent sampling for $\Delta\beta\gtrsim 0.08$ or $\Delta(C_{\rm str}C_{\rm len}^{3})\gtrsim$ (0.3 fm$)^{3}$ and rapidly degrading beyond that. Thus, even though the same set of events can be used in principle to study many points in parameter space, the range of parameter space that can be efficiently explored is limited. So it can not be used as a general replacement for the shifting-nucleons method unless one wants only to explore a quite small parameter space. However, it may be very useful in conjunction with other methods. For example, larger jumps in parameter space obtained from shifting nucleons can be filled out with reweighting around each point in parameter space that is sampled and simulated. In that way, one can obtain a more precise interpolation between sampled parameter values, without computing any extra simulations.

8 Conclusions

We have introduced methods to dramatically decrease statistical demands when studying how heavy-ion collision observables depend on the properties of the colliding nucleons. These properties include the average nucleon density (parameterized by a Woods-Saxon radius $R$ and diffusiveness $a$ along with any number of angular deformation coefficients $\beta_{\ell,m}$) as well as an arbitrary short-range correlation function $C({\bf x}_{1},{\bf x}_{2})$.

In general, the efficiency gain depends on the specific context, as well as the physical model used for simulations, with the largest benefit corresponding to the study of small changes in nuclear structure. Nevertheless, these methods always reduces the statistical requirements, and quite dramatic improvements are possible, with some cases seeing a reduction of the necessary number of simulations by multiple orders of magnitude.

Armed with these techniques, a large number of detailed and systematic studies of nuclear structure in the context of relativistic heavy-ion collisions will be possible.

Acknowledgments

We thank Govert Nijs and Dean Lee for suggesting to try a reweighting method. We thank the Institute for Nuclear Theory at the University of Washington for its hospitality and the Department of Energy for partial support during the completion of this work. ML acknowledges the support of the São Paulo Research Foundation (FAPESP) under grants 2021/08465-9, 2018/24720-6, and 2017/05685-2, as well as the support of the Brazilian National Council for Scientific and Technological Development (CNPq). We thank support from the “Emilie du Châtelet” visitor programme and from the GLUODYNAMICS project funded by the “P2IO LabEx (ANR-10-LABX-0038)” in the framework “Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR), France. MH was supported by the National Science Foundation (NSF) within the framework of the MUSES collaboration, under grant number OAC-2103680.