# Revealing initial state properties through ultra-central symmetric heavy-ion
collisions

S. M. A. Tabatabaee Mehr  School of Particles and Accelerators, Institute for
Research in Fundamental Sciences (IPM), P.O. Box 19395-5531, Tehran, Iran  S.
F. Taghavi  TUM School of Natural Sciences, Technische Universität München,
Garching, Germany

###### Abstract

Heavy-ion experiments provide a new opportunity to gain a deeper understanding
of the structure of nuclei. To achieve this, it is crucial to identify
observables under circumstances that are minimally affected by the process
that leads to the initial state of heavy-ion collisions from nuclear
wavefunction. In this study, we demonstrate that when assuming scale-
invariance, the effect of this stage on the initial energy or entropy density
moments in ultra-central symmetric collisions is negligible for nucleon sizes
of approximately 0.7 fm or larger for large nuclei. By borrowing cluster
expansion method from statistical physics and using scale-invariance
assumption, we calculate the average ellipticity of initial density at the
presence of short-range correlation. We compare our calculations to Monte
Carlo studies and assess the accuracy of various methods of short-range
correlation sampling. Additionally, we find that the isobar ratio can
constrain the initial state parameters, in addition to deformation. Our study
indicates that the isobar ratios in ultra-central collisions are especially
sensitive to the fluctuation in the weight of the nuclei constituents and the
two-body correlation among nucleons. This insight is crucial for drawing
conclusions about nuclear deformations based on isobar ratios.

###### Contents

  1. I Introduction
  2. II Ultra-central symmetric collisions and the impact of scale-invariance
  3. III Many event average initial state in ultra-central symmetric collisions
    1. III.1 Two-body correlations from a statistical mechanics analogy
    2. III.2 Note on the one-body density and two-body correlation and Monte Carlo event sampling
    3. III.3 Radius and skin thickness redefinition
    4. III.4 Constituent weight fluctuation
    5. III.5 From source fluctuation of a single nucleus to initial density fluctuation
  4. IV Ellipticity fluctuation and isobar ratio
  5. V Constraining the nuclear structure via isobar ratios
  6. VI Summary
  7. A Convergence of triangularity, symmetric cumulants and symmetry plane correlation
  8. B Effect of nucleon size and impact parameter on T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ reduced thickness function
  9. C Event-by-event fluctuation and cluster-expansion
    1. C.1 Decomposition of constituents contribution and source position
    2. C.2 Source fluctuation distribution and moments average
    3. C.3 Initial state fluctuation from cluster expansion method
  10. D Analytical moments from cluster expansion
  11. E Two-body correlation in a 1D toy model
  12. F Effect of substructure, finite cross-section, and deformation
    1. F.1 Estimation of the nucleon substructure contribution
    2. F.2 Effect of finite cross-section
    3. F.3 Effect of nucleus deformation
    4. F.4 Validation via T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model

##  I Introduction

The study of nuclear structure has become an intense area of research in
recent years due to its close relation to our understanding of quantum
chromodynamics (QCD) in a non-perturbative regime [1], as well as the link
between the neutron skin of heavy nuclei and the equation of state in neutron-
rich matter, similar to that of neutron stars [2, 3]. The study of heavy-ion
collisions has also gained popularity as a method of extracting information
about the properties of the colliding nuclei [4, 5]. To accomplish this, there
is a need for a thorough comprehension of collectivity during the process.
Alternatively, identifying observables that are minimally reliant on
collective expansion can offer valuable insights into the initial state. For
instance, the ratio between collision observables of isobar elements, such as
96Ru–96Ru and 96Zr–96Zr, can provide useful information as these observables
are particularly sensitive to the deformation of these nuclei in central
collisions [6, 7, 8, 9].

It is essential to incorporate the specific characteristics of each nucleus,
such as its overall geometry, and the correlations between nucleons, into the
initial state of ion collisions. This is done through heavy-ion initial state
models like MC-Glauber [10, 11], MC-KLN [12, 13],
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ [14, 15],
and IP-Glasma [16, 17]. The initial entropy of a system is a result of the
interaction between constituents of colliding nuclei. However, the process of
entropy production is quite complex and not yet fully understood, leading to
substantial differences between initial state models.
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
takes a phenomenological approach by introducing reduced thickness function as
a function of participant’s initial thickness functions of the target and
projectile,

$T_{R}(T_{A},T_{B})=\mathcal{N}\left(\frac{T_{A}^{p}+T_{B}^{p}}{2}\right)^{1/p}.$
(1)

The reduced thickness function is proportional to the initial entropy density.
The participant thickness functions is obtained from participant nuclear
matter density as $T_{A(B)}(x,y)=\int dz\rho^{\text{part}}(\vec{x})$ where $z$
direction is along the beam axis. In the most central collisions and in the
case of large nucleon inelastic cross-section, $\rho^{\text{part}}(\vec{x})$
is equivalent to a one-body density of nuclei. In
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model and
in the present work, the one-body density is parameterized by Woods-Saxon (WS)
distribution,

$\rho_{\text{WS}}(|\vec{x}|,\theta,\phi)=\frac{\rho_{0}}{1+\exp\left[(|\vec{x}|-R(\theta,\phi))/a_{0}\right]},$
(2)

where
$R(\theta,\phi)=R_{0}(1+\beta_{2}Y_{2}^{0}(\theta,\phi)+\beta_{3}Y_{3}^{0}(\theta,\phi)+\cdots)$.
Here, $\rho_{0}$ is the nucleon density constant, $R_{0}$ is the radius,
$a_{0}$ is the skin thickness, $\beta_{2}$ is quadruple, and $\beta_{3}$ is
the octupole deformations.

The definition of Eq. (1) has scale-invariance characteristic, meaning that
when we scale $T_{A,B}\to c\,T_{A,B}$, the reduced thickness function scales
the same way, $T_{R}\to c\,T_{R}$. This property is supported by elliptic flow
measurements in ultra-central uranium-uranium collisions at RHIC [18, 19]. In
our study, we demonstrate that in ultra-central symmetric collisions (UCSC),
the actual form of $T_{R}$ has a negligible impact on the geometrical
properties of the produced entropy and its fluctuation as long as the nucleons
are large enough. Even more complicated functions for $T_{R}$ lead to the same
conclusion. When the nucleons are large enough, the details of their positions
become less important, and in UCSC, we can assume that $T_{A}(x,y)\approx
T_{B}(x,y)$. In the specific example of the scale-invariant $T_{R}$ in Eq.
(1), this leads to independence on the parameter $p$. More generally, by only
assuming the scale-invariance $T_{R}(cT_{A},cT_{B})=cT_{R}(T_{A},T_{B})$, we
can express the reduced thickness function in UCSC as
$T_{R}(T_{A},T_{A})\equiv g(T_{A})$. Here, the scale-invariance is given by
$g(cT_{A})=cg(T_{A})$, which implies that $g(x)$ is a one-dimensional, degree
one, homogeneous function. From Euler’s homogeneous function theorem, we
conclude that a linear function is the only continuously differentiable
function with this property, i.e., $g(T_{A})\propto T_{A}$. The constant of
proportionality is absorbed into the overall normalization. Therefore, the
actual functionality of the reduced thickness function is irrelevant in this
case. For large nuclei, such as Pb, Au, Ru, and Zr, the nucleon size $w\sim
0.7\,$[fm] is large enough to employ this approximation. Here, $w$ is the
width of a Gaussian profile for nucleons.

In order to access nuclear structure information, assuming scale-invariance
can help to reduce theoretical uncertainty related to entropy production in
the evolution of ultra-central heavy-ion collisions. With this assumption, we
can connect the initial density geometry and its event-by-event fluctuation to
the parameters of the nuclei density, including WS parameters and the two-body
correlations. In UCSC, the average elliptic shape of the initial entropy,
represented by $\epsilon_{2}^{2}\\{2\\}=\langle\epsilon_{2}^{2}\rangle$,
relies solely on the fluctuation for two spherical nuclei. The fluctuation
caused by point-like sources results in $\epsilon_{2}^{2}\\{2\\}=5/(7A)$
(refer to section IV). The impact of the repulsive central core of nuclei with
a radius of $d_{\text{min}}/2$ can be described as follows: The volume of
central core need to be subtracted from the available space for fluctuation.
Therefore, for small values of $d_{\text{min}}/R_{0}$ one expects decreasing
trend for $\epsilon_{2}^{2}\\{2\\}$ proportional to
$1-A(d_{\text{min}}/R_{0})^{3}$. For larger $A$ and $d_{\text{min}}$, we have
smaller space for fluctuation. By an analogy between initial state fluctuation
and ensembles of a thermal system, we use the cluster expansion method
developed in statistical mechanics to make the above estimation more rigorous
in section III.1.

Another correction arises from the contribution of nucleon interactions in
entropy production. It is assumed in Ref. [20] that the weight of each
nucleon’s contribution in the initial entropy density fluctuates based on a
Gamma function with a unit mean value and width $\sigma_{\text{fluc}}$ to
explain the large multiplicity fluctuation in proton-lead collisions. Since
$\epsilon_{2}^{2}\\{2\\}$ is a second-order moment, one would expect that the
source fluctuation $\epsilon_{2}^{2}\\{2\\}$ is proportional to the second
moment of the Gamma function, which is $1+\sigma_{\text{fluc}}^{2}$ (see
section III.4). For a fixed value of $\epsilon_{2}^{2}\\{2\\}$, the two
competing contributions lead to the equation $\sigma_{\text{fluc}}^{2}\sim
A(d_{\text{min}}/R_{0})^{3}+\text{constant}$. A similar correlation between
$d_{\text{min}}$ and $\sigma_{\text{fluc}}$ has been observed in Bayesian
analysis studies [21, 22].

After investigating the impact of short-range correlation and constituent
weight fluctuation on $\epsilon_{2}\\{2\\}$ in UCSC, we focus to the isobar
ratio. This ratio, which compares $\epsilon_{2}\\{2\\}$ between two isobar
nuclei, such as Ru and Zr, helps to minimize the impact of the collective
expansion stage. This is because of the linear hydrodynamic response
approximation, where the elliptic flow $v_{2}$, the second Fourier coefficient
of the final particle distribution in the azimuthal direction [23, 24], is
proportional to ellipticity $\epsilon_{2}$ [25]. This linear approximation is
more accurate for large nuclei and central collisions. The ratio
$v_{2}^{\text{Ru}}\\{2\\}/v_{2}^{\text{Zr}}\\{2\\}$ has been measured by STAR
collaboration [26]. Utilizing the isobar ratio, we can minimize inaccuracies
caused by assuming a finite value for nucleon width and analyze the
contribution of initial state parameters to measurements. In particular, we
find that isobar ratio is sensitive to $\sigma_{\text{fluc}}$ and
$d_{\text{min}}$, although the effect on $\sigma_{\text{fluc}}$ is more
pronounced.

The paper is structured as follows: first, in Section. II, we examine the
moments of the initial entropy density for models with and without scale-
invariance property for central and non-central and symmetric and asymmetric
collisions. We discuss an analytical approach to study event-by-event
fluctuation in section III. In Section IV, we discuss the eccentricity average
and its isobar ratio. Using the isobar ratio, we study the constraints on the
nuclear structure in Section V. Finally, we summarize our study in Section VI.

##  II Ultra-central symmetric collisions and the impact of scale-invariance

In this section, we provide evidence to support the claim that the form of the
reduced thickness function has a small impact on the geometry of the initial
entropy density in UCSC, by assuming the scale-invariance.

The scale-invariance assumption in
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
relies on elliptic flow observations in ultra-central U-U collisions [14].
Uranium nuclei have a large quadruple deformation, which means that during a
collision, the two nuclei may collide on their tips, bodies, or somewhere in
between. It can be argued that in tip-tip collisions, there are more binary
collisions compared to body-body collisions. This leads to an increase in
entropy and, consequently, more particles. This concept is implemented in MC-
Glauber, where the reduced thickness function is defined as

$T_{R}(T_{A},T_{B})\propto\frac{1}{2}\left(T_{A}+T_{B}\right)+\alpha
T_{A}T_{B},$ (3)

where $\alpha$ controls the contribution of binary collisions. In the presence
of binary collisions, events with the highest multiplicities should be
populated by tip-tip collisions with smaller ellipticity, which goes against
the STAR data from the RHIC experiment [18, 19]. It is concluded in Ref. [14]
that U-U collisions do not favor the entropy production with binary
collisions. On the other hand, by using the combination
$(N_{A}+N_{B})/2+\alpha N_{col}$, it is possible to link the collision
geometry to the final particle production. This approach has been studied in
Refs. [27, 28]. The conclusion drawn is that the centrality dependence of the
charged-particle multiplicity in Pb-Pb collisions can only be described when
$\alpha$ is nonvanishing. This implies that one need to replace
$(N_{A}+N_{B})/2$ with an unspecified function $f(N_{A},N_{B})$ to open room
for a function compatible with both U–U collision measurements and centrality
dependence of the charged-particle multiplicity in Pb–Pb.

Since tip-tip and body-body U–U collisions lead to same multiplicity, it is
possible to assume that all the collision configurations are equivalent to the
collision of $N_{A}$ one-to-one nucleon. The latter configuration suggests
that the number of produced particles should be directly proportional to the
number of nucleons, i.e., $N_{\text{ch}}\propto N_{A}$, where the constant of
proportionality depends on the one-to-one nucleon interactions. Based on this
argument, one can suggest that the scaling property holds true for the
collision of any configuration of $N_{A}$ and $N_{B}$ nucleons:

$N_{\text{ch}}\propto f(N_{A},N_{B})$ (4)

where $f$ has a scale-invariant property, $f(cN_{A},cN_{B})=cf(N_{A},N_{B})$,
leading to linear scaling $N_{\text{ch}}\propto N_{A}$ when $N_{A}=N_{B}$. It
is important to note that the reverse logic is not valid, and we can have
functions $f(x,y)$ that are not scale-invariant but collapse into a linear
function at $f(x,x)$. Although microscopic theories based on CGC suggest
reduced thickness function with scale-invariance property (which we will
discuss shortly), there is still room for exploration regarding the scale-
invariance in non-central collisions and the case where $N_{A}\neq N_{B}$ from
an experimental point of view.

It is a known approximation that the final multiplicity is proportional to the
initial entropy, $N_{\text{ch}}\propto S$. Also, by definition, the
participant thickness functions,

$T_{A(B)}(\mathbf{x})=\kappa\sum_{i=1}^{N_{A(B)}}\int
dz\,\rho_{\text{nucleon}}(\vec{x}-\vec{x}_{i}),$ (5)

scales similar to the number of participant $N_{A(B)}$ where $\mathbf{x}$ is
the coordinate in the transverse space. Here, we have assumed $\int
d^{3}x\rho_{\text{nucleon}}$ is dimensionless and normalized to unity, and
$\kappa$ is a dimensionless constant for overall normalization. Therefore, by
defining the reduced thickness function as the entropy density in the
transverse direction, $T_{R}(\mathbf{x})=\int dz\,s(\vec{x})$, one suggests a
differential form of Eq. (4) as $T_{R}(T_{A},T_{B})\propto f(T_{A},T_{B})$. In
Eq. (1), an specific form of $f(T_{A},T_{B})$ is chosen.

It is important to note that the argument mentioned above assumes that the
system has reached its local thermal equilibrium from its initial state. In
order to extend the argument to the out-of-equilibrium period of the
evolution, it is necessary to establish a connection between the entropy
density, $s_{\text{hyd}}$, at the start of local thermal equilibrium,
$\tau_{\text{hyd}}$, and the energy density, $e_{0}$, at an earlier time,
$\tau_{0}\gtrsim 0^{+}$. In Ref. [29], the concept of hydrodynamic attractors
was used to obtain,
$s_{\text{hyd}}\tau_{\text{hyd}}\propto(e_{0}\tau_{0})^{2/3}$ where the right-
hand side is in the out-of-equilibrium period. This formula can be used to
calculate the final multiplicity, which is given by

$\frac{dN_{\text{ch}}}{d\eta}\propto\frac{dS}{d\eta}=A_{\perp}(s_{\text{hyd}}\tau_{\text{hyd}})\propto
A_{\perp}(e_{0}\tau_{0})^{2/3},$ (6)

where $A_{\perp}$ is the transverse size of the system at the time
$\tau_{\text{hyd}}$. At the pre-equilibrium stage, we anticipate a slight time
dependence for $A_{\perp}$. To achieve scale-invariance in Eq. (4), the energy
reduced thickness function needs to scale through a homogeneous function with
degree 3/2,

$T_{R}^{\text{energy}}(T_{A},T_{B})=\mathcal{N}\frac{\tau_{\text{hyd}}^{3/2}}{\tau_{0}}f^{3/2}(T_{A},T_{B}).$
(7)

Assuming $T_{A(B)}$ as the energy thickness function, the constant $\kappa$ in
Eq. (5) has a dimension of energy too.

In the context of the CGC framework, Ref. [30] obtains the initial energy
density using the saturation scale $Q_{s,A(B)}^{2}$. Since
$Q_{s,A(B)}^{2}\propto T_{A(B)}$, the reduced thickness function can be
defined as,

$\begin{split}T_{R}^{\text{energy}}&(T_{A},T_{B})\propto\\\
&\frac{T_{A}T_{B}}{(T_{A}+T_{B})^{5/2}}\left[2T_{A}^{2}+7T_{A}T_{B}+2T_{B}^{2}\right].\end{split}$
(8)

The scaling of $T_{A,B}\to c\,T_{A,B}$ leads to $T_{R}^{\text{energy}}\to
c^{3/2}T_{R}^{\text{energy}}$, and consequently, $(dN_{\text{ch}}/d\eta)\to
c\,(dN_{\text{ch}}/d\eta)$, indicating scale-invariance.

In Ref. [22], a modified version of reduced thickness function is introduced,

$T_{R}^{\text{energy}}(T_{A},T_{B})=\mathcal{N}E_{\text{ref}}^{2-2q}\left[(T_{A}^{p}+T_{B}^{p})/2\right]^{q/p},$
(9)

which violates the scale-invariance for $q\neq 3/2$. This definition with
$q=3/2$ and $E_{\text{ref}}=\tau_{0}(\tau_{\text{hyd}}\kappa)^{-3/2}$ is
equivalent to Eq. (7). The global Bayesian analysis in Ref. [22] suggests
$q=1.34\begin{subarray}{c}+0.14\\\ -0.18\end{subarray}$ which is marginally
compatible with $q=3/2$. In the mentioned study, the relation $e\propto
s^{4/3}$ is used to convert entropy to energy density, which is different from
our approach based on Eq. (6). Additionally, the study analyzes the scaling
property of the IP-Glasma and shows that $e_{0}\propto(T_{A}T_{B})^{q/2}$,
where the value of $q$ falls in the range of 0.54 to 1.71.

In summary, scale-invariance refers to the case that the initial state scaling
results a same amount of scaling in the final multiplicity. If the initial
entropy density is used, the reduced thickness function need to be a
homogeneous function of degree one. If the energy density is used, the degree
of homogeneity should be 3/2. For this study, we assume that the initial state
is modeled using the entropy density.

Figure 1: The convergence of initial entropy density moments is affected by
impact parameter and collision system symmetry. When the impact parameter is
zero, and the collision system is symmetric, the initial state moments
converge to the same curve as $w$ increases. The moments do not converge to
the same curve for events with non-zero impact parameters or asymmetric
collision systems.

To investigate the relationship between the nucleon size, impact parameter,
and asymmetric collisions, as well as the impact of scale-invariance on the
moments of the initial entropy density in UCSC, we use
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model and
the model based on CGC, which have scale-invariance properties. We also use
the MC-Glauber model and the generalized
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model to
examine scenarios with broken scale-invariance. The MC-Glauber model, CGC
model, and generalized
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model are
implemented in
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ event
generator via Eqs. (3), (8), and (9). To quantitatively study $T_{R}$, we
calculate the following moments,

$\\{r^{2}\\},\qquad\epsilon_{n}e^{in\phi_{n}}=-\frac{\\{r^{n}e^{in\varphi}\\}}{\\{r^{n}\\}},\quad
n=2,3,$ (10)

where $(r,\varphi)$ is the polar coordinate in the transverse space
$\mathbf{x}$, and $\\{\cdots\\}=\int d\mathbf{x}\cdots T_{R}(\mathbf{x})$.
Using T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$
event generator, we generate Au–Au collisions with impact parameters
$b=0\,$[fm] and $4\,$[fm], and Au–Cu with $b=0\,$[fm] at
$\sqrt{s_{\text{NN}}}=200\,$GeV as an example of an asymmetric collision.
Here, the fluctuation parameter and minimum distance between nucleons are set
as follows: $\sigma_{\text{fluc}}=1$, $d_{\text{min}}=0.6\,$[fm]. We calculate
the following many event average quantities, $\left\langle
dS/d\eta\right\rangle$, $\epsilon_{2}\\{2\\}$, and $\langle\\{r^{2}\\}\rangle$
and plot them in terms of nucleon width $w$, shown in Fig. 1. As seen from the
figure, when the collision system is symmetric, the moments for different
values of $p$ tend to converge to the same values for large $w$ during most
central collisions while it is not the case for non-central collisions and
also for asymmetric collisions, meaning the functionality of reduced thickness
function has a small impact on ultra-central symmetric collisions. In the
appendix A, we examined $\epsilon_{3}\\{2\\}$,
$\langle\epsilon_{2}^{2}\epsilon_{3}^{2}\rangle-\langle\epsilon_{2}^{2}\rangle\langle\epsilon_{3}^{2}\rangle$,
and
$\langle\epsilon_{2}^{2}\epsilon_{4}\cos\left(4\varphi_{2}-4\varphi_{4}\right)\rangle$
which also shows similar behavior.

In appendix B, we study the effect of impact parameter for $p=0$ and $p=1$ in
more detail for two limits $\pi w^{2}\gg\pi R_{0}^{2}/A$ and $\pi w^{2}\ll\pi
R_{0}^{2}/A$. We summarized our result in table in the appendix B. For large
$w$, the total entropy for
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ reduced
thickness function in Eq. (1) with $p=1$ is proportional to $(N_{A}+N_{B})/2$
while for $p=0$ it is given by $\sqrt{N_{A}N_{B}}$. These two values approach
the same value for the symmetric collisions $N_{A}\approx N_{B}$. Another
example is the rms radius $\\{r^{2}\\}$. For $p=1$, we find

$\begin{split}\left.\\{r^{2}\\}\right|_{p=1}=&\frac{N_{A}}{N_{A}+N_{B}}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}\\\
&\hskip
28.45274pt+\frac{N_{B}}{N_{A}+N_{B}}\\{x_{b}^{2}+y_{b}^{2}\\}_{s}+2w^{2},\end{split}$
(11)

and for $p=0$,

$\begin{split}&\left.\\{r^{2}\\}\right|_{p=0}=\frac{1}{2}\Big{(}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\\{x_{b}^{2}+y_{b}^{2}\\}_{s}+2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}\\\
&\hskip 99.58464pt+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}\Big{)}+2w^{2},\end{split}$
(12)

where $\mathbf{x}_{a(b)}$ is the location of the target (projectile) sources
and $\\{\cdots\\}_{s}$ refers to the average over sources location. For the
central collisions, we have
$\\{\mathbf{x}_{a}\\}_{s}\approx\\{\mathbf{x}_{b}\\}_{s}\approx 0$ and for
symmetric collisions $N_{A}\approx N_{B}$, leading to
$\left.\\{r^{2}\\}\right|_{p=1}\approx\left.\\{r^{2}\\}\right|_{p=0}$. We find
similar behavior for $\\{y^{2}-x^{2}\\}$ and $\\{2xy\\}$.

To examine whether this convergence is the result of the scale-invariance, we
have plotted $\epsilon_{2}\\{2\\}$ for Au–Au collisions where the reduced
thickness function are calculated from MC-Glauber model ($\alpha=0.28$), CGC
model, and (generalized)
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model for
$p=0$ with $q=1,1/2,2$. To translate the CGC model with initial energy density
into scale-invariant entropy density, $T_{R}^{\text{energy}}(T_{A},T_{B})$ is
raised to the power 2/3. Also, Eq. (9) is assumed to be entropy density. In
this scenario, $q=1$ represents the scale-invariant case, not $q=3/2$. We have
also assumed $E_{\text{ref}}=1\,$[fm-1]. The parameter $\alpha$ in the MC-
Glauber model is intentionally chosen to be larger than usual to emphasize the
violation of scale invariance. As observed in Fig. 1, the CGC model and
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model with
$q=1$ converge to the same values, while the MC-Glauber and
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model with
$q=2$ and $1/2$ do not, indicating the violation of scale invariance in the
latter models. The same behavior is observed for $\epsilon_{3}\\{2\\}$ as
shown in appendix A. It is evident that a small violation of scale-invariance
in the MC-Glauber model leads to a small convergence violation. This suggests
that a small contribution of the reduced thickness function definition in UCSC
is valid also in cases where scale-invariance is only approximately
maintained.

##  III Many event average initial state in ultra-central symmetric collisions

In this section, we study the effect of short-range correlations and
constituent weight fluctuation on the event-by-event source fluctuations. The
technical points and details of this section can be found in the appendix C.

###  III.1 Two-body correlations from a statistical mechanics analogy

The position of the sources is determined by the wavefunction of the colliding
nuclei. Specifically, the one-body density, denoted as $p_{A}(\vec{x})$, and
the two-body correlation, denoted as $c_{2}(\vec{x}_{1},\vec{x}_{2})$, have a
significant impact on the measurement of anisotropic flow. To establish a
relationship between the anisotropic flow observables and $p_{A}(\vec{x})$ and
$c_{2}(\vec{x}_{1},\vec{x}_{2})$ at UCSC, we make use of the approximation
discussed in Section II. In this scenario, the reduced thickness function is
assumed to be the sum of $T_{A}(\mathbf{x})$ and $T_{B}(\mathbf{x})$. It can
be demonstrated that in this case, the moments of $T_{R}$ can be divided into
two parts: one part is influenced by the location of the sources, and the
other by the distribution of the constituents, such as nucleons or quarks
(appendix C.1). This can be observed in Eqs. (11) and (12), where the
contribution of $2w^{2}$ is separated from the average source location. In
UCSC and assuming a sufficiently large nucleon-nucleon cross-section, it is
assumed that all nucleons participate in the collision. This implies that the
location of the initial state sources must be similar to those in each nucleus
at the time of the collision. Consequently, the fluctuation of sources in a
single nucleus can be translated into fluctuations of the entire initial state
(appendix C.2).

The process of sampling the nucleon sources becomes more complicated when
considering correlations, especially in the context of deformed nuclei. In
subsection III.2, we will explore this further. The approximation introduced
in UCSC allows us to take a more direct approach to understanding the effect
of short-range correlations by using classical techniques from statistical
mechanics. We will explore this concept in the following.

The cornerstone concept of classical statistical mechanics is that to obtain
thermodynamic quantities in an equilibrated system, the time average over
physical quantities can be replaced by averaging over an ensemble of systems.
The state of each system is distributed based on statistical distributions,
mostly depending on variables such as temperature or chemical potential. In
the study of fluctuation in heavy-ion physics, the ensemble of systems has a
clear counterpart, which is the collection of independent events. The concept
of an equilibrated system is translated to the fact that all events in heavy-
ion collisions are equivalent to each other and are achieved under the same
experimental conditions, such as the center of mass energy or the same species
collisions. The “state” of an event is the position of the sources, which are
distributed based on the nucleus wavefunction. To make the analogy clearer,
assume $P_{A}(\vec{x}_{1},\ldots,\vec{x}_{A})$ as the distribution of sources
in a nucleus,

$P_{A}(\vec{x})=\frac{1}{\mathcal{Z}_{A}}F_{A}(\vec{x}),$ (13)

where $\mathcal{Z}_{A}$ is $P_{A}$ normalization factor. If sources are
distributed based on a Gaussian profile with width $\sigma$, one needs to
define
$F_{A}(\mathbf{x})=e^{-\frac{1}{2\sigma^{2}}\sum_{i=1}^{A}|\mathbf{x}_{i}|^{2}}$,
analogous to a free non-relativistic gas with $F_{N}(\vec{p})=e^{-\beta
E}=e^{-\frac{1}{2mkT}\sum_{i=1}^{N}|\vec{p}_{i}|^{2}}$.

This analogy becomes particularly useful when attempting to account for short-
range correlations. The cluster expansion method is utilized to determine
thermodynamic quantities in the presence of short-range correlation (see Ref.
[31]). To apply this method to our problem, we define

$F_{A}(\vec{x})=\left[\prod_{i=1}^{A}p_{A}(\vec{x}_{i})e^{\frac{1}{A}\sum_{a}\lambda_{a}f_{a}(\vec{x}_{i})}\right]\left[\prod_{1\leq
i<j\leq A}c_{2}(\vec{x}_{j},\vec{x}_{j})\right].$ (14)

In statistical mechanics, the partition function typically serves as the
generating function. The parameter $\beta$ is used to extract the
thermodynamic variables and their fluctuations. In our approach, we introduce
an auxiliary term $e^{\frac{1}{A}\sum_{a}\lambda_{a}f_{a}(\vec{x}_{i})}$,
where $\lambda_{a}$ plays a similar role to $\beta$, but it needs to be set to
zero once the desired average is obtained by differentiating the partition
function. In order to find the desired moment, it is necessary to select an
appropriate function $f_{a}(\vec{x})$.

The cluster expansion provides a closed form for the grand-partition function
(details in appendix C.3),

$\begin{split}\log\mathscr{Q}&(\langle A\rangle,\lambda_{a})=\\\
b_{1}(\lambda_{a})&\;\langle
A\rangle+\left[b_{2}(\lambda_{a})-2b_{2}(0)b_{1}(\lambda_{a})\right]\;\langle
A\rangle^{2}+\cdots.\end{split}$ (15)

The $\ell$-cluster integrals, $b_{\ell}(\lambda_{a})$, connect partition
function to $p_{A}(\vec{x})$ and $c_{2}(\vec{x}_{1},\vec{x}_{2})$,

$\displaystyle b_{1}(\lambda_{a})$ $\displaystyle=\int
d^{3}x\;p_{A}(\vec{x})\;e^{\frac{1}{A}\sum_{a}\lambda_{a}f_{a}(\vec{x})},$
(16a) $\displaystyle b_{2}(\lambda_{a})$ $\displaystyle=\frac{1}{2}\int
d^{3}x_{1}d^{3}x_{2}\;p_{A}(\vec{x}_{1})p_{A}(\vec{x}_{2})\;$ (16b)
$\displaystyle\hskip
85.35826pte^{\frac{1}{A}\sum_{a}\lambda_{a}(f_{a}(\vec{x}_{1})+f_{a}(\vec{x}_{2}))}\;f_{12}.$

where the so-called Mayer function is defined as
$f_{12}=c_{2}(\vec{x}_{1},\vec{x}_{2})-1$. Here, we assume that $\langle
A\rangle\equiv A$ is large enough that we can disregard the fluctuation of $A$
and approximate $\mathscr{Q}(\langle A\rangle,\lambda_{a})\propto
Q_{A}(\lambda_{a})$, where $Q_{A}(\lambda_{a})$ is the canonical partition
function at a fixed value of $A$. It’s important to note that the calculated
observables from canonical and grand-canonical partition functions might have
differences in subleading orders in the $1/A$ expansion, which are ignored in
the present study.

To find moments such as
$\langle\\{f_{1}\\}^{m_{1}}\\{f_{2}\\}^{m_{2}}\rangle$, one should calculate
$\partial_{\lambda_{1}}^{m_{1}}\partial_{\lambda_{2}}^{m_{2}}\mathscr{Q}(\langle
A\rangle,\lambda_{a})/\mathscr{Q}(\langle A\rangle,\lambda_{a})$ and set
$\lambda_{a}=0$. To obtain $\langle\\{r^{2}\\}\rangle$, we choose
$f_{1}(\vec{x})=|\mathbf{x}|^{2}$, $m_{1}=1$, and $m_{2}=0$. The functions
$f_{1}(\vec{x})=f_{2}^{*}(\vec{x})=|\mathbf{x}|^{2}e^{2i\varphi}$ with
$m_{1}=m_{2}=1$ corresponds to the moment average
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$.

The method described above is general. For a specific choice of
$p_{A}(\mathbf{x})$ as $\rho_{\text{WS}}(\mathbf{x})$ with radius
$\tilde{R}_{0}$ and skin thickness $\tilde{a}_{0}$, and
$c_{2}(\mathbf{x}_{1},\mathbf{x}_{2})=\Theta(|\mathbf{x}_{1}-\mathbf{x}_{2}|-d_{\text{min}})$,
which corresponds to nucleons’ hard inner core with a radius of
$d_{\text{min}}/2$, we obtain an analytical estimate for the radius as

$\begin{split}\left\langle\\{r^{2}\\}_{s,A}\right\rangle=&\frac{2\tilde{R}_{0}^{2}}{5}\left[1+\frac{7\pi^{2}}{3}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{2}\right.\\\
&\left.\quad+\frac{7\pi^{2}}{3}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{2}\left(\frac{d_{\text{min}}}{\tilde{R}_{0}}\right)^{3}+\cdots\right].\end{split}$
(17)

Here, $\langle\\{\cdots\\}_{s,A}\rangle$ stands for the average over the
position of nucleons in a single nucleus $A$. The elliptic shape of the
initial density can also be found as

$\begin{split}&\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle=\frac{1}{A}\frac{8\tilde{R}_{0}^{2}}{35}\left[1+6\pi^{2}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{2}\right.\\\
&\left.\hskip
28.45274pt+\frac{31\pi^{4}}{3}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{4}+d^{\prime}_{\varphi}\left(\frac{d_{\text{min}}}{\tilde{R}_{0}}\right)^{3}+\cdots\right],\end{split}$
(18)

where:

$d^{\prime}_{\varphi}=-A\left[1-6\pi^{2}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{2}-\frac{31\pi^{4}}{3}\left(\frac{\tilde{a}_{0}}{\tilde{R}_{0}}\right)^{4}\right].$
(19)

When calculating the expression $d^{\prime}_{\varphi}$, we replaced WS
distribution with a step function $\Theta(\tilde{R_{0}}-|\vec{x}|)$ to
facilitate the analytical integration. We compared our analytical estimate for
$d^{\prime}_{\varphi}$ with the numerical result and found that the numerical
value of $d^{\prime}_{\varphi}$ is approximately half of the value obtained
using the step function, indicating the significance of the skin thickness in
determining this quantity. We confirmed our calculation using Monte Carlo
sampling in sections III.2 and IV.

###  III.2 Note on the one-body density and two-body correlation and Monte
Carlo event sampling

Short-range correlation changes the one-body distribution. For example,
nucleons with a hard central core will be pushed into the outer regions,
altering the one-body density assumed by the WS distribution. The two-body
density is given by

$\begin{split}&p_{A}^{(2)}(\vec{x}_{1},\vec{x}_{2})=\\\
&\frac{\rho_{\text{WS}}(\vec{x}_{1})\rho_{\text{WS}}(\vec{x}_{2})\Theta(|\vec{x}_{1}-\vec{x}_{2}|-d_{\text{min}})}{\int
d\vec{x}^{\prime}_{1}d\vec{x}^{\prime}_{2}\rho_{\text{WS}}(\vec{x}^{\prime}_{1})\rho_{\text{WS}}(\vec{x}^{\prime}_{2})\Theta(|\vec{x}^{\prime}_{1}-\vec{x}^{\prime}_{2}|-d_{\text{min}})},\end{split}$
(20)

where $\Theta(x)$ is the Heaviside step-function. It is important to note that
$p_{A}(\vec{x}_{1})=\int d\vec{x}_{2}p_{A}^{(2)}(\vec{x}_{1},\vec{x}_{2})$ is
not coincident with the WS distribution, although we might approximate it with
the WS distribution with modified $R_{0}$ and $a_{0}$ at small values of
$d_{\text{min}}$ (see section. III.3). The two-body correlation also deviates
from the step-function as discussed in Ref. [32].

To overcome the problem of one-body density modification, in
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ event
generator, the radial direction of the sources is sampled first based on WS
distribution to keep the radial size of the nucleus fixed. Afterward, the
spherical angles of the sources are sampled iteratively. Sources are sorted in
the radial direction and moving from smaller radii, the $i$th source is placed
at $\vec{x}_{i}$ if its distance from all other $i-1$ sources is greater than
$d_{\text{min}}$, otherwise the $i$th source angel are rejected and resample
the angles until it is placed somewhere in the space. Although the pre-
sampling of the radial direction is performed to keep the WS structure under
control, it systematically leads to inaccuracy in placing the nucleons with
minimum distance. Considering two nucleons are sampled at $r_{1}$ and $r_{2}$
radial direction, the maximum distance of these two nucleons is
$|r_{1}-r_{2}|$, meaning if this distance is smaller than $d_{\text{min}}$, it
is impossible to place the second source at any given spherical angles. In
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ event
generator, the angles are sampled 1000 times and if it fails, the source
leaves it at its last iteration. The failure percentage is negligible for
small values of $d_{\text{min}}$, but this systematic inconsistency should be
considered as a hint that this method may lead to an inaccurate two-body
correlation. More importantly, for deformed nuclei, the nucleus radius depends
on the orientation, and the method employed by
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ event
generator is no longer applicable.

There is also a caveat in placing the sources iteratively. One may assume that
at each iteration, all three directions are generated based on the WS
distribution, and the minimum distances with former generated sources are
examined iteratively. However, in the iterative method starting from WS with
no information from short-range behavior, the $i$th source generated in the
$i$th iteration is only correlated with all $(i-1)$ sources generated in the
previous iterations. The first source, in fact, contains no information from
short-range correlation. On the other hand, if we start from $p_{A}(\vec{x})$
and $c_{2}(\vec{x}_{1},\vec{x}_{2})$, the information of the correlation is
already encoded in the sampling from the first iteration. We have examined
this claim in a 1D toy model with only two sources, where each source is
generated based on a Gaussian distribution, and the correlation is imposed
using minimum distance. We can compare our toy Monte Carlo with analytical
calculations using this simplified toy model. The study mentioned in the
appendix E suggests that the iterative method works correctly only if we use
the true one-body density that contains the short-range effect information.
This indicates that in more realistic cases, the iterative method works well
if we use the one-body density and two-body correlation extracted from Eq.
(20). Short-range correlations can affect the radius of spherically symmetric
nuclei. In the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model, the
radius is kept fixed, so one considers that the short-range correlation is
accounted for in this method, although pre-generating the radius can cause
issues with source placement and not applicable for deformed nuclei.

Figure 2: Comparing three different methods of source sampling of point-like
sources in one Au nucleus at the presence of short-range correlation.

An alternative way of generating sources using WS distribution and minimum
distance among nucleons is to generate $A$ random sources based on WS and then
compare all the mutual distances. If there is any pair with condition
$|\vec{x}_{1}-\vec{x}_{2}|<d_{\text{min}}$, the whole set of $A$ sources
should be discarded, and a new set is generated. We call this method as direct
method. In this way, we find consistent one-body density (which is different
from the original WS) and two-body correlation. There is no problem with the
position of the source when using this method, and it can be used for any
deformed nuclei. However, it should be noted that this method takes longer to
execute than the iterative method, especially for larger values of
$d_{\text{min}}$. In Ref. [33], the sources with appropriate two-body
correlation are generated in a box using Metropolis random search, and then
the WS probability distribution is imposed to choose those sources relevant to
the one-body density. This method has been employed by SMASH initial state
model in Ref. [34]. We have examined a similar strategy in our 1D toy model by
generating sources with some minimum distances and then imposing the Gaussian
constraint. The result is coincident with direct method and leads to correct
one-body density and two-body correlation. There is also a recent method (we
call it as shifting method) in imposing the two-body correlations by
generating independent sources first based on WS distribution and then
shifting away the sources using a carefully constructed function based on the
two-body correlation [32].

In Fig. 2, we compare three different methods: direct method, iterative
method, and the shifting method, to appreciate their differences. The
iterative method generates all three coordinates of the sources, not just the
angular part, which makes it different from the method implemented in
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$. To
eliminate the complications arising from entropy production and nucleon
participation, we consider only one Au nucleus to generate point-like sources
and impose a two-body correlation using a step-function. We calculate
$\langle\\{r^{2}\\}_{s,A}\rangle$ and
$\langle\left|\\{r^{2}e^{i2\varphi}\\}_{s,A}\right|^{2}\rangle$ for different
values of $d_{\text{min}}$. In the direct method, we did not go beyond
$d_{\text{min}}=0.4\,$[fm] due to the inefficiency of the direct method. The
figure shows that all three sets of events give the same values when there are
zero short-range correlations. However, for $d_{\text{min}}>0$, the main
difference appears in $\langle\\{r^{2}\\}_{s,A}\rangle$, which quantifies the
size of the one-body density. Short-range correlation modifies the one-body
density, but the modification from the three methods is different. The dashed
(solid) lines are analytical (numerical) estimations of
$d^{\prime}_{\varphi}$. The cluster expansion calculation is compatible with
direct method of source sampling. The difference in
$\langle\left|\\{r^{2}e^{i2\varphi}\\}_{s,A}\right|^{2}\rangle$ is less
pronounced. The direct method and cluster expansion calculation rely on using
consistent one-body density and two-body correlation. Iteratively sampling
sources without considering the short-range correlation into the one-body
density can lead to inconsistency. The authors still need to investigate in
greater detail the observed difference between the direct method and the
shifting method.

###  III.3 Radius and skin thickness redefinition

As seen from Eq. (17), the short-range correlation modifies the size of the
nucleus as expected. The physical values of the radius and skin thickness
contain the effect of the short-range correlation. For that reason, we
redefine “bare” values of $\tilde{R}_{0}$ and $\tilde{a}_{0}$ to absorb the
term corresponding to the correlation in
$\left\langle\\{r^{2}\\}_{s,A}\right\rangle$ and write the final result in
terms of “physical” values $R_{0}$ and $a_{0}$. Using Eq. (17), we cannot fix
both radius and skin thickness. We assume two cases: first, the skin thickness
is fixed, $\tilde{a}_{0}=a_{0}$, but the radius is changed, and second, the
radius is fixed, $\tilde{R}_{0}=R_{0}$, but the skin thickness is modified.
For the former case, we find the value of $\tilde{R}_{0}$ in terms of physical
quantities $R_{0}$ as follows:

$\begin{split}\tilde{R}_{0}^{2}&=R_{0}^{2}\left[1-A\frac{7\pi^{2}}{3}\left(\frac{a_{0}}{R_{0}}\right)^{2}\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}+\cdots\right],\end{split}$
(21)

while $\tilde{a}_{0}$ for the latter reads as

$\begin{split}\tilde{a}_{0}&=a_{0}\left[1-\frac{A}{2}\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}+\cdots\right].\end{split}$
(22)

Figure 3: The Au nucleus density sampled with $R_{0}=6.38\,$[fm],
$a_{0}=0.535\,$[fm], and $d_{\text{min}}=0.4$[fm] over WS distribution with
the same radius and skin thickness (blue), larger radius fixed skin thickness
(red, dashed), fixed radius larger skin thickness (green, dot-dashed). The
sampled density is close to WS with a larger radius but fixed skin thickness.

In Fig. 3, we show the modification of WS distribution in the presence of
correlation. Here, we sampled a 200 k Au nucleus with
$\tilde{R}_{0}=6.38\,$[fm], $\tilde{a}_{0}=0.535\,$[fm], and
$d_{\text{min}}=0.4\,$[fm]. We divided the radial direction in the transverse
direction into 50 bins and calculated the average number of sources in each
bin. This average represents the density of the nucleus in the presence of the
correlation. Then, we compared the density from the sampled sources with the
original WS distribution. As seen from the blue curve in the figure, the
density from the sampled sources differs from the WS distribution, which has
the same radius and skin thickness. The density ratio goes above unity at the
tail of the distribution, indicating that sources are pushed away due to the
short-range correlation. We also calculated $R_{0}$ and $a_{0}$ from Eqs. (21)
and (22) and sampled distribution using redefined values. We observe that WS
with a modified radius and fixed skin thickness can explain the new
distribution while keeping $R_{0}$ fixed, but modifying $a_{0}$ does not match
the sampling. In the present study, we follow the first strategy and modify
the radius such that the final nucleus radius matches the physical radius.

###  III.4 Constituent weight fluctuation

In addition to the fluctuation in the position of the sources, the
contribution of each source to the initial entropy density can fluctuate. Ref.
[20] discusses that the final multiplicity fluctuation in p–Pb collision can
be described by assuming that the participant contribution to the initial
entropy density fluctuation follows a Gamma distribution,

$p_{\text{fluc}}(\gamma)=\frac{k^{k}}{\Gamma(k)}\gamma^{k-1}e^{-k\gamma},\qquad\sigma_{\text{fluc}}=1/\sqrt{k}.$
(23)

To take the effect of fluctuation into account, we replace
$p_{A}(\vec{x}_{i})$ with $p_{\text{fluc}}(\gamma_{i})p_{A}(\vec{x}_{i})$ and
modify the auxiliary function to
$e^{\frac{1}{A}\sum_{a}\lambda_{a}\gamma_{i}f_{a}(\vec{x}_{i})}$ in Eq. (14).
We have examined the contribution of
$\sigma_{\text{fluc}}^{2}(d_{\text{min}}/R_{0})^{3}$ in moment
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle$.
This contribution is zero for the case WS distribution is replaced by step-
function. Simulation based on
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model also
confirms that this contribution is small. Therefore, we calculate the
contribution of $\sigma_{\text{fluc}}$ and $d_{\text{min}}$ separately. By
ignoring the effect of two-body correlation, we find that the canonical
partition function is written in terms of 1-cluster integral,
$Q_{A}(\lambda)\propto b_{1}^{A}(\lambda)$. Employing
$\partial_{\lambda}Q_{A}(\lambda)/Q_{A}(\lambda)|_{\lambda=0}$ leading to

$\langle\\{f\\}_{s,A}\rangle=\int d^{3}x\,f(\vec{x})\,p_{A}(\vec{x}).$ (24)

In the case of calculating second-order moments average, from
$\partial_{\lambda_{1}}\partial_{\lambda_{2}}Q_{A}(\lambda)/Q_{A}(\lambda)|_{\lambda_{i}=0}$,
we find

$\begin{split}\langle\\{f_{1}\\}_{s,A}\\{f_{2}\\}_{s,A}\rangle&=\frac{1+\sigma_{\text{fluc}}^{2}}{A}\int
d^{3}x\,f_{1}(\vec{x})f_{2}(\vec{x})\,p_{A}(\vec{x})\\\
&+\frac{A-1}{A}\langle\\{f_{1}\\}_{s,A}\rangle\langle\\{f_{2}\\}_{s,A}\rangle,\end{split}$
(25)

where we have used the fact that $\int
d\gamma\,\gamma\,p_{\text{fluc}}(\gamma)=1$ and $\int
d\gamma\,\gamma^{2}\,p_{\text{fluc}}(\gamma)=1+\sigma_{\text{fluc}}^{2}$.
Approximating $Q_{A}(\lambda)$ by its grand-canonical counterpart, as
discussed in section III.1, would have led to replacing $A-1$ with $A$ in the
second term of Eq. (25). This modification is negligible for large $A$ and
ignored here.

To find the values of $\langle\\{r^{2}\\}\rangle$ and
$\langle\\{r^{2}\\}^{2}\rangle$, we need to substitute
$f(\vec{x})=|\mathbf{x}|^{2}$ and
$f_{1}(\vec{x})=f_{2}(\vec{x})=|\mathbf{x}|^{2}$ into Equations (24) and (25)
respectively. It is evident from Equation (25) that the second term
contributes the most when $A$ is large, which implies that
$\langle\\{r^{2}\\}^{2}\rangle\approx\langle\\{r^{2}\\}\rangle^{2}$. This
shows that the value of $\langle\\{r^{2}\\}^{2}\rangle$ is mainly influenced
by the size and has a small contribution from fluctuation. On the other hand,
to determine
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle$, we
choose $f_{1}(\vec{x})=f_{2}^{*}(\vec{x})=|\mathbf{x}|^{2}e^{2i\varphi}$. In
this case, the second term in Equation (25) is zero, which means that the
primary contribution is proportional to $(1+\sigma_{\text{fluc}}^{2})/A$,
indicating that it solely originates from the fluctuations.

In Equation (25), the presence of $\sigma_{\text{fluc}}^{2}$ arises from the
calculation of second order moments. Therefore, the dependence on
$\sigma_{\text{fluc}}^{2}$ can be accurate for values greater than unity. When
there is fluctuation present, we use the expression
$\\{f\\}_{s}=\sum_{i}f(\vec{x}_{i})/N$ instead of
$\\{f\\}_{s}=\sum_{i}\gamma_{i}f(\vec{x}_{i})/\sum_{i}\gamma_{i}$ to calculate
many event averages. The approximation we make is that the numerator and
denominator fluctuate independently event by event. We have provided a more
detailed discussion in Appendix C. We validated this approximation for
$\sigma_{\text{fluc}}$ up to 2 using a toy Monte Carlo.

###  III.5 From source fluctuation of a single nucleus to initial density
fluctuation

So far, we have calculated the impact of source fluctuation on a single
nucleus, $A$. The first-order moment for combining sources from two nuclei is
directly related to the single nucleus moment,
$\langle\\{f\\}_{s}\rangle=\langle\\{f\\}_{s,A}\rangle$. For a general second-
order moment, we find

$\begin{split}\langle\\{f_{1}\\}_{s}\\{f_{2}\\}_{s}\rangle&=\frac{1}{2}\left[\left\langle\\{f_{1}\\}_{s,A}\\{f_{2}\\}_{s,A}\right\rangle+\right.\\\
&\hskip
56.9055pt\left.\left\langle\\{f_{1}\\}_{s,A}\right\rangle\left\langle\\{f_{2}\\}_{s,A}\right\rangle\right].\end{split}$
(26)

In addition, at a finite cross-section, there are some sources that do not
participate in the interaction. This leads to a correction in the
decomposition of fluctuations into independent target and projectile.
Moreover, for the deformed nuclei in UCSC, an intrinsic value is introduced to
the elliptic shape, and their random orientation leads to the fluctuation [6].
The details of these corrections, along with the effect of nucleon
substructure with $n_{c}$ constituents and width $v$, are discussed in
appendix F. The final result for the initial state size is given as follows:

$\begin{split}&\langle\\{r^{2}\\}\rangle^{2}=\frac{4R_{0}^{4}\,\omega^{2}}{25}\Bigg{[}\omega^{2}+\alpha_{r}\left(\frac{a_{0}}{R_{0}}\right)^{2}\\\
&+s_{r}\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right)+b_{2,r}\beta_{2}^{2}+b_{3,r}\beta_{3}^{2}+\cdots\Bigg{]}.\end{split}$
(27)

Accordingly, the elliptical shape is obtained as

$\begin{split}&\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle=\frac{1}{A}\frac{4R_{0}^{4}}{35}\left[1+\alpha_{\varphi}\left(\frac{a_{0}}{R_{0}}\right)^{2}\right.\\\
&\hskip
48.36958pt\left.+s_{\varphi}\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right)+d_{\varphi}\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}\right.\\\
&\hskip
48.36958pt\left.+v_{\varphi}\frac{\sigma_{\text{fluc}}^{2}}{n_{c}}+b_{2,\varphi}\beta_{2}^{2}+\cdots\right].\end{split}$
(28)

The coefficients are shown in table 1. We note that $(w^{2}-v^{2})/R_{0}^{2}$
is the expansion parameter, and it should be small, while $(w/R_{0})^{2}$ is
not. In fact, the estimation is more accurate for larger values of $w$.

coefficient | value
---|---
$\omega^{2}$ | $1+5\left(w/R_{0}\right)^{2}$
$\alpha_{r}$ | $14\pi^{2}/3$
$\alpha_{\varphi}$ | $6\pi^{2}$
$s_{r}$ | $-6$
$s_{\varphi}$ | $51/4$
$d_{\varphi}$ | $d^{\prime}_{\varphi}-\alpha_{r}A\left(a_{0}/R_{0}\right)^{2}$
$v_{\varphi}$ | $1+10\left(w^{2}-v^{2}\right)/R_{0}^{2}+\alpha_{\varphi}\left(a_{0}/R_{0}\right)^{2}$
$b_{2,r}$ | $29/8\pi$
$b_{3,r}$ | $7/2\pi$
$b_{2,\varphi}$ | $21\,A/20\pi$
$\alpha_{\epsilon}$ | $\alpha_{\varphi}-\alpha_{r}/\omega^{2}$
$s_{\epsilon}$ | $s_{\varphi}-s_{r}/\omega^{2}$
$b_{2,\epsilon}$ | $b_{2,\varphi}-b_{2,r}\omega^{2}$
$b_{3,\epsilon}$ | $-b_{3,r}/\omega^{2}$

Table 1: Table of coefficients for $\langle\\{r^{2}\\}\rangle$,
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$, and
$\epsilon_{2}\\{2\\}$ estimations.

##  IV Ellipticity fluctuation and isobar ratio

Figure 4:  Comparing the analytical estimation (curves) for
$\langle\\{r^{2}\\}\rangle$,
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$, and
$\epsilon_{2}\\{2\\}$ with
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
outcome (dots) for Au–Au collisions with
$\sigma_{\textbf{inel}}^{NN}=4.2\,$[fm2]. The compatibility improves for
larger nucleons size due to the improvement of approximation based on scale
invariance.

Ellipticity fluctuation. In the previous section, we explicitly estimated two
moments of the initial state entropy density. In this section, we will use
these moments to derive an analytical estimate for the average ellipticity.
Additionally, we will validate this estimation, along with the moments stated
in Eqs. (27) and (28), using the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model.

The quantity $\epsilon_{2}^{2}\\{2\\}$ can be obtained by dividing Eq. (28) by
Eq. (27),

$\epsilon_{2}^{2}\\{2\\}=\langle\epsilon_{2}^{2}\rangle=\left\langle\frac{\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}}{\\{r^{2}\\}^{2}}\right\rangle\approx\frac{\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle}{\left\langle\\{r^{2}\\}\right\rangle^{2}}.$
(29)

To justify the last approximation equality, we utilize the Power distribution,
which can be expressed as
$p(\epsilon_{2})=(N-1)\epsilon_{2}(1-\epsilon_{2}^{2})^{(N-3)/2}$ [35]. The
Power distribution characterizes the distribution of $\epsilon_{2}$ for $N$
independent point-like sources that are distributed with a Gaussian profile.
It is determined that $\langle\epsilon_{2}^{2}\rangle=2/(N+1)$ using the Power
distribution. Furthermore, we can compute $\langle\\{r^{2}\\}^{2}\rangle$ and
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$ for
independent point-like sources distributed according to a Gaussian
distribution, using the methods outlined in the preceding section. We obtain
the same result for
$\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\rangle/\langle\\{r^{2}\\}^{2}\rangle$.
Discrepancies arise for finite-size sources. The comparison between
$\langle\epsilon_{2}^{2}\rangle$ from the Power distribution with finite-size
sources (not provided here) and
$\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\rangle/\langle\\{r^{2}\\}^{2}\rangle$
reveals a disparity on the order of $1/N^{3}$ which is neglected here. Hence,

$\begin{split}&\epsilon_{2}^{2}\\{2\\}=\frac{5}{7}\frac{1}{A\,\omega^{4}}\left[1+\alpha_{\epsilon}\left(\frac{a_{0}}{R_{0}}\right)^{2}\right.\\\
&\hskip
48.36958pt\left.+s_{\epsilon}\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{s_{\perp}}{\sigma^{\text{NN}}_{\text{inel}}}\right)+d_{\varphi}\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}\right.\\\
&\hskip
48.36958pt\left.+v_{\varphi}\frac{\sigma_{\text{fluc}}^{2}}{n_{c}}+b_{2,\epsilon}\beta_{2}^{2}+b_{3,\epsilon}\beta_{3}^{2}+\cdots\right].\end{split}$
(30)

To validate our estimate, we compare it with the realistic
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model in
the Au–Au collision with $\sigma_{\textbf{inel}}^{NN}=4.2\,$[fm2]. The results
are shown in Fig. 4. From the figure, we can see that our approximation
improves as the nucleon width increases. However, it still provides reasonable
agreement even for nucleon widths as small as $w=0.7\,$[fm]. The validation
related to the substructure of nucleons and the finite value of the cross-
section can be found in appendix F.

Based on the figure, we observe that the moment $\langle\\{r^{2}\\}\rangle$ is
not greatly impacted by the parameter $\sigma_{\text{fluc}}$, which aligns
with our calculations. This parameter primarily affects the fluctuation, while
the value of $\langle\\{r^{2}\\}\rangle$ is mainly determined by the nucleus
size. Additionally, the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ outcome is
not significantly influenced by $d_{\text{min}}$ because of the pre-generation
of the radial direction of the nucleons’ position, which prevents the
modification of the WS radius. In our case, we accomplished this by redefining
the radius. For each value of $d_{\text{min}}$, the radius of the WS
distribution is adjusted to maintain a fixed $\langle\\{r^{2}\\}\rangle$.
Furthermore, increasing $a_{0}$ enlarges the total nucleus size, subsequently
increasing $\langle\\{r^{2}\\}\rangle$.

---
Figure 5: Comparison of the isobar ratio between deformed nuclei Ru and Zr
from T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
(dots) with our analytical estimate (curves). The inaccuracy caused by the
small nucleon size is eliminated when calculating the ratio. However, the
impact of octupole deformation of Zr results in an inaccurate estimation of
the ratio, as shown in the middle panels.

The value of the quantity
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$ is
influenced only by the fluctuation in ultra-central Au-Au collisions. The
figure illustrates that as $\sigma_{\text{fluc}}$ increases, the value of the
quantity also increases. The size of the system also affects this quantity by
providing more room for fluctuation. Increasing $a_{0}$ leads to a larger
system, while an opposite trend is observed as $d_{\text{min}}$ increases.
This is because increasing $d_{\text{min}}$ reduces the available space for
fluctuation by removing the volume of the hard inner core of the nucleons from
the total volume in a nucleus.

We now discuss the behavior of $\epsilon_{2}\\{2\\}$ when we change the
initial state parameters. The ellipticity average can be thought of as a
normalized version of moment
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$.
Therefore, the dependence of $\epsilon_{2}\\{2\\}$ on the parameters is mostly
a competition between
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$ and
$\left\langle\\{r^{2}\\}\right\rangle$. In this case,
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$
dominates. Our calculation shows a smaller difference compared to the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ outcome
since some errors are canceled out in calculating the ratio. However,
inaccuracies due to small $w$ still affect $\epsilon_{2}\\{2\\}$ because
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$ has a
smaller dependence on $w$ compared to $\left\langle\\{r^{2}\\}\right\rangle$.
This dependence is reduced when we calculate the isobar ratio.

We have compared the dependence of $\epsilon_{2}\\{2\\}$ on $d_{\text{min}}$
between our numerical calculation of $d^{\prime}_{\varphi}$ and the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model. We
found that there is a very good agreement for $w=1.5\,$[fm]. We recall the
discussion about the implementation of short-range correlation for Monte Carlo
sampling in section III.2. We deduce that the pre-generation of radial
direction is effectively similar to implementing the effect of short-range
correlation into the WS distribution. This is because small values of
$d_{\mathrm{min}}$ do not modify the functionality of WS much. In three
dimensions, the issue for source placement is negligible, and the inaccuracy
of the iterative method is small compared to what we observe in the 1D toy
model shown in Fig. 9 (right) in appendix E. One may wonder if these issues
become more critical in other moments.

Figure 6:  Constrainting on $(\beta_{2}^{\text{Ru}})^{2}$ and
$(\beta_{2}^{\text{Zr}})^{2}$ for a given $R_{\epsilon}$, with different
values of $d_{\text{min}}$, $\sigma_{\text{fluc}}$, $\delta R$, and $\delta
a$. Here, the WS parameters of Zr nucleus is chosen
$R_{0}^{(\text{Zr})}=5.02\,$[fm] and $a_{0}^{(\text{Zr})}=0.52\,$[fm]. For Ru
nucleus, parameters are given as
$R_{0}^{(\text{Ru})}=R_{0}^{(\text{Zr})}+\delta R$ and
$a_{0}^{(\text{Ru})}=a_{0}^{(\text{Zr})}+\delta a$. The other parameters are
chosen as follows: $n_{c}=1$, $w=0.7\,$[fm], and
$\beta_{3}^{(\text{Zr})}=0.2$. The data point shows the location of
$\beta_{2}^{\text{Zr}}$ and $\beta_{2}^{\text{Ru}}$ inferred from electric
transition rate [36]. The results show that the most significant dependence
comes form $\sigma_{\text{fluc}}$ parameter and the least is the effect of
$R_{0}$ and $a_{0}$.

Isobar ratio. We also investigate our approximation for the isobar ratio of
deformed nuclei, where the WS parameters of two nuclei are different.
Specifically, we have considered the isobar nuclei Ru and Zr. Nucleus Ru, with
a radius of $R_{0}^{(\text{Ru})}$ and a skin thickness of
$a_{0}^{(\text{Ru})}$, is considered to have only $\beta_{2}^{(Ru)}$. On the
other hand, nucleus Zr, with $R_{0}^{(\text{Zr})}$ and $a_{0}^{(\text{Zr})}$,
exhibits non-vanishing $\beta_{2}^{(\text{Zr})}$ and $\beta_{3}^{(\text{Zr})}$
deformations. All other parameters of the two nuclei are assumed to be the
same. The ratio between eccentricities of two nuclei is shown by
$R_{\epsilon}$,

$R_{\epsilon}=\frac{\epsilon_{2}^{(\text{Ru})}\\{2\\}}{\epsilon_{2}^{(\text{Zr})}\\{2\\}}.$
(31)

We compare our analytical calculation with the results obtained from the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model. The
comparison is shown in Fig. 5. Our predictions of $\epsilon_{2}\\{2\\}$
underestimate the simulation for smaller values of nucleon width ($w$), but
the agreement improves as we move to larger values of $w$. However, the $w$
dependence has a small contribution to the ratio calculation, both in the
simulation and our analytical estimation. In the middle panel of Fig. 5, we
observe a systematic discrepancy between the simulation and the analytical
estimation. In the right panel, we ignored the effect of octupole deformation
of Zr and found excellent agreement. This observation suggests that our
estimation of the contribution of deformation in the appendix F.3 requires
improvement for nuclei that exhibit both quadrupole and octupole deformation
simultaneously.

##  V Constraining the nuclear structure via isobar ratios

It is known that in ultra-central heavy ion collisions, $\epsilon_{2}$ is
linearly correlated with the final elliptic flow, $v_{2}$ [25]. However, this
relationship depends on the specific observable being studied and requires
careful examination through hydrodynamic simulations. For example, in Ref.
[37], the isobar ratio in both the initial and final states was calculated
using the T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$
model for the initial state and the Trajectum hydrodynamic model for the
collective expansion. For the 0-5$\%$ centrality class, it was found that

$\frac{v_{2}^{(\text{Ru})}\\{2\\}}{v_{2}^{(\text{Zr})}\\{2\\}}=k_{\text{hyd}}\frac{\epsilon_{2}^{(\text{Ru})}\\{2\\}}{\epsilon_{2}^{(\text{Zr})}\\{2\\}},$
(32)

where $k_{\text{hyd}}\gtrsim 1$ is a constant. However, it should be noted
that the 0-5$\%$ centrality class is not completely populated by ultra-central
collisions. Therefore, it is important to examine the effect of hydrodynamic
response in ultra-central collisions to estimate a more accurate value for
$k_{\text{hyd}}$. In addition, we have found that our calculation at the
presence of octupole deformation of Zr shows a systematic difference compared
to simulation,

$\frac{\epsilon_{2}^{(\text{Ru})}\\{2\\}}{\epsilon_{2}^{(\text{Zr})}\\{2\\}}=k_{\text{corr}}\left.\frac{\epsilon_{2}^{(\text{Ru})}\\{2\\}}{\epsilon_{2}^{(\text{Zr})}\\{2\\}}\right|_{\text{Analytic}},$
(33)

where $k_{\text{corr}}\approx 0.96$, in our particular example. In Ref. [26],
it was reported that the value of
$v_{2}^{(\text{Ru})}\\{2\\}/v_{2}^{(\text{Zr})}\\{2\\}$ in the 0-5$\%$
centrality class is $1.0247\pm 0.0017$. Our calculation is more accurate for
events in a smaller width of the centrality class (less than 1$\%$), in which
a larger value for the ratio is expected [8]. Here, we assume
$k_{\text{corr}}k_{\text{hyd}}\approx 0.97$

In the ultra-central events, the primary contribution to $R_{\epsilon}$ arises
from the quadrupole deformation of Ru and Zr, denoted by
$(\beta_{2}^{\text{Ru}})^{2}-(\beta_{2}^{\text{Zr}})^{2}$, as referenced in
Ref. [7]. To assess the impact of various parameters in the UCSC (ultra-
central collisions), we modify the parameters while keeping the value of
$R_{\epsilon}$ constant. In Fig. 6, we demonstrate
$(\beta_{2}^{\text{Ru}})^{2}$ against $(\beta_{2}^{\text{Zr}})^{2}$ for a
fixed $R_{\epsilon}$. The values of the parameters $d_{\text{min}}$,
$\sigma_{\text{fluc}}$, $\delta R=R_{0}^{\text{Ru}}-R_{0}^{\text{Zr}}$, and
$\delta a=a_{0}^{\text{Ru}}-a_{0}^{\text{Zr}}$ are presented in the panels of
the figure. The following parameters are held constant in the plot:
$R_{0}^{(\text{Zr})}=5.02\,$[fm], $a_{0}^{(\text{Zr})}=0.52\,$[fm],
$\beta_{3}^{(\text{Zr})}=0.2$, $n_{c}=1$, and $w=0.7$[fm]. The data point in
the figure denotes the positions of $\beta_{2}^{\text{Zr}}$ and
$\beta_{2}^{\text{Ru}}$ inferred from electric transition rate [36]. It is
important to note that the inferred deformation values from transition rates
are not directly comparable with those obtained from heavy-ion collisions.
This depiction is presented to compare the model’s sensitivity with the
measurement accuracy. It is evident that the most influential parameters are
$\sigma_{\text{fluc}}$ and $d_{\text{min}}$. The effect of difference in the
radius is in the order of $d_{\text{min}}$ only for large values of $\delta
R$.

##  VI Summary

In this paper, we demonstrated that the functionality of the reduced thickness
function has minimal effect on the moments of initial state entropy density in
UCSC when the nucleon width is large enough, and the initial state has scale-
invariance property. We validated this argument by studying a few initial
state moments using the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ Monte
Carlo event generator (see Figs. 1 and 7).

In our study on UCSC, we utilized the cluster expansion method to determine
the impact of the one-body density and two-body correlation of the nucleus on
moments $\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$,
$\left\langle\\{r^{2}\\}^{2}\right\rangle$, and $\epsilon_{2}\\{2\\}$. To
carry out these calculations, we employed the WS distribution for the
distribution of nucleons and hard inner core for nucleons to enforce
correlation among them. Our findings indicate that even though the minimum
distance between nucleons modifies the one-body density, it is still possible
to describe this density with WS, with a larger radius up to at least
$d_{\text{min}}=0.4\,$[fm]. We also compared different methods of Monte Carlo
sampling of sources in the presence of short-range correlations and
demonstrated that our approach can be used to assess the accuracy of short-
range correlations sampling. This revealed that incorporating information on
short-range correlations in the one-body density is essential for obtaining
accurate results. We calculated the effect of cross-section and nucleon
substructure. To validate, we compared results with
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
outcome.

We investigated the impact of initial state parameters on the elliptic flow
ratio for Ru–Ru and Zr–Zr in ultra-central collisions. Our findings confirmed
that the quadrupole deformation of Ru and Zr has the most significant
contribution. We also found that the effect of constituent weight fluctuation
and two-body correlation must be taken into account when inferring the nuclei
deformation.

Although we compared our calculation with the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model in
this study, our result should be valid for any other initial state models that
satisfy the scale-invariance.

## Acknowledgment

The authors are thankful to Jiangyong Jia for discussion and comment, we also
thank Matthew Luzum, and Jean-Yves Ollitrault for their helpful discussions.
S.F.T. is supported by the Deutsche Forschungsgemeinschaft (DFG) through the
grant TA 1975/1-1.

APPENDICES

##  Appendix A Convergence of triangularity, symmetric cumulants and symmetry
plane correlation

We examine three additional initial state moments (cumulants):
$\epsilon_{3}\\{2\\}$,
$\langle\epsilon_{2}^{2}\epsilon_{3}^{2}\rangle-\langle\epsilon_{2}^{2}\rangle\langle\epsilon_{3}^{2}\rangle$,
and
$\langle\epsilon_{2}^{2}\epsilon_{4}\cos\left(4\varphi_{2}-4\varphi_{4}\right)\rangle$.
These moments exhibit the same convergence behavior as discussed in Section II
and are depicted in Figure 7. The symmetric cumulant
$\langle\epsilon_{2}^{2}\epsilon_{3}^{2}\rangle-\langle\epsilon_{2}^{2}\rangle\langle\epsilon_{3}^{2}\rangle$
and the symmetric plane correlation
$\langle\epsilon_{2}^{2}\epsilon_{4}\cos\left(4\varphi_{2}-4\varphi_{4}\right)\rangle$
show significantly larger values at non-central and asymmetric collisions. In
UCSC, convergence occurs at the same order of $w$ as the other moments
discussed in Section II.

Figure 7: More moments similar to those shown in Fig. 1. For symmetric
cumulant, convergence seems to occur even in non-central and asymmetric
collisions.

##  Appendix B Effect of nucleon size and impact parameter on
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ reduced
thickness function

Figure 8: The isolated participants in
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ reduced
thickness function for $p=-1,0,1$.

Isolated participants. We start with analyzing an idealized scenario. In Fig.
8, we depict one isolated participant from each target and projectile. Dashed
blue and black circles show the location of these participants in the
transverse space, and the radius of the circles indicates the width $w$. Here,
we assume two nucleons are separated by a distance $d$ along the $x$ axis.
Assuming these nucleons are isolated from the rest of the target and
projectile thickness functions (visualized by the shaded region on the upper
left part of each panel), we can decompose $T_{A(B)}(\mathbf{x})$ as

$T_{A(B)}(\mathbf{x})=\overline{T}_{A(B)}(\mathbf{x})+T_{A(B)}^{\text{isol}}(\mathbf{x}),$
(34)

such that

$\overline{T}_{A(B)}(\mathbf{x})\,T_{A(B)}^{\text{isol}}(\mathbf{x})\approx
0,$ (35)

for all points in the transverse plane. This assumption leads to a
decomposition in the reduced thickness function as well,

$T_{R}(\mathbf{x})\approx\overline{T}_{R}(\mathbf{x})+T_{R}^{\text{isol}}(\mathbf{x}).$
(36)

Focusing on the isolated target and projectile participants, located at
$x_{a}=-x_{b}=d/2$, $y_{a}=y_{b}=0$, one finds the reduced thickness function
for $p=1$,

$\displaystyle T^{\text{isol}}_{R}(1;\mathbf{x})=$
$\displaystyle\frac{\mathcal{N}}{2}\left(\gamma_{t}\frac{e^{-\frac{(x-d/2)^{2}+y^{2}}{2w^{2}}}}{2\pi
w^{2}}+\gamma_{p}\frac{e^{-\frac{(x+d/2)^{2}+y^{2}}{2w^{2}}}}{2\pi
w^{2}}\right),$ (37a)

and for $p=0$,

$\displaystyle
T^{\text{isol}}_{R}(0;\mathbf{x})=\left(\mathcal{N}\sqrt{\gamma_{t}\gamma_{p}}e^{-\frac{d^{2}}{8w^{2}}}\right)\frac{e^{-\frac{x^{2}+y^{2}}{2w^{2}}}}{2\pi
w^{2}}.$ (38a)

In Fig. 8, for two distant nucleons in $p=1$,
$T^{\text{isol}}_{R}(1;\mathbf{x})$ is the summation of two Gaussian elongated
along the connecting line between the center of two nucleons (Fig. 8 right
panel). The total entropy in this case, is a constant given by

$dS^{\text{\,isol}}/d\eta=\mathcal{N}\left(\frac{\gamma_{a}+\gamma_{b}}{2}\right).$
(39)

For $p=0$, the distribution is always centered at the middle of two nucleons
with a round Gaussian distribution with width $w$ even in case
$\gamma_{t}\neq\gamma_{p}$ as it is shown in Fig. 8 middle panel. The total
entropy reads as

$dS^{\text{\,isol}}/d\eta=\mathcal{N}\sqrt{\gamma_{a}\gamma_{b}}\,\exp\left[-\frac{d^{2}}{8w^{2}}\right].$
(40)

The exponential term leads to the fact that the more distant the two sources
are, the less contribution the isolated part delivers to the total entropy. A
similar analysis can be done for $p=-1$. The distribution is elongated
perpendicular to the connecting line as it is seen from Fig. 8 left panel. For
small values of $w$, the chance of finding isolated participants increases.

Moments of $T_{R}(1;\mathbf{x})$ for arbitrary value of $w$.  In more
realistic scenarios, finding isolated sources similar to the previous part is
unlikely unless we have $w\ll\bar{d}$ where $\bar{d}^{2}$ stands for the
inverse source density in the transverse space,

$\bar{d}^{2}\sim\pi R_{0}^{2}/N_{A(B)}.$ (41)

Assume a typical event with $N_{A}$ and $N_{B}$ number of sources at
$\\{\mathbf{x}_{a,i}\\}$ and $\\{\mathbf{x}_{b,i}\\}$ in the transverse space.
For $p=1$, at each point,
$\mathbf{x}_{i}\in\\{\mathbf{x}_{a,i}\\}\cup\\{\mathbf{x}_{b,i}\\}$, a
Gaussian distribution is located. Hence,

$T_{R}(1;\mathbf{x})=\frac{1}{S(1)}\frac{1}{2}\sum_{i=1}^{N_{A}+N_{B}}\frac{1}{2\pi
w^{2}}\exp\left[-\frac{|\mathbf{x}-\mathbf{x}_{i}|^{2}}{2w^{2}}\right],$ (42)

where we ignore the effect of normalization fluctuation $\gamma_{i}$ from now
on. Here, the total entropy is given by

$dS/d\eta=\frac{1}{2}\sum_{i=1}^{N_{A}+N_{B}}1=\frac{N_{A}+N_{B}}{2}.$ (43)

The characteristic function can be calculated from Eq. (42),

$\begin{split}\\{e^{i\bm{k}\cdot\mathbf{x}}\\}&=\frac{1}{2}\sum_{i=1}^{N_{A}+N_{B}}\frac{1}{2\pi
w^{2}}\exp\left[i\bm{k}\cdot\mathbf{x}-\frac{|\mathbf{x}-\mathbf{x}_{i}|^{2}}{2w^{2}}\right]\\\
&=e^{-\frac{k^{2}w^{2}}{2}}\frac{1}{2}\sum_{i=1}^{N_{A}+N_{B}}\,e^{i\bm{k}\cdot\mathbf{x}_{i}}\\\
&=e^{-\frac{k^{2}w^{2}}{2}}\left.\\{e^{i\bm{k}\cdot\mathbf{x}}\\}\right|_{w\to
0}.\end{split}$ (44)

From above, we find that

$\displaystyle\\{x^{2}+y^{2}\\}=2w^{2}+\left.\\{x^{2}+y^{2}\\}\right|_{w=0},$
(45a) $\displaystyle\\{y^{2}-x^{2}\\}=\left.\\{x^{2}+y^{2}\\}\right|_{w=0},$
(45b) $\displaystyle\\{2xy\\}=\left.\\{2xy\\}\right|_{w=0}.$ (45c)

This means that for $p=1$, we can assume the nucleons are point-like sources,
and the effect of nucleon width can be included trivially. The average of the
point-like source second-order moments are given by

$\displaystyle\left.\\{x^{2}+y^{2}\\}\right|_{w=0}=$ (46a)
$\displaystyle\hskip
22.76228pt\frac{N_{A}}{N_{A}+N_{B}}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{x_{b}^{2}+y_{b}^{2}\\}_{s},$
$\displaystyle\left.\\{y^{2}-x^{2}\\}\right|_{w=0}=$ (46b)
$\displaystyle\hskip
22.76228pt\frac{N_{A}}{N_{A}+N_{B}}\\{y_{a}^{2}-x_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{y_{b}^{2}-x_{b}^{2}\\}_{s},$
$\displaystyle\left.\\{2xy\\}\right|_{w=0}=$ (46c) $\displaystyle\hskip
22.76228pt\frac{N_{A}}{N_{A}+N_{B}}\\{2x_{a}y_{a}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{2x_{b}y_{b}\\}_{s},$

where we have used the following notation:

$\\{\cdots\\}_{s}=\frac{1}{N_{A(B)}}\sum_{i=1}^{N_{A(B)}}\cdots.$ (47)

| $w\ll\bar{d}$ | $w\gg\bar{d}$
---|---|---
$p=1$ | $S$ | $(N_{A}+N_{B})/2$ | $(N_{A}+N_{B})/2$
| $\\{x^{2}+y^{2}\\}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{x_{b}^{2}+y_{b}^{2}\\}_{s}+2w^{2}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{x_{b}^{2}+y_{b}^{2}\\}_{s}+2w^{2}$
| $\\{y^{2}-x^{2}\\}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{y_{a}^{2}-x_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{y_{b}^{2}-x_{b}^{2}\\}_{s}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{y_{a}^{2}-x_{a}^{2}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{y_{b}^{2}-x_{b}^{2}\\}_{s}$
| $\\{2xy\\}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{2x_{a}y_{a}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{2x_{b}y_{b}\\}_{s}$ | $\frac{N_{A}}{N_{A}+N_{B}}\\{2x_{a}y_{a}\\}_{s}+\frac{N_{B}}{N_{A}+N_{B}}\\{2x_{b}y_{b}\\}_{s}$
$p=0$ | $S$ | $N_{A}\,N_{B}\\{\exp\left[-d_{ij}^{2}/8w^{2}\right]\\}_{ij,s}$ | $\sqrt{N_{A}\,N_{B}}$
| $\\{x^{2}+y^{2}\\}$ | $\left\\{(x_{ij}^{2}+y_{ij}^{2})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}/\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}+2w^{2}$ | $\frac{1}{2}\Big{(}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\\{x_{b}^{2}+y_{b}^{2}\\}_{s}+2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}\Big{)}+2w^{2}$
| $\\{y^{2}-x^{2}\\}$ | $\left\\{(y_{ij}^{2}-x_{ij}^{2})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}/\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}$ | $\frac{1}{2}\Big{(}\\{y_{a}^{2}-x_{a}^{2}\\}_{s}+\\{y_{b}^{2}-x_{b}^{2}\\}_{s}+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}-2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}\Big{)}$
| $\\{2xy\\}$ | $\left\\{(2x_{ij}y_{ij})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}/\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}$ | $\\{x_{a}y_{a}\\}_{s}+\\{x_{b}y_{b}\\}_{s}+\\{x_{a}\\}_{s}\\{y_{b}\\}_{s}+\\{y_{a}\\}_{s}\\{x_{b}\\}_{s}$
Table 2: Total entropy and second order moments of reduced thickness function
for $p=0$ and $p=1$ in large and small $w$.

Limit of $T_{R}(0;\mathbf{x})$ for small and large $w$.  For $p=0$, we start
with rewriting the reduced thickness function as follows,

$\begin{split}&T_{R}(0;\mathbf{x})=\\\
&\frac{1}{S}\left[\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\left(\frac{\lambda_{ij}}{2\pi
w^{2}}\exp\left[-\frac{|\mathbf{x}-\mathbf{x}_{ij}|^{2}}{2w^{2}}\right]\right)^{2}\right]^{1/2},\end{split}$
(48)

where $\lambda_{ij}(w)=\exp\left[-\frac{d_{ij}^{2}}{8w^{2}}\right]$, and

$\mathbf{x}_{ij}=\frac{\mathbf{x}_{a,i}+\mathbf{x}_{b,j}}{2},\qquad
d_{ij}=|\mathbf{x}_{a,i}-\mathbf{x}_{b,j}|.$ (49)

This way of writing is inspired by Eq. (38a), which indicates that at $p=0$,
there is effectively a source at the center of a target and projectile sources
weighted exponentially by their distance. However, due to the power of $1/2$,
we cannot find a full analytical result for the characteristic function. For
that reason, we find an analytical estimation only in two limits, small and
large $w$.

In the $w\ll\bar{d}$ limit, as we argued for a single pair of isolated
participants, we employ the fact that

$(\alpha^{2}_{1}(\mathbf{x})+\alpha^{2}_{2}(\mathbf{x})+\cdots+\alpha^{2}_{n}(\mathbf{x}))^{1/2}\approx\alpha_{1}(\mathbf{x})+\cdots+\alpha_{n}(\mathbf{x})$
(50)

if $\alpha_{i}(\mathbf{x})\alpha_{j}(\mathbf{x})\approx 0$ for all $i\neq j$
pairs at any given $\mathbf{x}$. The sources are well-separated in the limit
$w\ll\bar{d}$ and fulfill this condition approximately. Using it, one can
simplify Eq. (48) as

$\begin{split}&\left.T_{R}(0;\mathbf{x})\right|_{w\ll\bar{d}}=\\\ &\hskip
28.45274pt\frac{1}{S}\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\frac{\lambda_{ij}(w)}{2\pi
w^{2}}\exp\left[-\frac{|\mathbf{x}-\mathbf{x}_{ij}|^{2}}{2w^{2}}\right].\end{split}$
(51)

The thickness function in this limit is reduced to that in case $p=1$, but
this time, we have $N_{A}N_{B}$ number of sources centered in the middle of
each pair of target and projectile participants shown by $\mathbf{x}_{ij}$.
The contribution of sources is weighted based on the distance $d_{ij}$ of the
participant sources. The total entropy in $w\ll\bar{d}$ is given by

$\begin{split}dS/d\eta&=\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\exp\left[-\frac{d_{ij}^{2}}{8w^{2}}\right]\\\
&=N_{A}N_{B}\\{\exp\left[-d_{ij}^{2}/8w^{2}\right]\\}_{ij,s}.\end{split}$ (52)

and second order moments are given by

$\displaystyle\\{x^{2}+y^{2}\\}=\frac{\left\\{(x_{ij}^{2}+y_{ij}^{2})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}}{\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}}+2w^{2},$
(53a)
$\displaystyle\\{y^{2}-x^{2}\\}=\frac{\left\\{(y_{ij}^{2}-x_{ij}^{2})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}}{\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}},$
(53b)
$\displaystyle\\{2xy\\}=\left\\{(2x_{ij}y_{ij})e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s}/\left\\{e^{-\frac{d_{ij}^{2}}{8w^{2}}}\right\\}_{ij,s},$
(53c)

where

$\\{\cdots\\}_{ij,s}=\frac{1}{N_{A}N_{B}}\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\cdots.$
(54)

We can estimate the asymptotic behavior of the total entropy in the
$w\gg\bar{d}$ limit. Considering
$\alpha_{i}(\mathbf{x})=\frac{\lambda_{ij}}{2\pi
w^{2}}e^{-\frac{|\mathbf{x}-\mathbf{x}_{ij}|^{2}}{2w^{2}}}$ in Eq. (48), one
assumes the total distribution inside the square root can be approximated with
the square of a single radially Gaussian distribution,

$\alpha^{2}_{1}(\mathbf{x})+\alpha^{2}_{2}(\mathbf{x})+\cdots+\alpha^{2}_{n}(\mathbf{x})\approx\left[\frac{\Lambda}{|\Sigma^{2}|}\exp\left[-\frac{|\mathbf{x}|^{2}}{2\Sigma^{2}(\varphi)}\right]\right]^{2},$
(55)

where $\varphi=\arctan(y/x)$, and we have assumed the average location of
sources is at the origin. Here, the normalization $|\Sigma^{2}|$ is obtained
via

$|\Sigma^{2}|:=\int_{0}^{2\pi}d\varphi\,\Sigma^{2}(\varphi).$ (56)

In the present study, only the ellipticity of the initial state is studied. In
this case, we approximate the $\varphi$ dependent width $\Sigma^{2}(\varphi)$
as follows,

$\frac{1}{\Sigma^{2}(\varphi)}=\bm{\varphi}^{T}\cdot\bm{\Sigma}^{-1}\cdot\bm{\varphi},$
(57)

where

$\bm{\varphi}=\begin{pmatrix}\cos\varphi\\\
\sin\varphi\end{pmatrix},\qquad\bm{\Sigma}=\begin{pmatrix}\Sigma_{xx}^{2}&\Sigma_{xy}^{2}\\\
\Sigma_{xy}^{2}&\Sigma_{yy}^{2}\end{pmatrix}.$ (58)

This leads to the square of a 2D-Gaussian for the right-hand side of Eq. (55).
By using the above estimation, we eventually find

$\begin{split}&T_{R}(0;\mathbf{x})\approx\\\
&\frac{1}{S}\left(\frac{\Lambda}{2\pi\sqrt{\det\bm{\Sigma}}}\right)\exp\left[-\frac{1}{2}\left(\mathbf{x}^{T}\cdot\bm{\Sigma}^{-1}\cdot\mathbf{x}\right)\right],\end{split}$
(59)

where

$\displaystyle\Lambda^{2}$ $\displaystyle=\sum_{ij}\lambda_{ij}^{2},$ (60a)
$\displaystyle\Sigma_{xx}^{2}$
$\displaystyle=w^{2}+\frac{2\sum_{ij}x_{ij}^{2}\lambda_{ij}^{2}}{\sum_{ij}\lambda_{ij}^{2}},$
(60b) $\displaystyle\Sigma_{yy}^{2}$
$\displaystyle=w^{2}+\frac{2\sum_{ij}y_{ij}^{2}\lambda_{ij}^{2}}{\sum_{ij}\lambda_{ij}^{2}},$
(60c) $\displaystyle\Sigma_{xy}^{2}$
$\displaystyle=\frac{2\sum_{ij}x_{ij}y_{ij}\lambda_{ij}^{2}}{\sum_{ij}\lambda_{ij}^{2}}.$
(60d)

In the above, we have used $\sum_{ij}:=\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}$
as a shorthand notation.

Noting the fact that $\lambda_{ij}\to 1$ for larege $w$, the total entropy is
given by

$dS/d\eta\approx\left[\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}1\right]^{1/2}=\sqrt{N_{A}N_{B}}.$
(61)

From Eqs. (47), (49), and (54), one finds,

$\begin{split}\left\\{x_{ij}^{2}\right\\}_{ij,s}&=\frac{1}{4}\frac{1}{N_{A}N_{B}}\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}(x_{i,a}+x_{j,b})^{2}\\\
&=\frac{1}{4}\frac{1}{N_{A}N_{B}}\sum_{i=1}^{N_{A}}\sum_{j=1}^{N_{B}}\left(x_{i,a}^{2}+x_{j,b}^{2}+2x_{i,a}\,x_{j,b}\right)\\\
&=\frac{1}{4}\left(\\{x_{a}^{2}\\}_{s}+\\{x_{b}^{2}\\}_{s}+2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}\right).\end{split}$
(62)

Similarly,

$\displaystyle\left\\{y_{ij}^{2}\right\\}_{ij,s}=\frac{1}{4}\left(\\{y_{a}^{2}\\}_{s}+\\{y_{b}^{2}\\}_{s}\right.$
(63a) $\displaystyle\hskip
76.82234pt\left.+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}\right),$
$\displaystyle\left\\{x_{ij}y_{ij}\right\\}_{ij,s}=\frac{1}{4}\left(\\{x_{a}y_{a}\\}_{s}+\\{x_{b}y_{b}\\}_{s}\right.$
(63b) $\displaystyle\hskip
76.82234pt\left.+\\{x_{a}\\}_{s}\\{y_{b}\\}_{s}+\\{x_{b}\\}_{s}\\{y_{a}\\}_{s}\right).$

Using above, second order moments of $T_{R}(0;\mathbf{x})$ in large $w$ limit
are given as

$\displaystyle\\{x^{2}+y^{2}\\}=\frac{1}{2}\Big{(}\\{x_{a}^{2}+y_{a}^{2}\\}_{s}+\\{x_{b}^{2}+y_{b}^{2}\\}_{s}$
(64a) $\displaystyle\hskip
42.67912pt+2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}\Big{)}+2w^{2},$
$\displaystyle\\{y^{2}-x^{2}\\}=\frac{1}{2}\Big{(}\\{y_{a}^{2}-x_{a}^{2}\\}_{s}+\\{y_{b}^{2}-x_{b}^{2}\\}_{s}$
(64b) $\displaystyle\hskip
42.67912pt+2\\{y_{a}\\}_{s}\\{y_{b}\\}_{s}-2\\{x_{a}\\}_{s}\\{x_{b}\\}_{s}\Big{)},$
$\displaystyle\\{2xy\\}=\\{x_{a}y_{a}\\}_{s}+\\{x_{b}y_{b}\\}_{s}+\\{x_{a}\\}_{s}\\{y_{b}\\}_{s}$
(64c) $\displaystyle\hskip 25.6073pt+\\{y_{a}\\}_{s}\\{x_{b}\\}_{s}.$

We have summarized the result in Table. 2. Assuming
$\\{x_{a}+x_{b}\\}_{s}=\\{y_{a}+y_{b}\\}_{s}=0$, we have
$\langle\\{x_{a}\\}_{s}\rangle\approx-\langle\\{x_{b}\\}_{s}\rangle=\tilde{b}/2$
and $\langle\\{y_{a}\\}_{s}\rangle\approx-\langle\\{y_{b}\\}_{s}\rangle=0$ for
non-central collisions where $\tilde{b}$ should be proportional to the impact
parameter of the event, $b$. For the central collisions of same nuclei, we
have $\tilde{b}=0$ and $N_{A}=N_{B}$, meaning the total entropy and second
order moments are the same for $p=0$ and $p=1$.

##  Appendix C Event-by-event fluctuation and cluster-expansion

In this section, we discuss the details of calculations based on the cluster
expansion.

###  C.1 Decomposition of constituents contribution and source position

We start with the collision of two nuclei, each with mass number A. Assuming
nucleons as structureless spherical objects, the reduced thickness function is
given by

$T_{R}(\mathbf{x})\approx\frac{\mathcal{N}}{2}\sum_{i=1}^{N_{\text{part}}}\gamma_{i}\int
dz\,\rho_{\text{nucleon}}(\vec{x}-\vec{x}_{i}),$ (65)

where we have used Eq. (1) with $p=1$, which is similar to any initial state
with the scale-invariance property at large nucleon size. Here, the location
of the nucleon is represented by $\vec{x}_{i}$, and the number of participants
is denoted by $N_{\text{part}}$. Following Refs.[14, 20], the parameter
$\gamma_{i}$ is randomly selected from a gamma distribution that has a mean
value of unity and a variance of $\sigma_{\text{fluc}}^{2}$. The shape of
$T_{R}(\mathbf{x})$ can be analyzed by calculating its characteristic
function,

$e^{W_{R}(\bm{k})}=\frac{\int
d\mathbf{x}e^{i\bm{k}\cdot\mathbf{x}}T_{R}(\mathbf{x})}{\int
d\mathbf{x}T_{R}(\mathbf{x})}.$ (66)

The expansion coefficient of $e^{W_{R}(\bm{k})}$ ($W_{R}(\bm{k})$) in terms of
$\bm{k}$ leads to the moments (cumulants) of $T_{R}(\mathbf{x})$ [38]. We
label the averages with superscript ${(c)}$ if we refer to cumulants.
Substituting $T_{R}$ from Eq. (65) into the above, we find that the
characteristic function can be written as
$e^{W_{\text{nucleon}}(\bm{k})}e^{W_{\text{source}}(\bm{k})}$. For Gaussian
nucleon with width $w$, one finds
$W_{\text{nucleon}}(\bm{k})=-|\bm{k}|^{2}w^{2}/2$. Therefore, the cumulants of
reduced thickness function depends on the nucleons width trivially only in
$\\{r^{2}\\}$. The non-trivial information, therefore, encoded in the position
of the sources,

$e^{W_{\text{source}}(\bm{k})}=\frac{\sum_{i=1}^{N_{\text{part}}}\gamma_{i}e^{i\bm{k}\cdot\mathbf{x}_{i}}}{\sum_{i=1}^{N_{\text{part}}}\gamma_{i}}.$
(67)

###  C.2 Source fluctuation distribution and moments average

Decomposing the sources into target and projectile participants. The initial
state event-by-event fluctuation comes from the fluctuation of the sources
(nucleons, constituent). Suppose we have a distribution of the location of
sources given by
$P_{s}(\gamma_{1},\vec{x}_{1},\ldots,\gamma_{N_{\text{tot}}},\vec{x}_{N_{\text{tot}}})$,
where $N_{\text{tot}}=n_{c}N_{\text{part}}$, and assume large
$\sigma_{\text{inel}}^{\text{NN}}$ such that all nucleons participate in the
collision. In this case, no matter how complex the distribution is, the scale-
invariance in UCSC allows us to break it down into two separate distributions
of target and projectile.

$\displaystyle
P_{s}(\gamma,\vec{x})=P_{A}(\gamma,\vec{x})P_{B}(\gamma,\vec{x}),$ (68)

where $P_{A(B)}$ is the distribution of participant sources from the target
(projectile). The correction to this decomposition is discussed in appendix.
F.2. The parameters $\gamma$ and $\vec{x}$ are used as collective coordinates
that run through the number of sources in the target (projectile) in
$P_{A(B)}$, or their combination in $P_{s}$. The average over position of
sources is given as follows (see Eq. (67)),

$\\{f\\}_{s}=\frac{\sum_{i=1}^{N}\gamma_{i}f(\mathbf{x}_{i})}{\sum_{i=1}^{N}\gamma_{i}},$
(69)

where $f(\mathbf{x}_{i})$ is a given function, defining our moment. The effect
of nucleon width appears only in the cumulants $\\{r^{2}\\}^{(c)}$. In the
present study, we are particularly interested in averages $\\{r^{2}\\}^{(c)}$
and $\\{r^{2}e^{i2\varphi}\\}^{(c)}$ where due to the rotational invariance of
WS at the leading order, these cumulants are equivalent with the moments
$\\{r^{2}\\}$ and $\\{r^{2}e^{i2\varphi}\\}$.

Assuming that the $\gamma$ variables fluctuate independently of the source
location, we can calculate the average of the moments using

$\langle\\{f\\}_{s}\rangle=\int
d^{N}\gamma\,d^{2N}x\,d^{N}z\,\\{f\\}_{s}\,P_{s}(\gamma,\vec{x}).$ (70)

In particular, one finds the following relation between moments and the
nucleons position averages:

$\displaystyle\langle\\{r^{2}\\}\rangle=\langle\\{r^{2}\\}_{s}\rangle+2w^{2},$
(71a)
$\displaystyle\langle\\{r^{2}\\}^{2}\rangle=\left\langle\left(\\{r^{2}\\}_{s}+2w^{2}\right)^{2}\right\rangle,$
(71b)
$\displaystyle\langle\left|\\{r^{2}e^{i2\varphi}\\}\right|^{2}\rangle=\langle\left|\\{r^{2}e^{i2\varphi}\\}_{s}\right|^{2}\rangle.$
(71c)

For nucleons with substructure, $w$ should be replaced by $v$. For symmetric
collisions $A=B$, one can assume
$\sum_{i=1}^{N_{A}}\gamma_{i}\approx\sum_{i=1+N_{A}}^{N_{A}+N_{B}}\gamma_{i}$,
and decompose Eq. (69) into two parts,
$\\{f\\}_{s}\approx\frac{1}{2}\left[\\{f\\}_{s,A}+\\{f\\}_{s,B}\right]$. Since
$P_{A}$ and $P_{B}$ are same distributions in Eq. (68), one immediately deduce
that

$\displaystyle\langle\\{f\\}_{s}\rangle=\langle\\{f\\}_{s,A}\rangle,$ (72a)
$\displaystyle\langle\\{f_{1}\\}_{s}\\{f_{2}\\}_{s}\rangle=\frac{1}{2}\left[\langle\\{f_{1}\\}_{s,A}\\{f_{2}\\}_{s,A}\rangle+\langle\\{f_{1}\\}_{s,A}\rangle\langle\\{f_{2}\\}_{s,A}\rangle\right],$
(72b)

The average on the left-hand side is taken over $P_{s}$, while on the right-
hand side, it is taken over $P_{A}$, given $\mathcal{O}_{A}$. Therefore, we
only need to calculate observables for one nucleus participant. The correction
due to the finite nucleon-nucleon inelastic cross-section will be discussed in
appendix F.2.

The effect of the constituent weight fluctuation on observables. As mentioned,
the constituent weight fluctuation is followed by the gamma distribution in
Eq. (23). Given that fluctuations are assumed independent, referring to Eq.
(69), the the averages such as $\langle\\{f\\}_{s}\rangle$ or
$\langle\\{f_{1}\\}_{s}\\{f_{2}\\}_{s}\rangle$ can be written as the average
of numerator divided by denominator. Using $p_{\text{fluc}}(\gamma)$ to find
the denominator average, one estimates that

$\displaystyle\langle\\{f\\}_{s}\rangle$
$\displaystyle\approx\left\langle\frac{1}{N}\sum_{i=1}^{N}\gamma_{i}f(\mathbf{x}_{i})\right\rangle,$
(73a) $\displaystyle\langle\\{f_{1}\\}_{s}\\{f_{2}\\}_{s}\rangle$
$\displaystyle\approx\left\langle\frac{1}{N^{2}}\sum_{i,j=1}^{N}\gamma_{i}\gamma_{j}f_{1}(\mathbf{x}_{i})f_{2}(\mathbf{x}_{j})\right\rangle.$
(73b)

There is an extra term proportional to $\sigma_{\text{fluc}}^{2}/N^{3}$ in Eq.
(73b). This is ignored in this study as it is a subleading contribution in the
$1/N$ expansion.

###  C.3 Initial state fluctuation from cluster expansion method

For independent nucleons without structure, the distribution of A sources is
given by:

$P_{A}(\gamma,\vec{x})=\prod_{i=1}^{A}p_{\text{fluc}}(\gamma_{i})p_{A}(\vec{x}_{i}),$
(74)

which allows us to decompose the average on the right-hand side of equation
(73a) and perform the integration over a single source distribution. This
leads to Eq. (24). In the case of calculating second-order moments average,
$\langle\\{f_{1}\\}_{s,A}\\{f_{2}\\}_{s,A}\rangle$, we find Eq. (25) from
equation (73b) where we have used $\int
d\gamma\,\gamma^{2}\;p_{\text{fluc}}(\gamma)=1+\sigma_{\text{fluc}}^{2}$. To
include the effect of short-range correlation, we employ the statistical
mechanics analogy as discussed in Sec. III where the distribution is defined
as $P_{A}(\vec{x})=F_{A}(\vec{x})/\mathcal{Z}_{A}$.

Although we usually assume a fixed number of sources, when studying
fluctuation in cluster expansion, closed-form results can be obtained in
grand-canonical systems. It should be noted that the calculated observables
from canonical and grand-canonical systems may differ with terms in the order
of $O(1/A)$. With this caveat in mind, we can define grand-canonical partition
functions as follows:

$\mathscr{Q}(z,V,\lambda_{a})=\sum_{n=0}^{\infty}z^{n}\mathcal{Q}_{n}(V,\lambda_{a}),$
(75)

where the canonical partition function is given as
$\mathcal{Q}_{A}(V,\lambda_{a})=\mathcal{Z}_{A}(V,\lambda_{a})/h^{3A}A!$. The
dimensionful parameter $h$ is included to make the partition functions
dimensionless. The Mayer function $f_{ij}$ is defined as
$c_{2}(\vec{x}_{i},\vec{x}_{j})=1+f_{ij}$, which leads to

$\mathcal{Q}_{A}(V,\lambda_{a})=\frac{1}{h^{3A}A!}\int\mathcal{D}^{3A}x\prod_{1\leq
i<j\leq A}(1+f_{ij}),$ (76)

where the measure is defined as

$\mathcal{D}^{3}x(\lambda_{a})=d^{3}x\;p_{A}(\vec{x})\;e^{\frac{1}{A}\sum_{a}\lambda_{a}f_{a}(\vec{x}_{i})}.$
(77)

In the limit of no correlation, we have $f_{ij}=0$, and

$\mathcal{Q}_{A}(V,\lambda_{a})=\mathcal{Q}_{1}^{A}(V,\lambda_{a})/A!=\left[\int\mathcal{D}^{3}x(\lambda_{a})\right]^{A}/h^{3A}A!.$
(78)

The averages $\langle\\{f_{1}\\}\rangle$ and
$\langle\\{f_{1}\\}\\{f_{2}\\}\rangle$ corresponds to $f_{1}(\mathbf{x})$ and
$f_{2}(\mathbf{x})$ in Eq. (24) and (25) can be obtained by appropriately
differentiating $\mathcal{Q}_{A}(V,\lambda_{a})$ with respect to $\lambda_{a}$
and setting $\lambda_{a}$ to zero in the end.

To calculate the canonical partition function in the presence of correlation,
one needs to evaluate integrals of the form
$\int\mathcal{D}^{3}x_{1}\cdots\mathcal{D}^{3}x_{N}f_{\alpha_{1}}f_{\alpha_{2}}\cdots
f_{\alpha_{k}}$, where $\alpha_{i}$ represents a pair of distinct numbers
between 1 and $A$. Depending on the values of the pairs $\alpha_{i}$, the
integral may or may not be decomposed into lower-dimensional integrals or
clusters. In the presence of correlations, finding the grand-canonical
partition function is an easier task, and it is given by:

$\log\mathscr{Q}(z,\lambda_{a})=\sum_{\ell=1}b_{\ell}(\lambda_{a})z^{\ell}.$
(79)

The first term in the above expansion corresponds to the grand-canonical
partition function between sources with no correlations. The coefficients
$b_{\ell}(\lambda_{a})$ are Mayer cluster integrals. The 1-cluster and
2-cluster integrals read as:

$\displaystyle b_{1}(\lambda_{a})$
$\displaystyle=\int\mathcal{D}^{3}x(\lambda_{a}),$ (80a) $\displaystyle
b_{2}(\lambda_{a})$
$\displaystyle=\frac{1}{2}\int\mathcal{D}^{3}x_{1}(\lambda_{a})\mathcal{D}^{3}x_{2}(\lambda_{a})\;f_{12}.$
(80b)

Now, we are able to connect the cluster integrals to the desired moments. The
average value of $A$ can be obtained as follows,

$\langle A\rangle=\left.z\frac{\partial\log\mathscr{Q}}{\partial
z}\right|_{\lambda_{a}\to 0}=\sum_{\ell=1}^{\infty}\ell b_{\ell}(0)z^{\ell}.$
(81)

In the absence of correlation, only 1-cluster integral is non-vanishing,
meaning $\langle A\rangle$ is equal to $b_{1}(0)z$. From the definition in Eq.
(77), we know that $b_{1}(0)=1$. Therefore, we obtain $z=\langle A\rangle$.
Also, assuming $A\equiv\langle A\rangle$, we can conclude from Eqs. (75) and
(78) that grand-canonical partition function is written in terms of canonical
partition function as $\mathscr{Q}=e^{A\mathcal{Q}_{1}}$. To account for the
effect of correlation, we need to find the correction to $z=\langle A\rangle$.
We do this by substituting $z=\sum_{\ell=1}^{\infty}a_{\ell}\langle
A\rangle^{\ell}$ into Eq. (81) and finding $a_{\ell}$ order by order.
Afterward, we substitute the expansion of $z$ into Eq. (79) to obtain the
partition function, which is written as

$\begin{split}\log\mathscr{Q}&(\langle A\rangle,\lambda_{a})=\\\
b_{1}(\lambda_{a})&\;\langle
A\rangle+\left[b_{2}(\lambda_{a})-2b_{2}(0)b_{1}(\lambda_{a})\right]\;\langle
A\rangle^{2}+\cdots.\end{split}$ (82)

From Eq. (75), we note that the dominant term in
$\sum_{n=0}^{\infty}z^{n}\mathcal{Z}_{n}(\lambda_{a})/n!$ is around $n\approx
z$ for large $z\approx\langle A\rangle$. Therefore, we can approximately
assume $\mathscr{Q}(\langle
A\rangle,\lambda_{a})\propto\mathcal{Q}_{A}(\lambda_{a})$. The moments average
can be obtained by appropriately differentiating the grand partition function
with respect to $\lambda_{a}$ as discussed in the main text.

The introduced method is generic, and a realistic one-body density and two-
body correlation can be used to perform the cluster integrals. For particular
choice of one-body density and two-body correlation, we can perform the
integrals analytically. In particular, we assume $p_{A}(\vec{x})$ keeps the WS
distribution structure. Considering that the value of $R_{0}$ and $a_{0}$ in
WS is fixed from experimental measurements, we start with unknown parameters
$\tilde{R}_{0}$ and $\tilde{a}_{0}$, and then adjust it so that $\\{r^{2}\\}$
leads to a correct nucleus radius. This is a reasonable approximation for
small values of $d_{\text{min}}$ (see subsection. III.3). To impose two-body
correlation, we assume a hard central core for each nucleon,

$c_{2}(\vec{x}_{1},\vec{x}_{2})=1+f_{12}(\vec{x}_{1},\vec{x}_{2})=\begin{cases}0&|\vec{x}_{1}-\vec{x}_{2}|\leq
d_{\text{min}}\\\ 1&|\vec{x}_{1}-\vec{x}_{2}|>d_{\text{min}}\end{cases}.$ (83)

The details of the moments calculation can be found in appendix D. The results
are presented in Eqs. (17) and (19).

##  Appendix D Analytical moments from cluster expansion

The initial state moments can be written in terms of the following integrals:

$\displaystyle I_{m}$ $\displaystyle=\int
d^{3}x\;|\mathbf{x}|^{2m}\;p_{A}(\vec{x}),$ (84a) $\displaystyle
J^{(2)}_{m,n,\ell}$ $\displaystyle=\int
d^{3}x_{1}d^{3}x_{2}\;(x_{1}^{2}+x_{2}^{2})^{m}(y_{1}^{2}+y_{2}^{2})^{n}(x_{1}y_{1}+x_{2}y_{2})^{\ell}p_{A}(\vec{x}_{1})p_{A}(\vec{x}_{2})f_{12}(\vec{x}_{1},\vec{x}_{2}).$
(84b) $\displaystyle J^{(4)}_{m,n,\ell}$ $\displaystyle=\int
d^{3}x_{1}d^{3}x_{2}\;(x_{1}^{4}+x_{2}^{4})^{m}(y_{1}^{4}+y_{2}^{4})^{n}(x_{1}^{2}y_{1}^{2}+x_{2}^{2}y_{2}^{2})^{\ell}p_{A}(\vec{x}_{1})p_{A}(\vec{x}_{2})f_{12}(\vec{x}_{1},\vec{x}_{2}).$
(84c)

The interested moment averages is obtained from
$\partial_{\lambda_{1}}^{m_{1}}\partial_{\lambda_{2}}^{m_{2}}\mathscr{Q}(\langle
A\rangle,\lambda_{a})/\mathscr{Q}(\langle A\rangle,\lambda_{a})$ leading to
the following result:

$\displaystyle\left\langle\\{r^{2}\\}_{s,A}\right\rangle$
$\displaystyle=I_{1}+\frac{A}{2}\left(J^{(2)}_{1,0,0}+J^{(2)}_{0,1,0}-2J^{(2)}_{0,0,0}I_{1}\right),$
(85a)
$\displaystyle\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle$
$\displaystyle=\frac{I_{2}}{A}+\frac{1}{2}\left(J^{(2)}_{2,0,0}+J^{(2)}_{0,2,0}-2J^{(2)}_{1,1,0}+4J^{(2)}_{0,0,2}-2J^{(2)}_{0,0,0}I_{2}\right)$
(85b)
$\displaystyle+\left[\frac{I_{2}}{A}+\left(J^{(4)}_{1,0,0}+J^{(4)}_{0,1,0}+2J^{(4)}_{0,0,1}-2J^{(2)}_{0,0,0}I_{2}\right)\right]\sigma_{\text{fluc}}^{2}+O(J^{2}),$
$\displaystyle\left\langle\\{r^{2}\\}^{2}_{s,A}\right\rangle-\left\langle\\{r^{2}\\}_{s,A}\right\rangle^{2}$
$\displaystyle=\frac{I_{2}}{A}+\frac{1}{2}\left(J^{(2)}_{2,0,0}+J^{(2)}_{0,2,0}+2J^{(2)}_{1,1,0}-2J^{(2)}_{0,0,0}I_{2}\right)$
(85c)
$\displaystyle+\left[\frac{I_{2}}{A}+\left(J^{(4)}_{1,0,0}+J^{(4)}_{0,1,0}+2J^{(4)}_{0,0,1}-2J^{(2)}_{0,0,0}I_{2}\right)\right]\sigma_{\text{fluc}}^{2}+O(J^{2}).$

Interestingly, we have the following relation,

$\left(\left\langle\\{r^{2}\\}^{2}_{s,A}\right\rangle-\left\langle\\{r^{2}\\}_{s,A}\right\rangle^{2}\right)-\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle=2(J^{(2)}_{1,1,0}-J^{(2)}_{0,0,2})+O(J^{2},1/A),$
(86)

which merely depends on the two-body density. Ignoring the correlation, one
can use canonical partition function to calculate this difference which leads
to an important contribution from terms proportional to $1/A$.

The one-body density is given by

$p_{A}(\vec{x})\approx\frac{1}{\rho_{0}\,}\frac{\rho_{\text{WS}}(\vec{x})}{-8a_{0}^{3}\text{Li}_{3}(e^{R_{0}/a_{0}})}.$
(87)

The moments of WS distribution can be obtained analytically as follows:

$\begin{split}I_{m}(R_{0},a_{0})&=\frac{4^{m}a_{0}^{2m}\,\Gamma(m+2)\Gamma(m+1)\;\text{Li}_{3+2m}\left(-e^{R_{0}/a_{0}}\right)}{\text{Li}_{3}\left(-e^{R_{0}/a_{0}}\right)}\\\
&=\frac{3}{4}\left(\frac{\Gamma(1/2)\Gamma(m+1)}{\Gamma(m+5/2)}\right)R_{0}^{2m}\left[1+\frac{2m(2m+5)\pi^{2}}{6}\left(\frac{a_{0}}{R_{0}}\right)^{2}+\cdots\right].\end{split}$
(88)

In particular, following moments are relevant for us in this study,

$\displaystyle I_{1}(R_{0},a_{0})$
$\displaystyle=\frac{2R_{0}^{2}}{5}\left[1+\frac{7\pi^{2}}{3}\left(\frac{a_{0}}{R_{0}}\right)^{2}+\cdots\right],$
(89a) $\displaystyle I_{2}(R_{0},a_{0})$
$\displaystyle=\frac{8R_{0}^{4}}{35}\left[1+6\pi^{2}\left(\frac{a_{0}}{R_{0}}\right)^{2}+\frac{31}{3}\pi^{4}\left(\frac{a_{0}}{R_{0}}\right)^{4}+\cdots\right].$

In the following, we suggest an analytical estimation for $J_{m,n,\ell}$
integral when a minimum distance between nucleons is imposed. To calculate
$J_{m,n,\ell}$ analytically, we still need to approximate WS distribution. It
is important to note that any approximation to WS will result in a slight
error to our estimation since the integral $J_{m,n,\ell}$ is the correction to
the observable. We replace WS with a step-function
$\Theta(\tilde{R}_{0}^{2}-|\vec{x}|^{2})$ and calculate the integral in Eq.
(84c) analytically by using the integral representation of the step-function.

In cases where there is no correlation between sources, the value of $f_{ij}$
becomes zero. In such cases, the absence of correlation causes the moments
$J_{m,n,\ell}$ to be equal to zero as well. To keep the leading contribution
in calculating $J_{m,n,\ell}$, we ignore the effect of skin thickness, and
therefore WS is approximated by a step-function,
$\rho_{\text{WS}}(|\vec{x}|)\propto\theta(R_{0}^{2}-|\vec{x}_{1}|^{2})$. To
calculate a generic integral

$J(R_{0},d_{\text{min}})=\frac{1}{(4/3\pi R_{0}^{3})^{2}}\int
d^{3}x_{1}d^{3}x_{2}\;f(\mathbf{x}_{1},\mathbf{x}_{2})\;\theta(R_{0}^{2}-|\vec{x}_{1}|^{2})\theta(R_{0}^{2}-|\vec{x}_{2}|^{2})f_{12}(\vec{x}_{1},\vec{x}_{2}),$
(90)

we employ the integral representation of the step-function

$\Theta(x)=\lim_{\begin{subarray}{c}\epsilon\to 0^{+}\\\
\Lambda\to\infty\end{subarray}}\frac{1}{2\pi
i}\int_{-\Lambda}^{\Lambda}\frac{e^{i\tau x}}{\tau-i\epsilon}d\tau.$ (91)

Performing the change of variable
$\vec{x}_{1}=(\vec{x}^{\prime}_{2}-\vec{x}^{\prime}_{1})/2$ and
$\vec{x}_{2}=(\vec{x}^{\prime}_{2}+\vec{x}^{\prime}_{1})/2$, the integral is
written as

$J(R_{0},d_{\text{min}})=\lim_{\begin{subarray}{c}\epsilon\to 0^{+}\\\
\Lambda\to\infty\end{subarray}}\lim_{\begin{subarray}{c}\epsilon^{\prime}\to
0^{+}\\\
\Lambda^{\prime}\to\infty\end{subarray}}\tilde{J}(R_{0},d_{\text{min}},\Lambda,\Lambda^{\prime},\epsilon,\epsilon^{\prime}),$
(92)

where

$\tilde{J}=\frac{1}{8}\frac{1}{(2\pi i)^{2}}\int d\tau
d\tau^{\prime}\,d^{3}x_{1}^{\prime}d^{3}x_{2}^{\prime}\;\frac{1}{\tau-i\epsilon}\frac{1}{\tau^{\prime}-i\epsilon^{\prime}}\;f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})\;e^{iR_{0}^{2}(\tau+\tau^{\prime})+\frac{i(\tau-\tau^{\prime})}{2}r_{1}r_{2}\cos\Omega_{12}-\frac{i(\tau+\tau^{\prime})}{4}(r_{1}^{2}+r_{2}^{2})}f_{12}(r_{1}).$
(93)

Here, we have used the $|\vec{x}^{\prime}_{i}|^{2}=r_{i}^{2}$, and

$\begin{split}|\vec{x}_{1}|^{2}&=\frac{1}{4}\left[r_{1}^{2}+r_{2}^{2}-2r_{1}r_{2}\cos\Omega_{12}\right],\qquad|\vec{x}_{2}|^{2}=\frac{1}{4}\left[r_{1}^{2}+r_{2}^{2}+2r_{1}r_{2}\cos\Omega_{12}\right].\end{split}$
(94)

In present study, the function
$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})$ has the following form:

$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})=r_{1}^{q_{1}}r_{2}^{q_{2}}g_{1}(\Omega_{1})g_{2}(\Omega_{2}),$
(95)

where $\Omega_{i}$ represents the collective azimuthal and polar angles. Using
the following integral,

$\int_{0}^{d}dx\,x^{\alpha}e^{-ax^{2}-bx}=\frac{d^{\alpha+1}}{\alpha+1}+\cdots,$
(96)

one can perform the radial integrals along $r_{1}$ and $r_{2}$ to find

$\begin{split}\tilde{J}=-\frac{1}{16}\frac{1}{(2\pi i)^{2}}\left(\int
d\Omega_{1}g_{1}(\Omega_{1})\right)\left(\frac{d_{\text{min}}^{q_{1}+3}}{q_{1}+3}\right)\left(\int
d\Omega_{2}g_{2}(\Omega_{2})\right)\Gamma\left(\frac{q_{2}+3}{2}\right)\left(\frac{4}{i}\right)^{\frac{q_{2}+3}{2}}\mathcal{I}_{\frac{q_{2}+3}{2}},\end{split}$
(97)

where

$\mathcal{I}_{n}=\int
d\tau^{\prime}\,\frac{e^{iR_{0}^{2}\tau^{\prime}}}{\tau^{\prime}-i\epsilon^{\prime}}\mathcal{J}_{n}(\tau^{\prime}),\qquad\mathcal{J}_{n}(\tau^{\prime})=\lim_{\eta\to
0}\int
d\tau\,\frac{e^{iR_{0}^{2}\tau}}{\tau-i\epsilon}\frac{1}{\left(\tau+\tau^{\prime}-i\eta\right)^{n}},\qquad\epsilon,\epsilon^{\prime},\eta>0.$
(98)

The parameter $n$ is non-integer in our case. Here, we assume it is an integer
and find a general formula for the integral and then analytically continue the
result into fraction values of $n$.

To calculate the integral $\mathcal{J}_{n}(\tau^{\prime})$, we should close
the contour from above. In this case, the integrand has pols at
$\tau=i\epsilon$ and $\tau=-\tau^{\prime}+i\eta$ in the upper half plane,
leading to

$\begin{split}\lim_{\begin{subarray}{c}\epsilon\to 0^{+}\\\
\Lambda\to\infty\end{subarray}}\mathcal{J}_{n}(\tau^{\prime})&=\frac{2\pi
i}{\tau^{\prime n}}-2\pi
i\frac{e^{-iR_{0}^{2}\tau^{\prime}}f_{n-1}(\tau^{\prime})}{\tau^{\prime
n}},\end{split}$ (99)

where $f_{n}(\tau^{\prime})$ is a polynomial of order $n$,
$f_{n}(\tau^{\prime})=\sum_{\ell=1}^{n}\frac{(iR_{0}^{2}\tau^{\prime})^{\ell}}{(\ell-1)!}$.
Therefore, the integral $\mathcal{I}_{n}$ is written as

$\begin{split}\lim_{\begin{subarray}{c}\epsilon^{\prime}\to 0^{+}\\\
\Lambda^{\prime}\to\infty\end{subarray}}\mathcal{I}_{n}(1,1,\eta)&=\int
d\tau^{\prime}\,\frac{e^{iR_{0}^{2}\tau^{\prime}}}{\tau^{\prime}-i\epsilon^{\prime}}\left(\frac{2\pi
i}{\tau^{\prime n}}-2\pi
i\frac{e^{-iR_{0}^{2}\tau^{\prime}}f_{n-1}(\tau^{\prime})}{\tau^{\prime
n}}\right).\end{split}$ (100)

The first part of the above integral is obtained as follows:

$\lim_{\begin{subarray}{c}\epsilon^{\prime},\epsilon^{\prime\prime}\to
0^{+}\end{subarray}}\int
d\tau^{\prime}\,\frac{e^{iR_{0}^{2}\tau^{\prime}}}{\tau^{\prime}-i\epsilon^{\prime}}\left(\frac{2\pi
i}{(\tau^{\prime}-i\epsilon^{\prime\prime})^{n}}\right)=-\frac{4i^{n}\pi^{2}R_{0}^{2n}}{\Gamma(n+1)}.$
(101)

For the second part of the integral in Eq. (100), we have

$\int_{-\Lambda^{\prime}}^{\Lambda^{\prime}}d\tau^{\prime}\,\frac{1}{\tau^{\prime}-i\epsilon^{\prime}}\frac{f_{n-1}(\tau^{\prime})}{(\tau^{\prime}-i\epsilon^{\prime\prime})^{n}}=-\frac{1}{\Lambda^{\prime}}\left(-2\frac{(iR_{0}^{2})^{n-1}}{(n-2)!}\right),$
(102)

which vanishes in the limit $\Lambda^{\prime}$ to infinity. Therefore

$\lim_{\begin{subarray}{c}\epsilon^{\prime}\to 0^{+}\\\
\Lambda^{\prime}\to\infty\end{subarray}}\mathcal{I}_{n}=-\frac{4i^{n}\pi^{2}R_{0}^{2n}}{\Gamma(n+1)}.$
(103)

In the following, we calculate the moments required in Eq. (85).

Moment $J^{(2)}_{0,0,0}$: In this case, the function
$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})=1$ and consequently
$q_{1}=q_{2}=0$, and $g_{1}(\Omega_{1})=g_{2}(\Omega_{2})=1$, leading to

$J^{(2)}_{0,0,0}=-\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}.$ (104)

Moment $J^{(2)}_{1,0,0}+J^{(2)}_{0,1,0}$: To calculate this moment
combination, one needs to choose the function
$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})$ as

$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})=x_{1}^{2}+x_{2}^{2}+y_{1}^{2}+y_{2}^{2}=\frac{1}{4}\left(r_{1}^{2}(1-\cos
2\theta_{1})+r_{2}^{2}(1-\cos 2\theta_{2})\right).$ (105)

Therefore, the integral has two parts. For the first part, we have

$q_{1}=2,\quad g_{1}(\Omega_{1})=\frac{1}{4}(1-\cos 2\theta_{1}),\qquad
q_{2}=0,\quad g_{2}(\Omega_{2})=1,$

corresponding to a term proportional to $d_{\text{min}}^{5}$, which is
subleading in our calculation. For the second part, one chooses

$q_{1}=0\quad g_{1}(\Omega_{1})=1,\qquad q_{2}=2,\quad
g_{2}(\Omega_{2})=\frac{1}{4}(1-\cos 2\theta_{2}),$

leading to term $d_{\text{min}}^{3}$. As a result, we finally find,

$J^{(2)}_{1,0,0}+J^{(2)}_{0,1,0}=-2\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}\left(\frac{2R_{0}^{2}}{5}\right).$
(106)

Moment $J^{(2)}_{2,0,0}+J^{(2)}_{0,2,0}-2J^{(2)}_{1,1,0}+4J^{(2)}_{0,0,2}$:
The function $f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})$ need to be
chosen as the following:

$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})=\left(\left(x_{2}+y_{1}\right)^{2}+\left(x_{1}-y_{2}\right)^{2}\right)\left(\left(x_{2}-y_{1}\right)^{2}+\left(x_{1}+y_{2}\right)^{2}\right)$
(107)

In $r_{i}$, $\Omega_{i}$ coordinate, this function is written as

$r_{1}^{4}\left(\frac{\sin^{4}\theta_{1}}{4}\right)+\frac{1}{2}r_{1}^{2}r_{2}^{2}\left(\cos\left(2\phi_{1}-2\phi_{2}\right)\sin^{2}\theta_{1}\sin^{2}\theta_{2}\right)+r_{2}^{4}\left(\frac{\sin^{4}\theta_{2}}{4}\right).$
(108)

The last term leads to the leading term $d_{\text{min}}^{3}$ which corresponds
to

$q_{1}=0,\quad g_{1}(\Omega_{1})=1,\qquad q_{2}=4,\quad
g_{2}(\Omega_{2})=\frac{\sin^{4}\theta_{2}}{4}.$

As a result, the contribution proportional to $d_{\text{min}}^{3}$ in this
case is as follows:

$J^{(2)}_{2,0,0}+J^{(2)}_{0,2,0}-2J^{(2)}_{1,1,0}+4J^{(2)}_{0,0,2}=-4\left(\frac{8R_{0}^{4}}{35}\right)\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}.$
(109)

Moment $J^{(4)}_{1,0,0}+J^{(4)}_{0,1,0}+2J^{(4)}_{0,0,1}$: these combination
of moments corresponds to the following function:

$f(\mathbf{x}^{\prime}_{1},\mathbf{x}^{\prime}_{2})=\left(x_{1}^{2}+y_{1}^{2}\right)^{2}+\left(x_{2}^{2}+y_{2}^{2}\right)^{2}.$
(110)

Similar to the previous case, in spherical coordinate the leading contribution
is coming from the term $r_{2}^{4}\left(\frac{\sin^{4}\theta_{2}}{8}\right)$.
This means that

$J^{(4)}_{1,0,0}+J^{(4)}_{0,1,0}+2J^{(4)}_{0,0,1}=-2\left(\frac{8R_{0}^{4}}{35}\right)\left(\frac{d_{\text{min}}}{R_{0}}\right)^{3}.$
(111)

The contribution from $\sigma_{\text{fluc}}^{2}$ at the presence of non-
vanishing $d_{\text{min}}$ is given by

$J^{(4)}_{1,0,0}+J^{(4)}_{0,1,0}+2J^{(4)}_{0,0,1}-2I_{2}J^{(2)}_{0,0,0}\approx
0.$ (112)

The result in Eqs. (17) and (19) can be obtained by substituing the moments
into Eqs (85).

##  Appendix E Two-body correlation in a 1D toy model

We aim to determine the fluctuation of a pair of sources on a one-dimensional
line distributed via a one-dimensional Gaussian distribution. The distribution
is given by the equation
$p(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{x^{2}}{2\sigma^{2}}}$. We also
assume that a minimum distance is imposed among these sources. The two-body
density is defined as follows:

$\rho_{2}(x_{1},x_{2})=\frac{p(x_{1})p(x_{2})\Theta(|x_{1}-x_{2}|-d_{\text{min}})}{\int
dx_{1}dx_{2}p(x_{1})p(x_{2})\Theta(|x_{1}-x_{2}|-d_{\text{min}})},$ (113)

where $\Theta$ is the Heaviside step function. We can redefine the two-body
density as
$\rho_{2}(x_{1},x_{2})=\rho_{1}(x_{1})\rho_{1}(x_{2})c_{2}(x_{1},x_{2})$,
where the one-body density is given by:

$\rho_{1}(x_{1})=\int
dx_{2}\rho_{2}(x_{1},x_{2})=p(x_{1})\left[1+\left(\frac{1}{\sqrt{\pi}}-2\sigma
p(x_{1})\right)\left(\frac{d_{\text{min}}}{\sigma}\right)+\cdots\right],$
(114)

and the two-body correlation reads as:

$c_{2}(x_{1},x_{2})=\Theta(|x_{1}-x_{2}|-d_{\text{min}})\left[1-\left(\frac{d_{\text{min}}}{\sigma}\right)\left(\frac{1}{\sqrt{\pi}}-2\sigma
p(x_{1})-2\sigma p(x_{2})\right)+\cdots\right].$ (115)

|  |
---|---|---
Figure 9:  This figure illustrates the implementation of two-body correlation
on a 1D toy model. The blue dots represent the results obtained using a toy
Monte Carlo based on the direct method, while the red dots correspond to the
iterative method with the step-function as a two-body correlation. The dashed
curves depict the analytical calculation.

In the following, we test three scenarios of generating two sources at
positions $x_{1}$ and $x_{2}$. We analytically calculate $\langle
x_{1}^{2}+x_{2}^{2}\rangle$ and $\langle(x_{1}^{2}+x_{2}^{2})^{2}\rangle$ and
compare our calculation with toy Monte Carlo simulations.

Method 1: Direct method. We first generate two random numbers and calculate
their distance. If the distance is smaller than $d_{\text{min}}$, we discard
the pair; otherwise, we accept it. This method matches with the following
calculation,

$\displaystyle\langle x_{1}^{2}+x_{2}^{2}\rangle$ $\displaystyle=\int
dx_{1}dx_{2}(x_{1}^{2}+x_{2}^{2})\rho_{2}(x_{1},x_{2})=2\sigma^{2}+\frac{\sigma
d_{\text{min}}}{\sqrt{\pi}}+\cdots,$ (116a)
$\displaystyle\langle(x_{1}^{2}+x_{2}^{2})^{2}\rangle$ $\displaystyle=\int
dx_{1}dx_{2}(x_{1}^{2}+x_{2}^{2})^{2}\rho_{2}(x_{1},x_{2})=8\sigma^{4}+\frac{5\sigma^{3}d_{\text{min}}}{\sqrt{\pi}}+\cdots.$
(116b)

Method 2: Iterative method with $p(x)$ as one-body density and step-function
as two-body correlation.  First, we generate a random source. Then, we
generate a second source and calculate its distance from the first one. If the
distance is less than the value of $d_{\text{min}}$, we discard the second
source and generate a new one. If the distance is greater than or equal to
$d_{\text{min}}$, we accept the second source. The analytical calculation for
this method is as follows: The second source’s position depends on the first
source’s position. Therefore, we first calculate $\langle
x_{2}^{2}\rangle(x_{1})\propto\int
dx_{2}x_{2}^{2}p(x_{2})\Theta(|x_{1}-x_{2}|-d_{\text{min}})$. Then, we
calculate $\langle x_{1}^{2}+\langle x_{2}^{2}\rangle(x_{1})\rangle=\int
dx_{1}(x_{1}^{2}+\langle x_{2}^{2}\rangle(x_{1}))p(x_{1})$ using the result
from the previous step. We follow the same process for
$(x_{1}^{2}+x_{2}^{2})^{2}$,

$\displaystyle\langle x_{1}^{2}+x_{2}^{2}\rangle$
$\displaystyle=2\sigma^{2}+\frac{\sigma d_{\text{min}}}{2\sqrt{\pi}}+\cdots,$
(117a) $\displaystyle\langle(x_{1}^{2}+x_{2}^{2})^{2}\rangle$
$\displaystyle=8\sigma^{4}+\frac{7\sigma^{3}d_{\text{min}}}{4\sqrt{\pi}}+\cdots.$
(117b)

Method 3: Iterative method with $\rho_{1}(x)$ as one-body density and and
$c_{2}(x_{1},x_{2})$ as two-body correlation.  This method is similar to the
previous one while we replace $p(x)$ with $\rho_{1}(a)$ (Eq. (114)) and step-
function with $c_{2}(x_{1},x_{2})$ (Eq. (115)). The result is given by

$\displaystyle\langle x_{1}^{2}+x_{2}^{2}\rangle$
$\displaystyle=2\sigma^{2}+\frac{\sigma d_{\text{min}}}{\sqrt{\pi}}+\cdots,$
(118a) $\displaystyle\langle(x_{1}^{2}+x_{2}^{2})^{2}\rangle$
$\displaystyle=8\sigma^{4}+\frac{5\sigma^{3}d_{\text{min}}}{\sqrt{\pi}}+\cdots.$
(118b)

The above analysis shows that method 2 produces inconsistent outcomes compared
to methods 1 and 3. It can be concluded that using an iterative method with
true $\rho_{1}(x)$ and $c_{2}(x_{1},x_{2})$ leads to a result similar to that
obtained by starting directly from the two-body density. The difference
between methods can be explained as follows: In the direct method (method 1),
the most likely position for uncorrelated sources is around the origin.
However, when correlation is present, finding sources very close to the origin
becomes unlikely, and the probability of finding their position is symmetric
for both sources. In method 2, the first source is generated independently and
is usually located around the center. On the other hand, the second source
cannot be too close to the center and needs to be placed further away from the
origin. In method 3, the one-body density modifies, resulting in a
distribution that differs from the Gaussian distribution. This distribution is
already more populated symmetrically in points on both sides of the origin,
making it similar to method 1.

In Fig. 9, we have compared the toy Monte Carlo of methods 1 and 2 with our
analytical calculation, where we define $\delta r$ as follows:

$\delta
r=\frac{\langle(x_{1}^{2}+x_{2}^{2})^{2}\rangle-2\langle(x_{1}^{2}+x_{2}^{2})\rangle^{2}}{\langle(x_{1}^{2}+x_{2}^{2})\rangle^{2}}.$
(119)

In equations (116) to (118), we presented the leading terms, while in the
figure, we showed the full calculation. The figure makes it clear that there
is a more drastic effect on $\langle x_{1}^{2}+x_{2}^{2}\rangle$, which
corresponds to the size of the one-body density, compared to $\delta r$, which
captures the two-body correlation.

##  Appendix F Effect of substructure, finite cross-section, and deformation

###  F.1 Estimation of the nucleon substructure contribution

It is possible to assume substructure constituents for the nucleons. The
reduced thickness function can then be written as follows:

$T_{R}(\mathbf{x})\approx\frac{\mathcal{N}}{2}\sum_{i=1}^{N_{\text{part}}}\sum_{j=1}^{n_{c}}\gamma_{i,j}\int
dz\,\rho_{\text{c}}(\vec{x}-\vec{x}_{i}-\vec{\xi}_{i,j}),$ (120)

where $\vec{\xi}_{i,j}$ indicates the location of $j$th constituent of the
$i$th nucleon, and $n_{c}$ is the number of constituents. Here, we assume the
constituents are distributed inside the nucleons via a normal distribution
with width $\sigma_{w}$, and constituents have a round Gaussian profile with
width $v$. We have similar decomposition as that in appendix C.1,
$e^{W_{R}}=e^{W_{c}}e^{W_{\text{source}}}$.

Considering that the substructure of nucleons has a small contribution, we
will exclude the two-body correlation from our calculations. Additionally, we
will assume that the density of the nucleus is Gaussian in shape with a width
of $\sigma_{G}\approx R_{0}/\sqrt{5}$. It is important to note that the
contribution of the nucleon substructure is local, and the overall shape of
the nucleus should not have a significant impact. These assumptions will allow
us to use Eqs (24) and (25).

Suppose a nucleon with a density distribution of
$\rho_{\text{nucleon}}(\vec{x})$ and width of $\sigma_{w}$ is located at
position $\vec{x}_{i}$. Hence, the constituents are distributed based on this
distribution at locations $\vec{\xi}_{i,1},\ldots,\vec{\xi}_{i,n_{c}}$. The
fluctuation of these constituents is given by
$\rho_{\text{nucleon}}(\vec{\xi}_{i,1}-\vec{x}_{i})\cdots\rho_{\text{nucleon}}(\vec{\xi}_{i,{n_{c}}}-\vec{x}_{i})$.
However, we impose a constraint that the average location of $n_{c}$
constituents must coincide with the location of the nucleon at $\vec{x}_{i}$.
This constraint is given by

$\rho_{\text{nucleon}}(\vec{\xi}_{i,1}-\vec{x}_{i})\cdots\rho_{\text{nucleon}}(\vec{\xi}_{i,{n_{c}}}-\vec{x}_{i})\;\delta([\vec{x}]_{i}-\vec{x}_{i}),$

where $[\vec{x}]_{i}=\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}\vec{\xi}_{i,j}$ is the
location of the nucleon’s constituent center. Hence, we obtain the following
expression for the source fluctuation distribution,

$\begin{split}P_{A}=\left[\prod_{i=1}^{A}g([\vec{x}]_{i})\prod_{j=1}^{n_{c}}\rho_{\text{nucleon}}(\vec{\xi}_{i,j}-[\vec{x}]_{i})p_{\text{fluc}}(\gamma_{i,j})\right],\end{split}$
(121)

which is the generalized version of Eq. (74) at the presence of nucleon
substructure. Here, $g(\vec{x})$ is a Gaussian distribution with width
$\sigma_{G}$. Evidently, the Gaussian distribution for nucleus density is a
poor approximation. For that reason, we only single out the contribution of
the substructure to the observables and incorporate it with those obtained in
the previous parts based on WS distribution.

To isolate the effect of the substructure, it is useful to define the
following measures:

$\begin{split}\mathcal{D}\xi_{i}&=\left(d^{3}\xi_{i,1}\cdots
d^{3}\xi_{i,n_{c}}g([\vec{x}]_{i})\prod_{j=1}^{n_{c}}\rho_{\text{nucleon}}(\vec{\xi}_{i,j}-[\vec{x}]_{i})\right),\\\
\mathcal{D}\gamma&=d\gamma_{i,1}p_{\text{fluc}}(\gamma_{i,1})\cdots
d\gamma_{i,n_{c}}p_{\text{fluc}}(\gamma_{i,n_{c}}).\end{split}$ (122)

We also define the following function:

$F(\gamma_{i,j},\vec{\xi}_{i,j})=\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}\gamma_{ij}f(\vec{\xi}_{i,j}).$
(123)

Using these definitions, the moments read as

$\mathcal{O}=\frac{1}{N_{A}}\sum_{i=1}^{N_{A}}F_{i}(\gamma_{i,j},\vec{\xi}_{i,j}),$
(124)

and one can straightforwardly find

$\displaystyle\langle\mathcal{O}\rangle$
$\displaystyle=\frac{\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F(\gamma_{1,j},\vec{\xi}_{1,j})}{\int\mathcal{D}\xi},$
(125a) $\displaystyle\langle\mathcal{O}_{1}\mathcal{O}_{2}\rangle$
$\displaystyle=\frac{1}{N_{A}}\frac{\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F_{1}(\gamma_{1,j},\vec{\xi}_{1,j})F_{2}(\gamma_{1,j},\vec{\xi}_{1,j})}{\int\mathcal{D}\xi}+\frac{N_{A}-1}{N_{A}}\left(\frac{\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F_{1}(\gamma_{1,j},\vec{\xi}_{1,j})}{\int\mathcal{D}\xi}\right)\left(\frac{\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F_{2}(\gamma_{1,j},\vec{\xi}_{1,j})}{\int\mathcal{D}\xi}\right).$
(125b)

Therefore, we need to calculate the integrals over substructures. These
integrals are $n_{c}$ dimensions and cannot be decomposed into smaller
dimensions due to the
constraint$[\vec{x}]_{i}=\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}\vec{\xi}_{i,j}$. We
still, however, relate the integration over $F(\gamma_{1,j},\vec{\xi}_{1,j})$
to the integration over $f(\vec{\xi}_{i,j})$,

$\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F(\gamma_{1,j},\vec{\xi}_{1,j})=\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}\frac{1}{n_{c}}\sum_{j=1}^{n_{c}}\gamma_{1,j}f(\vec{\xi}_{1,j})=\int\mathcal{D}\xi_{1}f(\vec{\xi}_{1,1}),$
(126)

where we used the fact that partons are indistinguishable in the second
equality. Similarly,

$\begin{split}\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}F_{1}(\gamma_{1,j},\vec{\xi}_{1,j})F_{2}(\gamma_{1,j},\vec{\xi}_{1,j})&=\frac{1}{n_{c}^{2}}\int\mathcal{D}\gamma_{1}\mathcal{D}\xi_{1}\left[\sum_{j=1}^{n_{c}}\gamma_{1,j}^{2}f_{1}(\vec{\xi}_{1,j})f_{2}(\vec{\xi}_{1,j})+\sum_{j\neq
k}^{n_{c}}\gamma_{1,j}\gamma_{1,k}f_{1}(\vec{\xi}_{1,j})f_{2}(\vec{\xi}_{1,k})\right]\\\
&=\frac{1}{n_{c}}\left(1+\sigma_{\text{fluc}}^{2}\right)\int\mathcal{D}\xi_{1}\;f_{1}(\vec{\xi}_{1,1})f_{2}(\vec{\xi}_{1,1})+\frac{n_{c}-1}{n_{c}}\int\mathcal{D}\xi_{1}\;f_{1}(\vec{\xi}_{1,1})f_{2}(\vec{\xi}_{1,2}).\end{split}$
(127)

Now, we calculate the integration over $f(\vec{\xi}_{1,2})$. We are
particularly interested in the following choices,

$f_{1}(\vec{\xi})=|\bm{\xi}|^{2},\qquad
f_{\pm}(\vec{\xi})=|\bm{\xi}|^{2}e^{\pm 2i\varphi},$ (128)

where $\vec{\xi}=(\bm{\xi},\xi_{z})$, and $\bm{\xi}=(\xi,\eta)$. The first
function corresponds to observable $\\{r^{2}\\}_{s,A}$,

$\begin{split}\left\langle\\{r^{2}\\}_{s,A}\right\rangle&=\frac{\int\mathcal{D}\xi\left|\bm{\xi}\right|^{2}}{\int\mathcal{D}\xi}=\frac{2\int\mathcal{D}\xi\;\xi^{2}}{\int\mathcal{D}\xi},\end{split}$
(129)

and the functions $f_{\pm}(\vec{\xi})$ correspond to observable
$\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}$,

$\begin{split}\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle&=\frac{\left(1+\sigma_{\text{fluc}}^{2}\right)}{N_{A}n_{c}}\frac{\int\mathcal{D}\xi\left|\bm{\xi}\right|^{4}}{\int\mathcal{D}\xi}+\frac{n_{c}-1}{N_{A}n_{c}}\frac{\int\mathcal{D}\xi_{1}\left|\bm{\xi}_{1,1}\right|^{2}\left|\bm{\xi}_{1,2}\right|^{2}e^{2i(\varphi_{1,1}-\varphi_{1,2})}}{\int\mathcal{D}\xi},\\\
&=\frac{2\left(1+\sigma_{\text{fluc}}^{2}\right)}{N_{A}n_{c}}\left[\frac{\int\mathcal{D}\xi\;\xi^{4}+\int\mathcal{D}\xi\;\xi^{2}\eta^{2}}{\int\mathcal{D}\xi}\right]\\\
&+\frac{2(n_{c}-1)}{N_{A}n_{c}}\left[\frac{\int\mathcal{D}\xi\;\xi_{1}^{2}\xi_{2}^{2}-\int\mathcal{D}\xi\;\xi_{1}^{2}\eta_{2}^{2}+2\int\mathcal{D}\xi\;\xi_{1}\xi_{2}\eta_{1}\eta_{2}}{\int\mathcal{D}\xi}\right],\end{split}$
(130)

where we have used the symmetry in $\bm{\xi}=(\xi,\eta)$ space to simplify
equations. The required integrals are listed in the following,

$\begin{split}&\frac{\int\mathcal{D}\xi\,\xi^{2}}{\int\mathcal{D}\xi}=\sigma_{R}^{2},\qquad\frac{\int\mathcal{D}\xi\,\xi^{4}}{\int\mathcal{D}\xi}=3\sigma_{R}^{4},\qquad\frac{\int\mathcal{D}\xi\,\xi^{2}\eta^{2}}{\int\mathcal{D}\xi}=\frac{\int\mathcal{D}\xi\,\xi_{1}^{2}\eta_{2}^{2}}{\int\mathcal{D}\xi}=\sigma_{R}^{4},\\\
&\frac{\int\mathcal{D}\xi\,\xi_{1}^{2}\xi_{2}^{2}}{\int\mathcal{D}\xi}=3\sigma_{R}^{4}-4\sigma_{R}^{2}\sigma_{w}^{2}+2\sigma_{w}^{4},\qquad\frac{\int\mathcal{D}\xi\,\xi_{1}\xi_{2}\eta_{1}\eta_{2}\eta_{2}}{\int\mathcal{D}\xi}=\left(\sigma_{R}^{2}-\sigma_{w}^{2}\right)^{2},\end{split}$
(131)

where

$\sigma_{R}^{2}=\sigma_{G}^{2}\left[1+\frac{n_{c}-1}{n_{c}}\left(\frac{\sigma_{w}}{\sigma_{G}}\right)^{2}\right].$
(132)

Based on the above calculation, one finds,

$\displaystyle\frac{1}{2\sigma_{G}^{2}}\left.\left\langle\\{r^{2}\\}_{s,A}\right\rangle\right|_{\text{substr}}=5\left(\frac{n_{c}-1}{n_{c}}\right)\left(\frac{\sigma_{w}}{R_{0}}\right)^{2},$
(133a)
$\displaystyle\frac{A}{8\sigma_{G}^{4}}\left.\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s,A}\right|^{2}\right\rangle\right|_{\text{substr}}=\frac{\sigma_{\text{fluc}}^{2}}{n_{c}}\left(1+10\frac{n_{c}-1}{n_{c}}\left(\frac{\sigma_{w}}{R_{0}}\right)^{2}\right)+\cdots.$
(133b)

The right-hand side is the correction that needs to be incorporated with the
cluster expansion calculation.

We will now discuss the different limits of Eqs. (133). First, we will ignore
the effect of skin thickness and short-range correlation in Eq. (17). Next, we
will incorporate Eq. (133a) with Eq. (17) and implement them into Eqs. (71a)
and (72a) where we replaced $w$ with $v$ as we are specifically discussing
nucleons with substructure,

$\left\langle\\{r^{2}\\}\right\rangle=\frac{2R_{0}^{2}}{5}+2\frac{n_{c}-1}{n_{c}}\sigma_{w}^{2}+2v^{2}.$
(134)

This equation suggests defining a new effective width for the nucleons as
follows: $w^{2}=\frac{n_{c}-1}{n_{c}}\sigma_{w}^{2}+v^{2}$. In
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model, $w$
is the input parameter not $\sigma_{w}$ which can be written in terms of other
parameters via [39]

$\sigma_{w}=\sqrt{\frac{w^{2}-v^{2}}{1-1/n_{c}}}.$ (135)

In the case where there is no substructure ($n_{c}=1$), the expected result is
$w=v$, which leads to the correction term in Eq. (133b) vanishing. According
to subsection C.1, only constituents that are round in shape contribute to
$\\{r^{2}\\}$. When $n_{c}$ approaches infinity, the fluctuation contributions
from nucleon constituents vanish, and the nucleons behave like smooth round
objects again. Therefore, the correction term in Eq. (133b) also vanishes in
this limit. The constituent weight fluctuation contributions also vanish in
this limit because they are only applied to nucleon constituents. As
$\sigma_{\text{fluc}}\to 0$, the substructure contribution appears at the
order of $(\sigma_{w}/R_{0})^{4}$ in Eq. (133b).

###  F.2 Effect of finite cross-section

Previously, we considered that all nucleons in nuclei were involved in the
collision, which means $N_{\text{part}}=2A$. However, this assumption does not
take into account the mutual correlation between the sources in the target and
the sources in the projectile. When the inelastic cross-section is finite, the
correlation between the nucleons in the target and the projectile is induced,
which means that the decompression in Eq. (68) requires a correction, and Eq.
(72) must be modified.

T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
calculates the probability of colliding sources forming the target and the
projectile. This probability is given by a distribution
$P_{\text{coll}}(\bm{d})$, which depends on the size of each nucleon and the
inelastic cross-section between two nucleons. In the black-disk approximation,
$P_{\text{coll}}(\bm{d})$ is assumed to be a step-function, and collisions
occur when the distance between two nucleons is less than
$\sqrt{\sigma_{\text{fluc}}^{\text{NN}}/\pi}$. After the position of the
sources is sampled, the distance between pairs from the target and the
projectile is calculated, and they are accepted as participants based on
$P_{\text{coll}}(\bm{d})$.