# Bayesian approach to long-range correlations and multiplicity fluctuations
in nucleus-nucleus collisions

Kianusch Vahid Yousefnia  Institut de physique théorique, Université Paris
Saclay, CNRS, CEA, F-91191 Gif-sur-Yvette, France  Atharva Kotibhaskar
Department of Physics, Fergusson College, F.C. Road, Pune 411004, India
Rajeev Bhalerao  Department of Physics, Indian Institute of Science Education
and Research (IISER), Homi Bhabha Road, Pune 411008, India  Jean-Yves
Ollitrault  Institut de physique théorique, Université Paris Saclay, CNRS,
CEA, F-91191 Gif-sur-Yvette, France

###### Abstract

The number of particles detected in a nucleus-nucleus collision strongly
depends on the impact parameter of the collision. Therefore, multiplicity
fluctuations, as well as rapidity correlations of multiplicities, are
dominated by impact parameter fluctuations. We present a method based on
Bayesian inference which allows for a robust reconstruction of fluctuations
and correlations at fixed impact parameter. We apply the method to ATLAS data
on the distribution of charged multiplicity and transverse energy. We argue
that multiplicity fluctuations are smaller at large rapidity than around
central rapidity. We suggest simple, new analyses, in order to confirm this
effect.

##  I Introduction

The effective theory of strong interactions predicts that in ultrarelativistic
nucleus-nucleus collisions, particle production occurs through the formation
of color flux tubes parallel to the collision axis Lappi:2006fp , which are
analogous to strings Andersson:1983ia . Since these tubes extend over a wide
range in rapidity, one expects that particle production at different
rapidities is strongly correlated Gelis:2008sz ; Lappi:2019kif . The
measurement of this rapidity-rapidity correlation Bzdak:2012tp  has remained
elusive so far because it is hidden by the trivial correlation induced by the
variation of impact parameter $b$ within a given centrality class, since a
more central collision produces more particles at all rapidities. Methods have
been proposed to overcome this limitation Olszewski:2017vyg , but in practice,
experimental results are essentially limited to the forward-backward
correlation PHOBOS:2006mfc ; STAR:2009goo ; ATLAS:2016rbh ; ALICE:2017mtc ,
which has reduced sensitivity to impact parameter fluctuations.

Figure 1: (Color online) Normalized histogram of the distribution of
$(E_{T},N_{\rm ch})$ measured by the ATLAS collaboration Aaboud:2019sma  in
Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV. $E_{T}$ denotes the
transverse energy in the pseudorapidity window $3.2<|\eta|<4.9$, and $N_{\rm
ch}$ the number of reconstructed tracks of charged particles in the
pseudorapidity window $|\eta|<2.5$. The black ellipse is one of the results
given by our Bayesian reconstruction. It is the 99% confidence ellipse for
collisions at zero impact parameter.

We introduce a simple Bayesian method which allows us to reconstruct
correlations at fixed impact parameter, and therefore eliminates the trivial
correlation due to variations of $b$. More specifically, if
$(N_{1},N_{2},\cdots,N_{p})$ denote multiplicities, or transverse energies, in
$p$ rapidity windows, we show that their $p\times p$ covariance matrix can be
reconstructed for $b=0$ (central collisions), and a $(p-1)\times(p-1)$
projection of this matrix for $b>0$. Note that Bayesian methods in the field
of heavy-ion collisions usually involve sophisticated models and extensive
calculations Bernhard:2016tnd ; Nijs:2020roc ; JETSCAPE:2020mzn ;
Parkkila:2021tqq . By contrast, the simple Bayesian approach implemented in
this article does not involve any specific model of the collision dynamics. It
could be implemented directly by experimental collaborations.

We apply the reconstruction to ATLAS data Aaboud:2019sma  on the distribution
of the transverse energy $E_{T}$ measured in a calorimeter located at forward
rapidity, and of the charged multiplicity $N_{\rm ch}$ measured in the central
rapidity region. The histogram of $(E_{T},N_{\rm ch})$ for Pb+Pb collisions is
represented in Fig. 1. The impact parameter of a Pb+Pb collision is not
measured, and this histogram is the superposition of the contributions of all
impact parameters. However, we are able to reconstruct accurately the location
of central collisions on this diagram, which is represented in the figure as a
black curve corresponding to the 99% confidence ellipse for collisions at
$b=0$. This is a specific example of what can be achieved through the Bayesian
reconstruction.

The method is explained in Sec. II. In Sec. III, it is validated using a
realistic model calculation in which the impact parameter is known, using the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model of initial conditions
Ke:2016jrd . In Sec. IV, we apply it to ATLAS data. In Sec V, we present our
results for the fluctuations of $E_{T}$ and $N_{\rm ch}$, and discuss what
they tell us about the early collision dynamics. In Sec. VI, we suggest
analyses which could be done easily with existing data, and which would shed
new light on the rapidity dependence of fluctuations, and on long-range
correlations.

##  II Method

Detectors of heavy-ion experiments typically measure multiplicities of charged
particles or energies deposited by these particles. We assume that there are
$p$ such observables in every event, which we denote by
$(N_{1},N_{2},\cdots,N_{p})$. Experimentally, one typically measures the
probability distribution of $(N_{1},N_{2},\cdots,N_{p})$ in minimum-bias
nucleus-nucleus collisions. Our goal is to reconstruct as much as we can of
the probability distribution of $(N_{1},N_{2},\cdots,N_{p})$ at fixed impact
parameter, without relying on any specific model of the collision dynamics.
The single-variable case, $p=1$, has been studied in Das:2017ned . The present
work generalizes this earlier study to several variables. We illustrate this
generalization on ATLAS data shown in Fig. 1. These data correspond to the
case $p=2$, where $N_{1}=E_{T}$ and $N_{2}=N_{\rm ch}$. Below, we will use
either the general notation $(N_{1},N_{2})$, or the specific notation
$(E_{T},N_{\rm ch})$, depending on the context.

The key assumption of our method is that for fixed impact parameter,
fluctuations of $N_{i}$ are Gaussian. This can be viewed as a consequence of
the central limit theorem. In order for it to hold, we first need $N_{i}$ to
be large enough. In the case of the ATLAS data represented in Fig. 1, the
multiplicity $N_{\rm ch}$ exceeds 3000 in central collisions, and an even
larger number of particles contribute to the transverse energy $E_{T}$, so
that this assumption is verified. Even though the particles seen in detectors
are emitted in the last stages of the collision dynamics Busza:2018rrf , there
is a consensus that the origin of fluctuations lies in the early collision
dynamics Miller:2007ri ; Broniowski:2007nz ; Moreland:2014oya ;
Loizides:2017ack ; Bierlich:2018xfw ; Bozek:2019wyr . Now, due to the strong
Lorentz contraction at ultrarelativistic energies, the processes through which
particle production occurs at various points in the transverse plane are
causally disconnected and hence independent from each other. Therefore, the
multiplicity or transverse energy $N_{i}$ in some detector can be seen as the
sum of a large number of independent contributions, which is the condition
under which the central limit theorem applies.

More specifically, we assume that the probability distribution of
$(N_{1},N_{2},\cdots,N_{p})$ at fixed impact parameter is a multivariate
normal distribution.

$P(N_{1},...,N_{p}|b)=\frac{\exp\left(-\frac{1}{2}(N_{i}-\bar{N}_{i}(b))\Sigma^{-1}_{ij}(b)(N_{j}-\bar{N}_{j}(b))\right)}{\sqrt{(2\pi)^{p}|\Sigma(b)|}},$
(1)

where, in the exponential, we use the Einstein summation convention over the
repeated indices $i$ and $j$. In this equation, $\bar{N}_{i}$ is the mean, or
average, value of $N_{i}$, and $\Sigma_{ij}$ is the symmetric covariance
matrix:

$\displaystyle\bar{N}_{i}$ $\displaystyle=$ $\displaystyle\langle
N_{i}\rangle$ (2) $\displaystyle\Sigma_{ij}$ $\displaystyle=$
$\displaystyle\langle(N_{i}-\bar{N}_{i})(N_{j}-\bar{N}_{j})\rangle$ (3)
$\displaystyle=$ $\displaystyle\langle
N_{i}N_{j}\rangle-\bar{N}_{i}\bar{N}_{j},$ (4)

where angular brackets denote an average over events with the same impact
parameter $b$. $\Sigma^{-1}$ denotes the inverse matrix and
$\left|\Sigma\right|$ the determinant. The validity of this Gaussian
approximation will be discussed further in Sec. III.

The measured distribution is integrated over all values of impact parameter
$b$. We carry out a simple change of variables and integrate instead over the
cumulative probability distribution of $b$,

$c_{b}\simeq\frac{\pi b^{2}}{\sigma_{\rm PbPb}},$ (5)

where $\sigma_{\rm PbPb}$ is the inelastic nucleus-nucleus cross section
Das:2017ned .111Deviations from Eq. (5) only appear for very peripheral
collisions which are excluded from our study. $c_{b}$ is referred to as the
centrality fraction, or just centrality in the heavy-ion literature. The
probability distribution of $c_{b}$ is uniform in the interval $[0,1]$.
Therefore, the measured distribution is a simple integral over $c_{b}$:

$P(N_{1},\cdots,N_{p})=\int_{0}^{1}P(N_{1},\cdots,N_{p}|c_{b})dc_{b}.$ (6)

We further assume that the parameters in Eq. (1), namely, the mean values
$\bar{N}_{i}$ and the elements of the covariance matrix $\Sigma_{ij}$, are
smooth positive functions of $c_{b}$, which we choose to parametrize as the
exponential of a polynomial Das:2017ned :

$\displaystyle\bar{N}_{i}(c_{b})$ $\displaystyle=$
$\displaystyle\bar{N}_{i}(0)\exp\left(-\sum_{n=1}^{n_{\rm
max}}a_{i,n}c_{b}^{n}\right)$ (7) $\displaystyle\Sigma_{ij}(c_{b})$
$\displaystyle=$ $\displaystyle\Sigma_{ij}(0)\exp\left(-\sum_{m=1}^{m_{\rm
max}}A_{i,j,m}c_{b}^{m}\right)$ (8)

where $\bar{N}_{i}(0)$, $a_{i,n}$, $\Sigma_{ij}(0)$, $A_{i,j,m}$ are free
parameters, and $n_{\rm max}$ and $m_{\rm max}$ are the degrees of the
polynomials used to parametrize the mean and the covariance. The parameters
are fitted in such a way that the distribution (6) matches data, and the
degree of the polynomial is adjusted so as to obtain a satisfactory fit. We
carry out a standard $\chi^{2}$ fit, keeping all non-empty boxes in the
histogram of $(E_{T},N_{\rm ch})$.

As we shall see in Secs. III and IV, the simple procedure defined by Eqs. (1),
(6) and (7) allows one to obtain a good fit to $P(N_{1},\cdots,N_{p})$. The
mean values $\bar{N}_{i}(c_{b})$ returned by the fit closely match those
obtained directly by fixing the impact parameter and averaging $N_{i}$ over
events, which implies that $\bar{N}_{i}(c_{b})$ can be accurately
reconstructed from data. By contrast, one cannot reconstruct the impact
parameter dependence of the whole covariance matrix $\Sigma_{ij}(c_{b})$, and
discrepancies are expected between the best fit value, and the value
calculated directly by fixing the impact parameter. We conclude this section
by explaining which information can be reconstructed, and which cannot.

First, consider the tip of the distribution, corresponding to the largest
values of $N_{i}$. In this region, the integral in Eq. (6) is dominated by the
contribution of central collisions, that is, $c_{b}=0$. Therefore, both
$\bar{N}_{i}(0)$ and $\Sigma_{ij}(0)$ can be reconstructed using the tip of
the distribution. Looking at Fig. 1, one clearly sees that the ellipse, which
is defined by $\bar{N}_{i}(0)$ and $\Sigma_{ij}(0)$, closely fits the tip of
the distribution of $(E_{T},N_{\rm ch})$.

Figure 2: (Color online) Illustration of Eq. (10). The line is the parametric
curve $(\bar{N}_{1}(c_{b}),\bar{N}_{2}(c_{b}))$, or ridge line (see text). For
a given point $(N_{1},N_{2})$ close to the ridge line, $c_{1}$ and $c_{2}$
defined by Eq. (10) are obtained by projections on the ridge line.

Consider now the bulk of the distribution, corresponding to less central
collisions. The mean value $\bar{N}_{i}(c_{b})$ can be reconstructed from data
in the following way. Define $c(N_{i})$ as the cumulative probability
distribution of $N_{i}$:

$c(N_{i})\equiv\int_{N_{i}}^{\infty}P(N_{i})dN_{i}.$ (9)

This is the usual definition of centrality in experiment Abelev:2013qoq .
Except for very peripheral and very central collisions, $c(N_{i})$ defined by
Eq. (9) is almost equal to the inverse function of $\bar{N}_{i}(c_{b})$
Broniowski:2001ei ; Das:2017ned . Deviations are typically well below 1%. The
functions $\bar{N}_{i}(c_{b})$ define a ridge line which is displayed in Fig.
2.

By contrast, the centrality dependence of the covariance matrix cannot be
fully reconstructed. This can be understood qualitatively as follows. A
displacement parallel to the ridge line can be due either to a change in
impact parameter, or to a fluctuation at fixed impact parameter. Therefore,
only the projection of the covariance matrix orthogonal to the ridge line can
be reconstructed.

We now explain this more quantitatively in the case $p=2$. The only points
where the probability is not negligible are those close to the ridge line (see
Fig. 1). For a point $(N_{1},N_{2})$ close to the ridge line, we denote by
$c_{1}$ and $c_{2}$ the centralities defined by $N_{1}$ and $N_{2}$ (see Fig.
2):

$\displaystyle\bar{N}_{1}(c_{1})$ $\displaystyle=$ $\displaystyle N_{1}$ (10)
$\displaystyle\bar{N}_{2}(c_{2})$ $\displaystyle=$ $\displaystyle N_{2}.$ (11)

On the ridge line, $c_{1}=c_{2}$. The only sizable contribution to the
integral in Eq. (6) is from values of $c_{b}$ close to $c_{1}$ and $c_{2}$. We
then replace $N_{i}$ with $\bar{N}_{i}(c_{i})$ in Eq. (1) and expand around
$c_{b}$:

$N_{i}-\bar{N}_{i}(c_{b})=\bar{N}_{i}(c_{i})-\bar{N}_{i}(c_{b})=(c_{i}-c_{b})\bar{N}_{i}^{\prime}(c_{b}),$
(12)

where $\bar{N}_{i}^{\prime}$ is the derivative of $\bar{N}_{i}$ with respect
to $c_{b}$. Since the integrand in Eq. (6) is sizable only in a narrow range,
one can neglect the variation of $\bar{N}_{i}^{\prime}(c_{b})$ and of
$\Sigma_{ij}^{-1}(c_{b})$ in this range. Then, the integral over $c_{b}$ is a
simple Gaussian integral, which is worked out in Appendix A. In the case
$p=2$, one obtains:

$P(N_{1},N_{2})\propto\exp\left(-\frac{\delta
c^{2}}{2\sigma_{\perp}^{2}(c_{b})}\right),$ (13)

where $\delta c\equiv c_{2}-c_{1}$ (see Fig. 2), and

$\sigma_{\perp}(c_{b})\equiv\sqrt{\frac{\Sigma_{11}(c_{b})}{\bar{N}_{1}^{\prime}(c_{b})^{2}}+\frac{\Sigma_{22}(c_{b})}{\bar{N}_{2}^{\prime}(c_{b})^{2}}-2\frac{\Sigma_{12}(c_{b})}{\bar{N}_{1}^{\prime}(c_{b})\bar{N}_{2}^{\prime}(c_{b})}}.$
(14)

In writing Eq. (13), we have neglected the pre-exponential factor, which is
independent of $\delta c$. $c_{b}$ can denote any point between $c_{1}$ and
$c_{2}$, since this is a narrow range in which we neglect the variation of
$\sigma_{\perp}(c_{b})$. Up to a normalization constant,
$\sigma_{\perp}^{2}(c_{b})$ is the projection of the covariance matrix
perpendicular to the ridge line.

In the case of the ATLAS analysis, $\delta c$ represents the difference
between the centralities determined according to $E_{T}$ and $N_{\rm ch}$, and
$\sigma_{\perp}(c_{b})$ represents the width of the distribution of
$(E_{T},N_{\rm ch})$, in units of this centrality difference. As we shall see
in Sec. III, the reconstruction of $\sigma_{\perp}(c_{b})$ from the
distribution of $E_{T}$ and $N_{\rm ch}$ is robust for all $c_{b}$, while the
elements of the covariance matrix $\Sigma_{ij}(c_{b})$ are well reconstructed
only for $c_{b}=0$.

##  III Validation

Figure 3: (Color online) Normalized histogram of $(E_{T},N_{\rm ch})$ for
$10^{7}$ minimum-bias Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV
generated using the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model Ke:2016jrd
. Note that the number of events is smaller by a factor $\sim 16$ than in the
data shown in Fig. 1, and that we use a slightly finer binning in
$(E_{T},N_{\rm ch})$. The grey square at the bottom left represents the part
of the histogram which is excluded from the fit, corresponding to the most
peripheral collisions. The ellipses correspond to the 99% confidence ellipses
for fixed values of $b=0,2,4,6$ fm, calculated directly (full lines, circles
at the center) or from the Bayesian reconstruction (dashed lines, crosses at
the center).

In order to validate the reconstruction outlined in Sec. II, we test it using
fake data generated using the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model
Ke:2016jrd . This is a state of the art Monte Carlo generator of the initial
state of proton-nucleus and nucleus-nucleus collisions, which has been tuned
to reproduce several rapidity-dependent observables. It returns, for each
collision event, an entropy density profile $s(x,y,\eta)$ in three dimensions
at an early time after the collision, where $(x,y)$ are cartesian coordinates
in the transverse plane and $\eta$ is the space-time rapidity. We generate
$10^{7}$ minimum-bias Pb+Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, and we
convert the model predictions into values of $E_{T}$ and $N_{\rm ch}$, that
can be compared with the ATLAS data in Fig. 1 in the following way. We assume
that $N_{\rm ch}$ is proportional to the initial entropy, obtained by
integrating $s(x,y,\eta)$ over $(x,y)$ and over the interval $|\eta|<2.5$. We
assume that $E_{T}$ is proportional to the initial energy, obtained by
integrating $s(x,y,\eta)^{4/3}$ over $(x,y)$ and over the interval
$3.2<|\eta|<4.9$.222We thereby assume that the equation of state is conformal,
which is approximately true at the high temperatures achieved in the early
stages of the collision. Finally, we normalize $E_{T}$ and $N_{\rm ch}$ in
such a way that their mean values for central collisions coincide with those
reconstructed from ATLAS data (Sec. IV). The resulting histogram of
$(E_{T},N_{\rm ch})$ is displayed in Fig. 3. It is roughly similar to the
experimental distribution in Fig. 1.

Figure 4: (Color online) Left: Histogram of the distribution of $(E_{T},N_{\rm
ch})$ from the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model of initial
conditions Ke:2016jrd  (see Sec. III for details) in Pb+Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV for fixed values of the impact parameter, from
top to bottom: $b=0,4,8,12$ fm. Right: two-dimensional Gaussian fits to these
distributions. We indicate the $\chi^{2}$ of each fit.

In order to check the validity of the assumption of Gaussian fluctuations, Eq.
(1), we also run the model for fixed impact parameter. We have generated
$10^{5}$ events for each value of $b$. The right panels of Fig. 4 display
Gaussian fits to the left panels. The quality of the fits, measured by their
$\chi^{2}$, decreases as the impact parameter increases. The Gaussian
approximation captures the distribution except for the bottom panel,
corresponding to $b=12$ fm, where deviations are sizable. For this value of
$b$, the Gaussian extends to negative values of $E_{T}$ and $N_{\rm ch}$ which
are unphysical.333In the single-variable case this could in principle be fixed
by replacing the Gaussian with a gamma distribution Rogly:2018ddx , but we
have not found a simple way of generalizing this approach to the multi-
variable case. Therefore, the assumption of Gaussian fluctuations breaks down
for peripheral collisions, which we exclude from the fit. Excluded events are
represented as a grey square at the bottom left of Fig. 3, corresponding to
events with $E_{T}<E_{T,{\rm min}}=0.289$ TeV and $N_{\rm ch}<N_{\rm
ch,min}=194.39$.

We then fit the remaining distribution using the method outlined in Sec. II.
We choose $n_{\rm max}=3$ and $m_{\rm max}=2$ in Eq. (7). Therefore, there are
4 parameters for $\bar{E}_{T}(c_{b})$ and $\bar{N}_{\rm ch}(c_{b})$, and 3
parameters for each element of the covariance matrix, which gives a total of
17 fit parameters. Note that we need to evaluate the Pb+Pb cross section in
order to convert the impact parameter into a centrality fraction using Eq. 5.
This is done by computing the fraction $f(b)$ of events whose impact parameter
is below some value $b$, small enough that the probability of collision is
very close to unity. The cross section is then evaluated as $\sigma_{\rm
PbPb}=\pi b^{2}/f(b)$ and we obtain $\sigma_{\text{PbPb}}=(800\pm
2)\,\text{fm}^{2}$.

Figure 5: (Color online) Left to right: $\bar{N}_{\rm ch}(b)$,
$\bar{E}_{T}(b)/\bar{N}_{\rm ch}(b)$ and $\sigma_{\perp}(b)$ (defined by Eq.
(14)) versus impact parameter. The vertical grey bands indicate the limits
above which the reconstruction can no longer be trusted because peripheral
events have been excluded from the fit (see text for details). The full lines
correspond to the direct calculation, where one averages over events at fixed
$b$. The derivatives $\bar{N}_{\rm ch}^{\prime}(c_{b})$ and
$\bar{E}_{T}^{\prime}(c_{b})$ in Eq. (14) are evaluated by fitting
$\bar{N}_{\rm ch}(c_{b})$ and $\bar{E}_{T}(c_{b})$ with a smooth curve and
taking the derivative of the fit. The dashed lines in the left and middle
panels correspond to the values reconstructed using the projected
distributions $P(E_{T})$ and $P(N_{\rm ch})$, following the same method as in
Ref. Das:2017ned . The dotted lines correspond to the values reconstructed
using the distribution $P(E_{T},N_{\rm ch})$, which is done as explained in
Sec. II. The bottom panels display the ratio between the reconstructed value
and the direct calculation.

In order to assess the quality of the reconstruction, we first evaluate the
mean values of $E_{T}$ and $N_{\rm ch}$ and their covariance matrix using Eq.
(2), where one fixes the impact parameter before averaging over events. We
compare these results with those obtained by fitting the distribution as
explained in Sec. II. Since we have excluded peripheral events, we do not
expect the reconstruction to be valid if $b$ is too large. However, the
exclusion of peripheral events is done according to the values of $E_{T}$ and
$N_{\rm ch}$, not to the value of impact parameter. We expect that the
reconstruction is robust for values of $b$ such that at least 90% of events
are included. The limiting value of $b$, above which this condition is no
longer satisfied, is represented as light grey bands in Fig. 5. The dark grey
band is defined by $\bar{E}_{T}(b)<E_{T,\rm min}$.

Figure 6: (Color online) Elements of the covariance matrix as a function of
impact parameter in the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model. Top to
bottom: Variance of $N_{\rm ch}$, variance of $E_{T}$, covariance of $N_{\rm
ch}$ and $E_{T}$. The full lines correspond to the direct calculation, where
one averages over events at fixed $b$. The dotted lines correspond to the
reconstructed values.

Figure 5 shows that the mean values $\bar{N}_{\rm ch}(b)$ and $\bar{E}_{T}(b)$
are accurately reconstructed (we plot the ratio $\bar{E}_{T}(b)/\bar{N}_{\rm
ch}(b)$ rather than $\bar{E}_{T}(b)$ for reasons that will appear in Sec. IV),
as well as the width $\sigma_{\perp}(b)$ defined by Eq. (14). We have argued
on general grounds that one cannot expect to reconstruct the elements of the
covariance matrix $\Sigma_{ij}$, except for central collisions. This is
confirmed by the results displayed in Fig. 6. All three elements of the
covariance matrix are reconstructed with an error smaller than 6% for $b=0$.
However, as $b$ increases, the reconstructed values quickly deviate from the
direct calculation. This can also be seen by looking at the ellipses in Fig.
3. They are the 99% confidence ellipses (that is, ellipses that contain 99% of
the events) for several values of $b$. They are defined by

$(N_{i}-\bar{N}_{i}(b))\Sigma^{-1}_{ij}(b)(N_{j}-\bar{N}_{j}(b))=-2\ln(1-0.99).$
(15)

For $b=0$, the reconstructed ellipse almost coincides with the direct
calculation, but the discrepancy quickly increases as $b$ increases.

##  IV Application to ATLAS data

Figure 7: (Color online) Ratio of the fit to ATLAS data in Fig. 1. As in Fig.
3, the grey square at the bottom left represents the part of the histogram
which is excluded from the fit, corresponding to the most peripheral
collisions. The lines correspond to contours with $\delta c=2.6\sigma_{\perp}$
in Eq. (13), so that the region between the lines contains $\sim 99\%$ of the
events. Specifically, the upper and lower lines are the parametric curves
$\left(\bar{E}_{T}(c_{b})+2.6\sigma_{\perp}(c_{b})\bar{E}_{T}^{\prime}(c_{b}),\bar{N}_{\rm
ch}(c_{b})\right)$, and $\left(\bar{E}_{T}(c_{b}),\bar{N}_{\rm
ch}(c_{b})+2.6\sigma_{\perp}(c_{b})\bar{N}_{\rm ch}^{\prime}(c_{b})\right)$.

We finally apply the Bayesian reconstruction to the ATLAS data shown in Fig.
1. Some of the low multiplicity events correspond to photo-nuclear events
which do not contribute to the inelastic cross section. Therefore, we evaluate
the total number of inelastic collisions $N_{\rm eve}$ using the same
centrality calibration as the ATLAS collaboration, which gives $N_{\rm
eve}=158,568,641$. We exclude peripheral events using the same cuts as in the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation. Roughly 40% of the
events are excluded. Figure 7 displays the ratio of the fit to data. The fit
is excellent, except far out in the tail of the distribution. The deviation
between fit and data is at the percent level between the two lines, which are
drawn is such a way that they encompass 99% of the events.

| ATLAS | T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D
---|---|---
$\bar{E}_{T}(0)$ | 4.424(33) TeV | 4.457(33) TeV
$\bar{N}_{\rm ch}(0)$ | 3104(1) | 3103(1)
$\Sigma_{E_{T}E_{T}}(0)$ |  $0.0224(11)$ TeV2 |  $0.0663(34)$ TeV2
$\Sigma_{N_{\rm ch}N_{\rm ch}}(0)$ | $19236(759)$ | $19664(776)$
$\Sigma_{E_{T}N_{\rm ch}}(0)$ |  $14.71(28)$ TeV |  $34.62(70)$ TeV
$a_{E_{T},1}$ | $4.00$ | $4.43$
$a_{N_{\rm ch},1}$ | $3.77$ | $4.29$
$a_{E_{T},2}$ | $-1.49$ | $-1.19$
$a_{N_{\rm ch},2}$ | $-0.91$ | $-1.78$
$a_{E_{T},3}$ | $3.97$ | $4.17$
$a_{N_{\rm ch},3}$ | $3.86$ | $4.15$
$A_{E_{T},E_{T},1}$ | $4.08$ | $17.33$
$A_{N_{\rm ch},N_{\rm ch},1}$ | $7.18$ | $11.71$
$A_{E_{T},N_{\rm ch},1}$ | $7.48$ | $15.94$
$A_{E_{T},E_{T},2}$ | $1.18$ | $5.55$
$A_{N_{\rm ch},N_{\rm ch},2}$ | $-1.89$ | $-6.90$
$A_{E_{T},N_{\rm ch},2}$ | $1.32$ | $-2.14$
Table 1:  Values of fit parameters, as defined by Eq. (7), for ATLAS data
shown in Fig. 1, and for the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D
calculation shown in Fig. 3. The uncertainties on the first five parameters
are estimated using the relative difference between the direct calculation and
its reconstruction in the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model. We
omitted error bars on the polynomial coefficients in Eq. (7) because the
individual uncertainties are not reflective of the uncertainty of the
corresponding physical observable, e.g. $\bar{E}_{T}$.

The values of the fit parameters are listed in Table 1. As explained in Sec.
II, the output of the fit consists of the following information:

  * •

Impact parameter dependence of the mean transverse energy $\bar{E}_{T}$ and of
the mean charged multiplicity $\bar{N}_{\rm ch}$.

  * •

Variances of $E_{T}$ and $N_{\rm ch}$, and covariance of $E_{T}$ and $N_{\rm
ch}$, for central collisions at $b=0$.

  * •

Impact parameter dependence of the width $\sigma_{\perp}$ defined by Eq. (14).

We estimate the error on these quantities as follows. First, we check the
robustness of the fit by varying the cutoffs in $E_{T}$ and $N_{\rm ch}$, and
we exclude a fraction of events varying from 35% to 45% of the events. The
variation of the results is significantly smaller than the difference between
the direct calculation and the reconstruction in the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation. Therefore, we estimate
the relative error on our result as the maximum relative error between the
direct calculation and the reconstructed value in the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation, shown in Figs. 5 and 6.

In the remainder of this section, we discuss our results for the mean values
of $N_{\rm ch}$ and $E_{T}$. Note that these results could also have been
obtained directly from the one-dimensional projections of the two dimensional
histogram in Fig. 1, i.e., from the histograms of $N_{\rm ch}$ and $E_{T}$
alone Das:2017ned . The most important results, which involve the covariance
matrix of $N_{\rm ch}$ and $E_{T}$, will be discussed in Sec. V.

Figure 8: (Color online) Variation of the mean charged multiplicity as a
function of impact parameter. The full line is the value reconstructed using
ATLAS data. The shaded band is our estimate of the error band (see text). The
dotted line is the result of the direct T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo
3D calculation, shown as a full line in Fig. 5.  Figure 9: (Color online) Same
as Fig. 8 for the ratio of the mean transverse energy to the mean charged
multiplicity.

Figure 8 displays the mean charged multiplicity as a function of impact
parameter. For the conversion between centrality fraction and impact
parameter, we have used the value of the inelastic Pb+Pb cross section
$\sigma_{\rm PbPb}=767$ fm2 extracted from a Glauber calculation ALICE:2015juo
.444Note that this value of $\sigma_{\rm PbPb}$ differs from that returned by
the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation. We compare model and
experiment at the same impact parameter, not at the same centrality fraction.
Note that $N_{\rm ch}$ counts reconstructed tracks, and that the efficiency of
the reconstruction quickly decreases for $p_{t}<0.8$ GeV/c Aaboud:2019sma .
Therefore, the value of $N_{\rm ch}$ seen by ATLAS is only $\sim 29\%$ of that
seen by ALICE ALICE:2016fbt  in the corresponding pseudorapidity interval. For
the sake of comparison, we also show the results from the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation (corresponding to the
full lines in Fig. 5). The variation of $\bar{N}_{\rm ch}$ as a function of
$b$ is well reproduced by the model, although not perfectly.

Figure 9 displays the ratio of the mean transverse energy to the mean charged
multiplicity, as a function of impact parameter. Let us first comment on the
order of magnitude of the result for $b=0$. The ratio
$\bar{E}_{T}/\bar{N}_{\rm ch}$ can be decomposed as $(\bar{E}_{T}/\bar{N}_{\rm
FCAL})(\bar{N}_{\rm FCAL}/\bar{N}_{\rm ch})$, where $N_{\rm FCAL}$ is the
hadron multiplicity (neutral and charged) falling into the calorimeter
acceptance. We evaluate $N_{\rm FCAL}$ by integrating the pseudorapidity
spectra of charged particles ALICE:2016fbt  over the acceptance covered by the
calorimeter, multiplying by a factor $\frac{3}{2}$ to take into account
neutral particles. We then evaluate its average value at $b=0$ by
extrapolating linearly the values in the centrality intervals $0-5\%$ and
$5-10\%$, that is, $\bar{N}_{\rm FCAL}(0)\simeq\frac{3}{2}N_{\rm
FCAL}(0-5\%)-\frac{1}{2}N_{\rm FCAL}(5-10\%)$. We obtain $\bar{N}_{\rm
FCAL}(b=0)\simeq 7930$. The contribution of a hadron to the transverse energy
is roughly the transverse mass $m_{t}=\sqrt{p_{t}^{2}+m^{2}}$. Using the value
of $\bar{E}_{T}(0)$ in Table 1, we obtain $\bar{E}_{T}(0)/\bar{N}_{\rm
FCAL}(0)=0.558$ GeV, which is the expected order of magnitude for the average
transverse mass.

We next comment on the dependence of the ratio $\bar{E}_{T}/\bar{N}_{\rm ch}$
on impact parameter. ATLAS data show a non-monotonic behavior, which is not
reproduced by the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model. This
variation can be understood by decomposing again $E_{T}/N_{\rm ch}$ as
$(E_{T}/N_{\rm FCAL})(N_{\rm FCAL}/N_{\rm ch})$. The ratio $E_{T}/N_{\rm
FCAL}$ is the average transverse mass, which is determined by the mean
transverse momentum $\langle p_{t}\rangle$, which itself decreases mildly as a
function of impact parameter ALICE:2018hza . This effect is responsible for
the decrease seen in the model, and in data at small $b$. The ratio $N_{\rm
FCAL}/N_{\rm ch}$ is the ratio of the multiplicity in the forward rapidity
region, covered by the calorimeter, and the central rapidity region, where
$N_{\rm ch}$ is measured. Now, the rapidity distribution becomes slightly
broader as the impact parameter increases ALICE:2016fbt . The intuitive
picture for this phenomenon is that the stopping between the nuclei is not as
strong as in central collisions. This implies that the ratio $N_{\rm
FCAL}/N_{\rm ch}$ increases as a function of impact parameter. This effect is
not reproduced by the T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model. It
overrides the decrease of $\langle p_{t}\rangle$ for large values of $b$,
leading to the increase seen in ATLAS data.

##  V Fluctuations at large and central rapidities

In this Section, we present results involving the covariance matrix of $E_{T}$
and $N_{\rm ch}$. These are our most important results, as they shed light on
the fluctuations in the central rapidity region probed by $N_{\rm ch}$, in the
large rapidity region probed by $E_{T}$, and on their mutual correlations. We
first isolate dynamical fluctuations of $E_{T}$ and $N_{\rm ch}$ by
subtracting out Poisson fluctuations in Sec. V.1. We then present our results
for central collisions in Sec. V.2, and finally the results on the impact
parameter dependence in Sec. V.3.

###  V.1 Subtraction of Poisson fluctuations

We generally consider quantities of the form $N=\sum_{i=1}^{M}x_{i}$, where
$M$ is the particle multiplicity. Our goal is to identify non-trivial
fluctuations and correlations of these quantities. The baseline is the case
where the probabilities of finding a particle in different regions of phase
space are independent variables DiFrancesco:2016srj . In this case, $x_{i}$
are independent variables and $M$ follows a Poisson distribution. A simple
calculation shows that the variance of $N$ is the expectation value of
$\sum_{i=1}^{M}x_{i}^{2}$. We refer to this contribution loosely as “the
variance of Poisson fluctuations”, and we subtract it from the observed
variance in order to isolate the non-trivial part.

For the charged multiplicity, one simply counts particles, so that $x_{i}=1$,
and the subtraction is straightforward:

$\Sigma_{N_{\rm ch}N_{\rm ch}}(b)\rightarrow\Sigma_{N_{\rm ch}N_{\rm
ch}}(b)-\bar{N}_{\rm ch}(b).$ (16)

There is however an uncertainty on this subtraction due to hadronic decays. If
a hadron decays into two charged particles, which both fall into the detector
acceptance, then the contribution of that hadron is $x_{i}=2$, not $x_{i}=1$.
The magnitude of this effect is modest, and cannot be evaluated accurately,
because it depends on the hadronization mechanism. In the case of a fluid-
dynamical model, it depends on the freeze-out temperature Mazeliauskas:2018irt
. We include this uncertainty in our error bar, by multiplying $\bar{N}_{\rm
ch}(b)$ in Eq. (16) by a coefficient which can vary between $1$ and $1.2$.

The transverse energy can be written as $E_{T}=\sum_{i}x_{i}$, where the sum
has $N_{\rm FCAL}$ terms, corresponding to all hadrons, neutral and charged,
falling in the calorimeter, and $x_{i}\simeq m_{t}$, where
$m_{t}=\sqrt{p_{t}^{2}+m^{2}}$ is the transverse mass. The average value of
$E_{T}$ at fixed $b$ is $\bar{E}_{T}(b)=\langle m_{t}\rangle\bar{N}_{\rm
FCAL}(b)$, and the variance of Poisson fluctuations is $\langle
m_{t}^{2}\rangle\bar{N}_{\rm FCAL}(b)$, where angular brackets denote an
average value over hadrons in the calorimeter. This quantity can be decomposed
as

$\langle m_{t}^{2}\rangle\bar{N}_{\rm FCAL}(b)=\frac{\langle
m_{t}^{2}\rangle}{\langle m_{t}\rangle^{2}}\frac{\bar{N}_{\rm
ch}(b)}{\bar{N}_{\rm FCAL}(b)}\frac{\bar{E}_{T}(b)^{2}}{\bar{N}_{\rm ch}(b)}.$
(17)

We evaluate the first factor on the right-hand side using ALICE data on
identified particle spectra ALICE:2013mez , which give $\langle
m_{t}^{2}\rangle/\langle m_{t}\rangle^{2}\simeq 1.52$ for central collisions
at $\sqrt{s_{\rm NN}}=2.76$ TeV. We neglect the dependence of this ratio on
impact parameter, rapidity and $\sqrt{s_{\rm NN}}$. The second factor can be
estimated for central collisions using the value of $\bar{N}_{\rm ch}(0)$ in
Table 1, and the estimate $\bar{N}_{\rm FCAL}(0)\simeq 7930$ obtained in Sec.
IV. We also neglect its dependence on rapidity and impact parameter. Putting
these factors together, the subtraction of Poisson fluctuations is done
according to the formula:

$\Sigma_{E_{T}E_{T}}(b)\rightarrow\Sigma_{E_{T}E_{T}}(b)-\alpha\frac{\bar{E}_{T}(b)^{2}}{\bar{N}_{\rm
ch}(b)},$ (18)

where $\alpha\simeq 0.59$. We assign an uncertainty of $\pm 20\%$ to this
factor $\alpha$.

Finally, since $E_{T}$ and $N_{\rm ch}$ are measured in separate rapidity
regions, no hadron contributes simultaneously to $E_{T}$ and $N_{\rm ch}$, and
Poisson fluctuations do not contribute to the covariance $\Sigma_{E_{T}N_{\rm
ch}}(b)$.

Using the values in Table 1, for central collisions, the subtracted quantities
in Eqs. (16) and (18) are $16\%$ and $17\%$ of the total. This implies that at
least $83\%$ of the variances of $E_{T}$ and $N_{\rm ch}$ can be attributed to
dynamical fluctuations.

###  V.2 Fluctuations of $N_{\rm ch}$ and $E_{T}$ in central collisions

| $E_{T}$ | $N_{\rm ch}$
---|---|---
| $1.14(6)\times 10^{-3}$ | $1.07(2)\times 10^{-3}$
$E_{T}$ | $\mathbf{0.95(7)\times 10^{-3}}$ | $\mathbf{1.07(2)\times 10^{-3}}$
| $\mathit{(3.57\times 10^{-3})}$ | $\mathit{(2.52\times 10^{-3})}$
| $1.07(2)\times 10^{-3}$ | $2.00(8)\times 10^{-3}$
$N_{\rm ch}$ | $\mathbf{1.07(2)\times 10^{-3}}$ | $\mathbf{1.68(9)\times 10^{-3}}$
| $\mathit{(2.52\times 10^{-3})}$ | $\mathit{(1.96\times 10^{-3})}$
Table 2:  Relative covariance matrix $\sigma_{ij}$ of $E_{T}$ and $N_{\rm
ch}$, defined by Eq. (19), for central Pb+Pb collisions ($b=0$). For each
element, the first line is the total covariance returned by the fit to data.
The second line is the dynamical covariance obtained after subtracting Poisson
fluctuations (see Sec. V.1). The third line is the value obtained in the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation by simulating events at
$b=0$, which differs by a few percent from the reconstructed value in Table 1.
The T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation returns a continuous
density profile, and does not include hadronization. Therefore, the
fluctuations calculated in this model are dynamical, and the subtraction of
Poisson fluctuations explained in Sec. V.1 does not apply.

We now discuss the fluctuations of $E_{T}$ and $N_{\rm ch}$ for central
collisions, and their mutual correlation. We use as a measure of these
fluctuations the relative covariance matrix Pruneau:2002yf , defined as

$\sigma_{ij}\equiv\frac{\langle N_{i}N_{j}\rangle}{\langle N_{i}\rangle\langle
N_{j}\rangle}-1=\frac{\Sigma_{ij}}{\bar{N}_{i}\bar{N}_{j}}.$ (19)

Its diagonal elements are the relative variances of $N_{i}$. Values extracted
from ATLAS data (Table 2) show that the relative fluctuations of $E_{T}$ are
smaller than those of $N_{\rm ch}$, both before and after isolating dynamical
fluctuations. This is the reason why $E_{T}$ is a better estimator of the
centrality than $N_{\rm ch}$ Aaboud:2019sma ; Zhou:2018fxx . It is a non-
trivial observation, which is not reproduced by the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation. The
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation reproduces the variance
of $N_{\rm ch}$ to a good approximation, which is not surprising since it was
fitted to the distribution of the charged multiplicity measured near mid-
rapidity. On the other hand, it overestimates the variance of $E_{T}$ by at
least a factor 3.

There is no trivial relation between the relative fluctuations of $E_{T}$ and
the relative fluctuations of the multiplicity in the same rapidity window. In
a hydrodynamic picture of the collision Ollitrault:2007du , however, one
expects them to be of the same order of magnitude, because the multiplicity is
proportional to the entropy, and the energy and the entropy are related
through the equation of state. In Sec. VI, we will suggest a specific analysis
in order to check this assumption. If it is correct, our results imply that
multiplicity fluctuations are smaller at large rapidity than around central
rapidity.

The values of the relative variance of $E_{T}$ in Table 2 can be compared with
those previously extracted from ATLAS data at $2.76$ TeV Das:2017ned . After
subtracting Poisson fluctuations with help of multiplicity densities measured
by ALICE ALICE:2013jfw , we obtain $\sigma_{11}=0.76\times 10^{-3}$. This is
smaller than the value at $5.02$ TeV in Table 2, which implies that
multiplicity fluctuations increase as a function of the rapidity gap between
incoming nuclei and the detector, which increases by $\sim 0.6$ between $2.76$
TeV and $5.02$ TeV. It is tempting to postulate that multiplicity fluctuations
depend on the rapidity difference with incoming nuclei, much as average
multiplicities themselves PHOBOS:2005zhy ; Gelis:2006tb .

Figure 10: (Color online) Full line: $99\%$ confidence ellipse of zero impact
parameter collisions, already shown in Fig. 1. The ellipse is defined by Eq.
(15). The dashed line represents the ellipse, after subtraction of the
contribution of Poisson fluctuations. The dotted line is the value from the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D calculation.

The second important observation is that $E_{T}$ and $N_{\rm ch}$ are strongly
correlated, even at fixed impact parameter. This is illustrated by the
elongated covariance ellipse represented in Fig. 10. The correlation is
significantly stronger after Poisson fluctuations have been subtracted out.

We now interpret these results within a simple model. We assume that a
nucleus-nucleus collision at impact parameter $b$ produces $N_{s}$ sources,
and that $E_{T}$ and $N_{\rm ch}$ are obtained by summing the contributions of
all sources in the event Bozek:2017qir ; Zhou:2018fxx ; Jia:2020tvb . We
assume that the sources are independent.555 Note that a similar picture
underlies our assumption of Gaussian fluctuations, as explained in Sec. II. We
denote the contribution of a single source to $N_{1}=E_{T}$ and $N_{2}=N_{\rm
ch}$ by $n_{1}$ and $n_{2}$, respectively. Under these assumptions, one easily
obtains:

$\sigma_{ij}=\frac{\langle N_{s}^{2}\rangle-\langle N_{s}\rangle^{2}}{\langle
N_{s}\rangle^{2}}+\frac{1}{\langle N_{s}\rangle}\frac{\langle
n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle}{\langle
n_{i}\rangle\langle n_{j}\rangle}$ (20)

The first term in the right-hand side is the relative variance of the number
of sources, which gives the same contribution to all elements of the
covariance matrix. It is a trivial correlation arising from a global change in
the system, much as correlation stemming from the variation of impact
parameter Kovchegov:2012nd . It does not shed light on rapidity correlations
in the emission process. The last term in the right-hand side, which involves
the covariance matrix of a single source, is physically more interesting.

It is not possible to uniquely determine the quantities in the right-hand side
of Eq. (20), because there are more unknowns than equations. But one can
obtain non-trivial information. For instance, we have noted that the relative
variance of $E_{T}$ is smaller than that of $N_{\rm ch}$. Since they get the
same contribution from the variance of the number of sources, it implies that
the difference between the relative variances of $E_{T}$ and $N_{\rm ch}$ is
larger at the level of a single source. The suppression of fluctuations at
large rapidity is by no means a small effect. One can also obtain non-trivial
information on the correlation between $n_{1}$ and $n_{2}$, which represent
the contributions of a single source to $E_{T}$ and $N_{\rm ch}$. For
instance, one can readily exclude that they are uncorrelated, because this
would imply $\sigma_{11}>\sigma_{12}$, at variance with the values in Table 2,
after Poisson fluctuations have been subtracted out. One can obtain a lower
bound $r_{\rm min}$ on the Pearson correlation coefficient $r$ between $n_{1}$
and $n_{2}$, defined by $\langle n_{1}n_{2}\rangle-\langle n_{1}\rangle\langle
n_{2}\rangle=r\sigma(n_{1})\sigma(n_{2})$, where
$\sigma(n_{i})\equiv\sqrt{\langle n_{i}^{2}\rangle-\langle n_{i}\rangle^{2}}$
is the standard deviation of $n_{i}$. Elementary algebra shows that
$r^{2}-r_{\rm min}^{2}$ is a perfect square, where $r_{\rm min}$ is defined
by:

$r_{\rm min}\equiv
2\frac{\sqrt{(\sigma_{12}-\sigma_{11})(\sigma_{22}-\sigma_{12})}}{\sigma_{22}-\sigma_{11}}$
(21)

With the values in Table 2, and taking into account the error bars (from the
reconstruction and from the subtraction of Poisson fluctuations) we obtain
$r_{\rm min}=0.72\pm 0.15$. Two effects may contribute to this strong
correlation: Some sources may be stronger than others, for instance those
situated in the center of the interaction region. Stronger sources yield
larger values of both $n_{1}$ and $n_{2}$, which induces a mutual correlation.
The second effect is the dynamical effect that one would like to isolate, that
particle production from a single source is strongly correlated across
rapidities. More detailed modeling will be necessary to disentangle the
relative contributions of these two effects.

Figure 11: (Color online) Same as Fig. 8 for the width of the distribution,
defined by Eq. (14). Full line: before subtracting Poisson fluctuations.
Dashed line: after subtracting Poisson fluctuations according to Eqs. (18) and
(16). The shaded area around the curve is our estimate of the error on the
reconstruction. Dotted line: value calculated in the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model.

###  V.3 Impact parameter dependence of fluctuations

As explained in Sec. II, one cannot reconstruct the impact parameter
dependence of the full covariance matrix. Only the specific linear combination
$\sigma_{\perp}$, defined by Eq. (14), can be reconstructed for all impact
parameters. Its impact parameter dependence is represented in Fig. 11, before
and after subtracting Poisson fluctuations (Sec. V.1). The
T$\mathrel{\raisebox{-2.1pt}{R}}$ENTo 3D model underpredicts $\sigma_{\perp}$
for all values of impact parameter, which is the reason why the distribution
in Fig. 3 is narrower than that in Fig. 1. This shows that the width of the
distribution, as measured by $\sigma_{\perp}$, contains non-trivial dynamical
information, which can be used to discriminate between models.

##  VI Conclusions and perspectives

We have introduced a simple Bayesian method which allows for a robust
reconstruction of multiplicity fluctuations and rapidity correlations in
nucleus-nucleus collisions at fixed impact parameter $b$. Starting from the
distribution $p(N_{1},\cdots,N_{p})$ of multiplicities (or transverse
energies, or number of hits) $N_{1},\cdots,N_{p}$ measured in different parts
of the detector, one can reconstruct the full covariance matrix of
$(N_{1},\cdots,N_{p})$ at $b=0$, as well as the impact parameter dependence of
a $(p-1)\times(p-1)$ projection of the matrix. We have applied the method to
ATLAS data on the joint distribution of transverse energy and charged
multiplicity. We have shown that dynamical fluctuations are smaller at large
rapidity than around central rapidity, and that particle production is
strongly correlated across rapidity, even at fixed impact parameter. We have
also shown that the width of the distribution and its impact parameter
dependence are not reproduced by our model calculation, so that this
observable can be used to rule out three-dimensional models of initial
fluctuations.

Rapidity correlations have so far been studied using multiplicities of charged
particle tracks STAR:2009goo ; ATLAS:2016rbh . However, particle tracks are
typically reconstructed only around central rapidity, so that the rapidity
coverage of such studies is limited. We have shown that track reconstruction
is not needed. The transverse energy in a calorimeter, the multiplicity of
pixel clusters CMS:2013bza  or the number of hits in a scintillator
ALICE:2013hur  also give a quantitative information on fluctuations, provided
that the detector is large enough that dynamical fluctuations dominate over
Poisson fluctuations. A few tests should however be made in order to interpret
correctly the information from calorimeters. Specifically, one should first
figure out how energy fluctuations relate to multiplicity fluctuations in the
same rapidity window. This can be done easily using one of the central
detectors of ALICE, CMS or ATLAS. These detectors measure the multiplicity of
charged particle tracks and their momenta. One can use the sum of transverse
momenta of charged particles, $P_{T}\equiv\sum_{N_{\rm ch}}p_{t}$, as a proxy
for their transverse energy. By comparing the tails of the distributions of
$N_{\rm ch}$ and $P_{T}$, one can readily check how their relative variances
compare for central collisions.

We conclude by listing analyses that would shed additional light on long-range
correlation. The transverse energy is usually defined as $E_{T}=E_{F}+E_{B}$,
where $E_{F}$ and $E_{B}$ are transverse energies in forward and backward
calorimeters CMS:2013bza ; ATLAS:2016rbh . If measured, the distribution
$p(E_{F},E_{B})$ would provide direct information on forward-backward rapidity
correlations. This analysis would be even simpler than that carried out in
this paper because backward and forward calorimeters are symmetric around mid-
rapidity, and this symmetry reduces the number of fit parameters in the
Bayesian analysis. Similarly, the ALICE collaboration could measure the
distribution $p(V_{0A},V_{0C})$, where $V_{0A}$ and $V_{0C}$ are the
multiplicities in backward and forward scintillators Abelev:2013qoq . These
analyses can be readily carried out by experimental collaborations using
existing data, by following the exact same steps as our analysis of ATLAS
public data.

More detailed information on the rapidity structure could be obtained by
extending the analysis to three variables. For instance, in the case of the
ATLAS data studied in this paper Aaboud:2019sma , one could split the
calorimeter into its forward and backward components and measure
$p(E_{F},E_{B},N_{\rm ch})$. Similar analyses could be done with the CMS and
ALICE detectors. They would yield detailed information on the long-range
rapidity structure of correlations. Finally, the same method could be applied
to heavy-ion collisions at lower energies INDRA:2020kyj ; Parfenov:2021ipw ;
Li:2021plq .

## Acknowledgements

We thank J. Jia for sending us the centrality calibration of data. AK thanks
Bhavya Bhatt for technical assistance in the early stages of this work. JYO
thanks F. Gelis and G. Giacalone for discussions.

##  Appendix A Reconstruction of the covariance matrix in non-central
collisions

In this Appendix, we specify which information about the covariance matrix
$\Sigma_{ij}(b)$ can be reconstructed for $b>0$. Inserting Eq. (12) into Eq.
(1) and neglecting the pre-exponential factor, one obtains:

$P(N_{1},...,N_{p}|b)\propto\exp\left(-\frac{1}{2}\sum\limits_{ij}(c_{i}-c_{b})\bar{N}^{\prime}_{i}\Sigma^{-1}_{ij}\bar{N}^{\prime}_{j}(c_{j}-c_{b})\right),$
(22)

where we omit the dependence of $\bar{N}^{\prime}_{i}$ and $\Sigma^{-1}_{ij}$
on $b$ for simplicity. Then, the integral in Eq. (6) is a Gaussian integral
over $c_{b}$. Evaluating this integral, and neglecting again the pre-
exponential factor, one obtains

$P(N_{1},...,N_{p})\propto\exp\left(-\frac{1}{2}\sum_{ij}c_{i}\Pi_{ij}c_{j}\right),$
(23)

where

$\Pi_{ij}\equiv\bar{N}^{\prime}_{i}\Sigma^{-1}_{ij}\bar{N}^{\prime}_{j}-\frac{\left(\sum_{\alpha}\bar{N}^{\prime}_{i}\Sigma^{-1}_{i\alpha}\bar{N}^{\prime}_{\alpha}\right)\left(\sum_{\beta}\bar{N}^{\prime}_{j}\Sigma^{-1}_{j\beta}\bar{N}^{\prime}_{\beta}\right)}{\sum\limits_{\alpha\beta}\bar{N}^{\prime}_{\alpha}\Sigma^{-1}_{\alpha\beta}\bar{N}^{\prime}_{\beta}}.$
(24)

$\mathbf{\Pi}$ is the quantity which can be reconstructed for all impact
parameters. It is a $(p\times p)$ symmetric matrix which verifies the property
$\sum_{i}\Pi_{ij}=0$ for all $j$, that is, all lines and columns sum up to
zero. Due to this property, one can rewrite Eq. (23) as:

$P(N_{1},...,N_{p})\propto\exp\left(-\frac{1}{2}\sum_{i>1,j>1}(c_{i}-c_{1})\Pi_{ij}(c_{j}-c_{1})\right).$
(25)

This form shows that the distribution of $(N_{1},...,N_{p})$ is a function of
$p-1$ variables $c_{i}-c_{1}$. Physically, $\mathbf{\Pi}$ represents the
projection of the covariance matrix onto the $(p-1)$-dimensional subspace
orthogonal to the ridge line. In the case $p=2$, Eq. (24) gives

$\mathbf{\Pi}=\frac{1}{\sigma_{\perp}^{2}}\begin{pmatrix}\phantom{-}1&-1\\\
-1&\phantom{-}1\end{pmatrix},$ (26)

where $\sigma_{\perp}$ is defined by Eq. (14). Inserting Eq. (26) into Eq.
(25), one recovers Eq. (13).

## References

  * (1) T. Lappi and L. McLerran, Nucl. Phys. A 772, 200-212 (2006) [arXiv:hep-ph/0602189 [hep-ph]].
  * (2) B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, Phys. Rept. 97, 31-145 (1983)
  * (3) F. Gelis, T. Lappi and R. Venugopalan, Phys. Rev. D 79, 094017 (2009) [arXiv:0810.4829 [hep-ph]].
  * (4) T. Lappi and A. Ramnath, Phys. Rev. D 100, no.5, 054003 (2019) [arXiv:1904.00782 [hep-ph]].
  * (5) A. Bzdak and D. Teaney, Phys. Rev. C 87, no.2, 024906 (2013) [arXiv:1210.1965 [nucl-th]].
  * (6) A. Olszewski and W. Broniowski, Phys. Rev. C 96, no.5, 054903 (2017) [arXiv:1706.02862 [nucl-th]].
  * (7) B. B. Back et al. [PHOBOS], Phys. Rev. C 74, 011901 (2006) [arXiv:nucl-ex/0603026 [nucl-ex]].
  * (8) B. I. Abelev et al. [STAR], Phys. Rev. Lett. 103, 172301 (2009) [arXiv:0905.0237 [nucl-ex]].
  * (9) M. Aaboud et al. [ATLAS], Phys. Rev. C 95, no.6, 064914 (2017) [arXiv:1606.08170 [hep-ex]].
  * (10) S. Acharya et al. [ALICE], Phys. Lett. B 781, 20-32 (2018) [arXiv:1710.07975 [nucl-ex]].
  * (11) M. Aaboud et al. [ATLAS], JHEP 01, 051 (2020) [arXiv:1904.04808 [nucl-ex]].
  * (12) J. E. Bernhard, J. S. Moreland, S. A. Bass, J. Liu and U. Heinz, Phys. Rev. C 94, no.2, 024907 (2016) [arXiv:1605.03954 [nucl-th]].
  * (13) G. Nijs, W. van der Schee, U. Gürsoy and R. Snellings, Phys. Rev. C 103, no.5, 054909 (2021) [arXiv:2010.15134 [nucl-th]].
  * (14) D. Everett et al. [JETSCAPE], Phys. Rev. C 103, no.5, 054904 (2021) [arXiv:2011.01430 [hep-ph]].
  * (15) J. E. Parkkila, A. Onnerstad and D. J. Kim, [arXiv:2106.05019 [hep-ph]].
  * (16) W. Ke, J. S. Moreland, J. E. Bernhard and S. A. Bass, Phys. Rev. C 96, no.4, 044912 (2017) [arXiv:1610.08490 [nucl-th]].
  * (17) S. J. Das, G. Giacalone, P. A. Monard and J. Y. Ollitrault, Phys. Rev. C 97, no.1, 014905 (2018) [arXiv:1708.00081 [nucl-th]].
  * (18) W. Busza, K. Rajagopal and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339-376 (2018) [arXiv:1802.04801 [hep-ph]].
  * (19) M. L. Miller, K. Reygers, S. J. Sanders and P. Steinberg, Ann. Rev. Nucl. Part. Sci. 57, 205-243 (2007) [arXiv:nucl-ex/0701025 [nucl-ex]].
  * (20) W. Broniowski, M. Rybczynski and P. Bozek, Comput. Phys. Commun. 180, 69-83 (2009) [arXiv:0710.5731 [nucl-th]].
  * (21) J. S. Moreland, J. E. Bernhard and S. A. Bass, Phys. Rev. C 92, no.1, 011901 (2015) [arXiv:1412.4708 [nucl-th]].
  * (22) C. Loizides, J. Kamin and D. d’Enterria, Phys. Rev. C 97, no.5, 054910 (2018) [erratum: Phys. Rev. C 99, no.1, 019901 (2019)] [arXiv:1710.07098 [nucl-ex]].
  * (23) C. Bierlich, G. Gustafson, L. Lönnblad and H. Shah, JHEP 10, 134 (2018) [arXiv:1806.10820 [hep-ph]].
  * (24) P. Bożek, W. Broniowski, M. Rybczynski and G. Stefanek, Comput. Phys. Commun. 245, 106850 (2019) [arXiv:1901.04484 [nucl-th]].
  * (25) B. Abelev et al. [ALICE], Phys. Rev. C 88, no.4, 044909 (2013) [arXiv:1301.4361 [nucl-ex]].
  * (26) W. Broniowski and W. Florkowski, Phys. Rev. C 65, 024905 (2002) [arXiv:nucl-th/0110020 [nucl-th]].
  * (27) R. Rogly, G. Giacalone and J. Y. Ollitrault, Phys. Rev. C 98, no.2, 024902 (2018) [arXiv:1804.03031 [nucl-th]].
  * (28) J. Adam et al. [ALICE], Phys. Rev. Lett. 116, no.22, 222302 (2016) [arXiv:1512.06104 [nucl-ex]].
  * (29) J. Adam et al. [ALICE], Phys. Lett. B 772, 567-577 (2017) [arXiv:1612.08966 [nucl-ex]].
  * (30) S. Acharya et al. [ALICE], Phys. Lett. B 788, 166-179 (2019) [arXiv:1805.04399 [nucl-ex]].
  * (31) P. Di Francesco, M. Guilbaud, M. Luzum and J. Y. Ollitrault, Phys. Rev. C 95, no.4, 044911 (2017) [arXiv:1612.05634 [nucl-th]].
  * (32) A. Mazeliauskas, S. Floerchinger, E. Grossi and D. Teaney, Eur. Phys. J. C 79, no.3, 284 (2019) [arXiv:1809.11049 [nucl-th]].
  * (33) B. Abelev et al. [ALICE], Phys. Rev. C 88, 044910 (2013) [arXiv:1303.0737 [hep-ex]].
  * (34) C. Pruneau, S. Gavin and S. Voloshin, Phys. Rev. C 66, 044904 (2002) [arXiv:nucl-ex/0204011 [nucl-ex]].
  * (35) M. Zhou and J. Jia, Phys. Rev. C 98, no.4, 044903 (2018) [arXiv:1803.01812 [nucl-th]].
  * (36) J. Y. Ollitrault, Eur. J. Phys. 29, 275-302 (2008) [arXiv:0708.2433 [nucl-th]].
  * (37) E. Abbas et al. [ALICE], Phys. Lett. B 726, 610-622 (2013) [arXiv:1304.0347 [nucl-ex]].
  * (38) B. B. Back et al. [PHOBOS], Phys. Rev. C 74, 021901 (2006) [arXiv:nucl-ex/0509034 [nucl-ex]].
  * (39) F. Gelis, A. M. Stasto and R. Venugopalan, Eur. Phys. J. C 48, 489-500 (2006) [arXiv:hep-ph/0605087 [hep-ph]].
  * (40) P. Bozek and W. Broniowski, Phys. Rev. C 97, no.3, 034913 (2018) [arXiv:1711.03325 [nucl-th]].
  * (41) J. Jia, C. Zhang and J. Xu, Phys. Rev. Res. 2, no.2, 023319 (2020) [arXiv:2001.08602 [nucl-th]].
  * (42) Y. V. Kovchegov and D. E. Wertepny, Nucl. Phys. A 906, 50-83 (2013) [arXiv:1212.1195 [hep-ph]].
  * (43) S. Chatrchyan et al. [CMS], JHEP 02, 088 (2014) [arXiv:1312.1845 [nucl-ex]].
  * (44) B. Abelev et al. [ALICE], Phys. Rev. C 88, no.4, 044909 (2013) [arXiv:1301.4361 [nucl-ex]].
  * (45) J. D. Frankland et al. [INDRA], [arXiv:2011.04496 [nucl-ex]].
  * (46) P. Parfenov, D. Idrisov, V. B. Luong and A. Taranenko, Particles 4, no.2, 275-287 (2021)
  * (47) F. Li, Y. Wang, Z. Gao, P. Li, H. Lv, Q. Li, C. Y. Tsang and M. B. Tsang, [arXiv:2105.08912 [nucl-th]].