# Controlling volume fluctuations for studies of critical phenomena in nuclear
collisions

Romain Holzmann<EMAIL_ADDRESS>GSI Helmholtzzentrum für
Schwerionenforschung, 64291 Darmstadt, Germany  Volker Koch<EMAIL_ADDRESS>Nuclear Science Division Lawrence Berkeley National Laboratory Berkeley, CA,
94720, USA  ExtreMe Matter Institute EMMI, GSI, 64291 Darmstadt, Germany  Anar
Rustamov<EMAIL_ADDRESS>GSI Helmholtzzentrum für Schwerionenforschung,
64291 Darmstadt, Germany  Joachim Stroth<EMAIL_ADDRESS>Institut für
Kernphysik, Goethe-Universität, 60438 Frankfurt am Main, Germany  GSI
Helmholtzzentrum für Schwerionenforschung, 64291 Darmstadt, Germany  Helmholtz
Research Academy Hesse for FAIR (HFHF), Campus Frankfurt , 60438 Frankfurt am
Main, Germany

###### Abstract

We generalize and extend the recently proposed method Rustamov  _et al._
(2023) to account for contributions of system size (or volume/participant)
fluctuations to the experimentally measured moments of particle multiplicity
distributions. We find that in the general case there are additional biases
which are not directly accessible to experiment. These biases are, however,
parametrically suppressed if the multiplicity of the particles of interest is
small compared to the total charged-particle multiplicity, e.g., in the case
of proton number fluctuations at top RHIC and LHC energies. They are also
small if the multiplicity distribution of charged particles per wounded
nucleon is close to the Poissonian limit, which is the case at low energy
nuclear collisions, e.g., at GSI/SIS18. We further find that mixed events are
not necessarily needed to extract the correction for volume fluctuations,
albeit it can help if event statistics is small, which is typically the case
for reconstructing the higher-order cumulants. We provide the formulas to
correct pure and mixed cumulants of particle multiplicity distributions up to
any order together with their associated biases.

††preprint: This line only printed with preprint option

##  I Introduction

One of the main goals of studying relativistic heavy-ion collisions is to
explore the structure of the QCD phase diagram. Fluctuations of observed
particles carrying quantum numbers of conserved charges, baryon number (B),
electric charge and strangeness, represent a powerful tool for this endeavor
as the cumulants of their distributions measure the derivatives of the grand-
canonical partition function, and thus the pressure ($P$), with respect to the
associated chemical potentials. For example, for a thermal system of volume
$V$ and temperature $T$, the cumulants of the net baryon number distribution,
within the Grand Canonical Ensemble (GCE), are given by Gross  _et al._ (2023)

$\kappa_{n}[B]=\frac{\partial^{n}(\ln
Z)}{\partial\,(\mu_{B}/T)^{n}}=\frac{V}{T}\frac{\partial^{n}P}{\partial\,(\mu_{B}/T)^{n}},$

where $Z$ is a GCE partition function and $\mu_{B}$ is a baryon chemical
potential. Any nontrivial structures in the equation of state such as a
possible phase transition Stephanov (2009, 2011); Bzdak  _et al._ (2020);
Asakawa  _et al._ (2009) will result in potentially large derivatives of the
pressure and thus in large values of the cumulants of conserved charges. In
addition, as cumulants are derivatives of the pressure, they are accessible
(at vanishing or small values of chemical potential) to Lattice QCD
calculations Bazavov  _et al._ (2017); Borsanyi  _et al._ (2018), which in
principle enables a direct comparison of results from ab initio QCD
calculations with experiment. For example, as pointed out in Ref. Friman  _et
al._ (2011), the measurement of higher-order cumulants close to vanishing
chemical potential may test the remnants of chiral criticality.

Measurements of fluctuations have meanwhile been carried out by many
experiments. The STAR collaboration has measured cumulants of the net-proton
number up to sixth order over the entire energy range available at RHIC
Abdallah  _et al._ (2021, 2023). The HADES experiment has measured cumulants
of proton number up to forth order at the low energy of
$\sqrt{s_{\mathrm{NN}}}=2.4{\rm\,GeV}$ Adamczewski-Musch  _et al._ (2020) and
ALICE has measured the second- and third-order net-proton number cumulants at
$\sqrt{s_{\mathrm{NN}}}=2.76$ and $5.02{\rm\,TeV}$ Acharya  _et al._ (2020,
2023).

When comparing cumulants measured in experiment with those obtained from
lattice QCD or other field theoretical calculations Fu  _et al._ (2016) one
needs to be aware of several key differences. While theoretical calculations
are typically done in the grand canonical ensemble where charges can be
exchanged with a heat bath and are only conserved on the average, in
experiment charges are explicitly conserved on event by event basis and one
has to account for global as well as local charge conservation Bleicher  _et
al._ (2000); Bzdak  _et al._ (2013); Braun-Munzinger  _et al._ (2019, 2021,
2023). Also, in experiments one usually is restricted to the measurement of
net protons whereas theory can only calculate cumulants of the net baryon
number. In the presence of many pions this difference can be corrected for
Kitazawa and Asakawa (2012). Finally, and this will be the topic of the
present paper, in experiment the size of the particle emitting system is not
constant. Even under the tightest centrality selection criteria, this gives
rise to so-called volume fluctuations Skokov  _et al._ (2013) or,
equivalently, fluctuations of the number of wounded nucleons Braun-Munzinger
_et al._ (2017). Moreover, centrality is a concept assuming a strict
correlation of event activity, i.e., charged particle multiplicity in a broad
region around mid rapidity, with the size of the system. This requires a
strict separation of the particles used to determine the fluctuation of
conserved charges from the ones used for centrality determination. The effects
due to volume fluctuations may be sizable, especially at lower energies where
the charged particle multiplicity is dominated by the primordial protons,
limiting the achievable resolution of the centrality selection. In Ref.
Rustamov  _et al._ (2023) a novel and promising method based on event mixing
has been proposed to experimentally determine and subtract the contributions
to the cumulants caused by volume fluctuations. In the present work we will
further elaborate on this topic, generalize the results, and provide the
formulas for corrections of any higher-order cumulants.

This paper is organized as follows. In the next section we define the
notation. We then present an analytical formulation of event mixing as
proposed in Ref. Rustamov  _et al._ (2023). We find that the cumulants of the
the mixed events have additional bias terms which were assumed to vanish in
the original work of (Rustamov  _et al._ , 2023), and we discuss the magnitude
of these corrections for various scenarios. Next we extend our study to
cumulants of higher order before we discuss and summarize our results.

##  II Notation

In this paper we will mostly work within the wounded-nucleon model Bialas  _et
al._ (1976) to discuss volume or participant fluctuations. We would like to
point out that this model has its limitations when applied to low collision
energies because of the moderate separation of projectile and target rapidity
and a multiplicity of created particles per wounded nucleon much smaller than
unity. However, as we shall show later, the formalism can be easily applied
also to the situation where one has generic volume fluctuations, as for
example discussed in Refs. Skokov  _et al._ (2013); Braun-Munzinger  _et al._
(2017). Let us start with the expression of the particle number cumulants
$\kappa_{j}[N]$ in the presence of wounded-nucleon fluctuations (for details
see Appendix A):

$\displaystyle\kappa_{1}[N]$ $\displaystyle=\left\langle
N_{w}\right\rangle\kappa_{1}[n]=\left\langle N_{w}\right\rangle\left\langle
n\right\rangle=\left\langle N\right\rangle$ (1) $\displaystyle\kappa_{2}[N]$
$\displaystyle=\left\langle N_{w}\right\rangle\kappa_{2}[n]+\left\langle
n\right\rangle^{2}\kappa_{2}[N_{w}]=\bar{\kappa}_{2}[N]+\left\langle
N\right\rangle^{2}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}$ (2) $\displaystyle\kappa_{3}[N]$
$\displaystyle=\left\langle N_{w}\right\rangle\kappa_{3}[n]+3\left\langle
n\right\rangle\kappa_{2}[n]\kappa_{2}[N_{w}]+\left\langle
n\right\rangle^{3}\kappa_{3}[N_{w}]=\bar{\kappa}_{3}[N]+3\left\langle
N\right\rangle\bar{\kappa}_{2}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+\left\langle
N\right\rangle^{3}\frac{\kappa_{3}[N_{w}]}{\left\langle
N_{w}\right\rangle^{3}}$ (3) $\displaystyle\kappa_{4}[N]$
$\displaystyle=\left\langle N_{w}\right\rangle\kappa_{4}[n]+4\left\langle
n\right\rangle\kappa_{3}[n]\kappa_{2}[N_{w}]+3\kappa_{2}^{2}[n]\kappa_{2}[N_{w}]+6\left\langle
n\right\rangle^{2}\kappa_{2}[n]\kappa_{3}[N_{w}]+\left\langle
n\right\rangle^{4}\kappa_{4}[N_{w}]$
$\displaystyle=\bar{\kappa}_{4}[N]+4\left\langle
N\right\rangle\bar{\kappa}_{3}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+3\bar{\kappa}_{2}^{2}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+6\left\langle
N\right\rangle^{2}\bar{\kappa}_{2}[N]\frac{\kappa_{3}[N_{w}]}{\left\langle
N_{w}\right\rangle^{3}}+\left\langle
N\right\rangle^{4}\frac{\kappa_{4}[N_{w}]}{\left\langle
N_{w}\right\rangle^{4}}$ (4)

Here $N$ refers to the particles of interest, say protons, and $n$ to the
number of these particles arising from one wounded nucleon; thus $\left\langle
n\right\rangle$ is the average number of particles per wounded nucleon. The
cumulants of the wounded-nucleon distribution are denoted by
$\kappa_{j}[N_{w}]$ while the cumulants for the distribution of particles
stemming from one wounded nucleon are $\kappa_{j}[n]$. The corresponding
relations for cumulants of any order can be obtained with the provided
software package tes .

The cumulants of interest are those at a fixed number of wounded nucleons.
They reflect the true density fluctuations in a system at constant volume. We
denote these cumulants for a system with fixed, i.e. non-fluctuating, number
of $\left\langle N_{w}\right\rangle$ wounded nucleons as

$\bar{\kappa}_{j}[N]=\left\langle N_{w}\right\rangle\kappa_{j}[n],$

Below we will also deal with factorial cumulants, which we shall denote by
$C_{j}$. Factorial cumulants, which measure the deviation from Poisson
statistics, tell us about the true correlations in the system. As discussed in
the Appendix B, they are linear combinations of the regular cumulants. For the
first four orders we have

$\displaystyle C_{1}[N]$ $\displaystyle=\kappa_{1}[N]=\left\langle
N\right\rangle,$ $\displaystyle C_{2}[N]$
$\displaystyle=-\kappa_{1}[N]+\kappa_{2}[N],$ $\displaystyle C_{3}[N]$
$\displaystyle=2\kappa_{1}[N]-3\kappa_{2}[N]+\kappa_{3}[N],$ $\displaystyle
C_{4}[N]$
$\displaystyle=-6\kappa_{1}[N]+11\kappa_{2}[N]-6\kappa_{3}[N]+\kappa_{4}[N].$

The expressions for the particle number factorial cumulants are similar to
Eqs. 1\- 4

$\displaystyle C_{1}[N]$ $\displaystyle=\left\langle N_{w}\right\rangle
C_{1}[n]=\left\langle N_{w}\right\rangle\left\langle
n\right\rangle=\left\langle N\right\rangle,$ (5) $\displaystyle C_{2}[N]$
$\displaystyle=\bar{C}_{2}[N]+\left\langle
N\right\rangle^{2}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}},$ (6) $\displaystyle C_{3}[N]$
$\displaystyle=\bar{C}_{3}[N]+3\left\langle
N\right\rangle\bar{C}_{2}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+\left\langle
N\right\rangle^{3}\frac{\kappa_{3}[N_{w}]}{\left\langle
N_{w}\right\rangle^{3}},$ (7) $\displaystyle C_{4}[N]$
$\displaystyle=\bar{C}_{4}[N]+4\left\langle
N\right\rangle\bar{C}_{3}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+3\bar{C}_{2}^{2}[N]\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+6\left\langle
N\right\rangle^{2}\bar{C}_{2}[N]\frac{\kappa_{3}[N_{w}]}{\left\langle
N_{w}\right\rangle^{3}}+\left\langle
N\right\rangle^{4}\frac{\kappa_{4}[N_{w}]}{\left\langle
N_{w}\right\rangle^{4}}.$ (8)

Similar to the cumulants, we denote by

$\bar{C}_{k}[N]=\left\langle N_{w}\right\rangle C_{k}[n]$

the factorial cumulants for a system at constant volume or number of wounded
nucleons, $\left\langle N_{w}\right\rangle$.

##  III Mixed events

In Ref. (Rustamov  _et al._ , 2023) a mixed event is constructed such that it
has the same total multiplicity as a given real event but each particle
(track) is drawn from a different event, so that, by construction, the mixed
events follow the same total multiplicity distribution as the original events.
This is done in order to preserve volume flucutations as in real events. Since
each particle (track) is chosen randomly from a random event, the distribution
of particle species will follow a multinomial distribution with the Bernoulli
probabilities $p_{\rm i}=\left\langle N_{i}\right\rangle/\left\langle
M\right\rangle$ for particles of type $i$. Here $\left\langle
N_{i}\right\rangle$ denotes the mean number of particles of type $i$ and
$\left\langle M\right\rangle$ the mean total multiplicity. Hence, the
probability to find $A$ particles (successes) of type A and $B$ particles of
type B is given by the trinomial probability $B_{3}(A,B,M;p_{\rm A},p_{\rm
B})$ and so on. Here $M$ denotes the multiplicity of the event under
consideration. Thus the distribution, $P_{mix}\left(A,B\right)$, of particles
of species A and B in the mixed events is obtained by folding the multiplicity
distribution $P_{M}(M)$ with a trinomial (in general multinomial)
distribution:

$P_{mix}\left(A,B\right)=\sum_{M}B_{3}(A,B,M;p_{\rm A},p_{\rm B})P_{M}(M)$

with

$p_{\rm A}=\frac{\left\langle A\right\rangle}{\left\langle
M\right\rangle},\;\;\;\;\;p_{\rm B}=\frac{\left\langle
B\right\rangle}{\left\langle M\right\rangle}.$

and

$\displaystyle B_{3}(A,B,M;p_{\rm A},p_{\rm
B})=\frac{M!}{A!B!(M-A-B)!}\,p_{\rm A}^{A}\,p_{\rm B}^{B}\,(1-p_{\rm A}-p_{\rm
B})^{M-A-B}$ (9)

The factorial-cumulant generating function for this distribution is

$\displaystyle g_{F,mx}\left(z_{\rm A},z_{\rm B}\right)$
$\displaystyle=\ln\left[\sum_{A,B}P_{mx}\left(A,B\right)(z_{\rm A})^{A}(z_{\rm
B})^{B}\right]$ $\displaystyle=\ln\left[\sum_{M}\left[h_{3}\left(z_{\rm
A},z_{\rm B}\right)\right]^{M}P_{M}(M)\right]$
$\displaystyle=G_{F,M}\left(h_{3}\left(z_{\rm A},z_{\rm B}\right)\right)$ (10)

where

$h_{3}\left(z_{\rm A},z_{\rm B}\right)=\sum_{A,B}B_{3}(A,B;M=1;p_{\rm
A},p_{\rm B})(z_{\rm A})^{A}(z_{\rm B})^{B}=\left(1-p_{\rm A}-p_{\rm B}+p_{\rm
A}z_{\rm A}+p_{\rm B}z_{\rm B}\right)$

is the factorial-moment generating function for the trinomial distribution
with one trial $\left(M=1\right)$, and $G_{F,M}\left(z\right)$ is the
factorial-cumulant generating function for the multiplicity distribution,
$P_{M}(M)$ (see Eq.(52)). The factorial cumulants are then obtained via

$C_{i,j}^{mix}[mix]=\left.\frac{\partial^{i}\partial^{j}}{\partial(z_{\rm
A})^{i}\partial(z_{\rm B})^{j}}G_{F,M}\left(h_{3}\left(z_{\rm A},z_{\rm
B}\right)\right)\right|_{z_{\rm A}=z_{\rm B}=0}=p_{\rm A}^{i}p_{\rm
B}^{j}C_{i+j}\left[M\right]$

with $C_{k}[M]$ being the $k^{th}$-order factorial cumulant. Using the
expression for the factorial cumulants of the multiplicity distribution
derived in Appendix C, Eq. 63, we get within the wounded-nucleon model

$\displaystyle C_{1}^{mix}[A]=\kappa_{1}^{mix}[A]$ $\displaystyle=p_{\rm
A}\left\langle N_{w}\right\rangle\left\langle m\right\rangle$ $\displaystyle
C_{2}^{mix}[A]$ $\displaystyle=p_{\rm A}^{2}C_{2}[M]=p_{\rm
A}^{2}\left[\kappa_{2}\left[N_{w}\right]\left\langle
m\right\rangle^{2}+\left\langle N_{w}\right\rangle C_{2}[m]\right]$
$\displaystyle C_{1,1}^{mix}\left[A,B\right]$ $\displaystyle=p_{\rm A}p_{\rm
B}C_{2}[M]=p_{\rm A}p_{\rm B}\left[\kappa_{2}\left[N_{w}\right]\left\langle
m\right\rangle^{2}+\left\langle N_{w}\right\rangle C_{2}[m]\right].$ (11)

For the corresponding cumulants up to second order we get accordingly

$\displaystyle\kappa_{1}^{mix}[A]=C_{1}^{mix}[A]$ $\displaystyle=p_{\rm
A}\left\langle N_{w}\right\rangle\left\langle m\right\rangle$
$\displaystyle\kappa_{2}^{mix}[A]=C_{2}^{mix}[A]+C_{1}^{mix}[A]$
$\displaystyle=p_{\rm A}^{2}\left[\kappa_{2}\left[N_{w}\right]\left\langle
m\right\rangle^{2}+\left\langle
N_{w}\right\rangle\left(\kappa_{2}[m]-\kappa_{1}[m]\right)\right]+p_{\rm
A}\left\langle N_{w}\right\rangle\left\langle m\right\rangle$
$\displaystyle=p_{\rm A}^{2}\left[\kappa_{2}\left[N_{w}\right]\left\langle
m\right\rangle^{2}+\left\langle
N_{w}\right\rangle\left(\kappa_{2}[m]-\left\langle
m\right\rangle\right)\right]+p_{\rm A}\left\langle
N_{w}\right\rangle\left\langle m\right\rangle$ $\displaystyle
cov^{mix}\left[A,B\right]=C_{1,1}^{mix}\left[A,B\right]$ $\displaystyle=p_{\rm
A}p_{\rm B}\left[\kappa_{2}\left[N_{w}\right]\left\langle
m\right\rangle^{2}+\left\langle
N_{w}\right\rangle\left(\kappa_{2}[m]-\left\langle
m\right\rangle\right)\right].$ (12)

With $\left\langle m\right\rangle$ denoting the mean number of total particles
emitted by a wounded nucleon, we get $\left\langle a\right\rangle=p_{\rm
A}\left\langle m\right\rangle$ and $\left\langle b\right\rangle=p_{\rm
B}\left\langle m\right\rangle$ for the mean number of particles per wounded
nucleon of type A and B, respectively, and recover the results of Ref.Rustamov
_et al._ (2023). For that we have to assume that the multiplicity distribution
per wounded nucleon is Poissonian, i.e. that
$C_{2}[m]=\kappa_{2}[m]-\left\langle m\right\rangle=0$. This has been an
implicit assumption in Ref. Rustamov  _et al._ (2023), which however is not
valid in general as we shall discuss below.

The main benefit of the event mixing is to be able to relate the factorial
cumulants of the various multiplicity distributions, as can be seen from
Eq.11. All that enters is the second-order factorial cumulant, $C_{2}[M]$. The
binomial probabilities, $p_{\rm A}$ and $p_{\rm B}$, are in the sense trivial
as they can be determined without any mixed events. Thus we may express the
fluctuations of the wounded nucleons in terms of the factorial cumulant of the
track multiplicity distribution

$\left\langle N\right\rangle^{2}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}=\frac{\left\langle N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\left(C_{2}[M]-\left\langle N_{w}\right\rangle
C_{2}[m]\right)=\frac{\left\langle N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\left(C_{2}[M]-\bar{C}_{2}[M]\right),$ (13)

where $\bar{C}_{2}[M]=\left\langle N_{w}\right\rangle C_{2}[m]$ is the second-
order factorial cumulant for a system of $\left\langle N_{w}\right\rangle$
wounded nucleons _without_ wounded nucleon fluctuations and $N$ stands now for
the multiplicity of the particles of interest, i.e., either $A$ or $B$. While
the factorial cumulant of the multiplicity distribution, $C_{2}[M]$, is
accessible to experiment, that of a non-fluctuating system, $\bar{C}_{2}[M]$,
is not. Let us, therefore define a bias term, $\Delta_{2}$, as

$\Delta_{2}\equiv\frac{\left\langle N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\bar{C}_{2}[M].$ (14)

In case of a Poissonian multiplicity distribution for one wounded nucleon the
bias term vanishes, i.e., $\Delta_{2}=0$, since $\bar{C}_{2}[M]=\left\langle
N_{w}\right\rangle C_{2}[m]=0$ in this case, and we recover the results of
Ref. Rustamov  _et al._ (2023). Let us furthermore define the corrected
cumulant, $\kappa_{2}^{corr}[N],$ which is based on measurable quantities only

$\kappa_{2}^{corr}[N]=\kappa_{2}[N]-\frac{\left\langle
N\right\rangle^{2}}{\left\langle M\right\rangle^{2}}C_{2}[M].$ (15)

Following Eq. 2 and using Eq. 6, the cumulant of the system _without_ wounded
nucleon fluctuations, $\bar{\kappa}_{2}[N]$, is given by

$\displaystyle\bar{\kappa}_{2}[N]=\kappa_{2}[N]-\left\langle
N\right\rangle^{2}\frac{\kappa_{2}[N_{W}]}{\left\langle
N_{W}\right\rangle^{2}}$ $\displaystyle=\kappa_{2}^{corr}[N]+\Delta_{2}.$ (16)

The bias, $\Delta_{2}$, while not directly measurable, may be constrained by a
fit to the track multiplicity distribution within the wounded-nucleon model
Bialas  _et al._ (1976), as it is commonly done Adamczewski-Musch  _et al._
(2018); Abdallah  _et al._ (2023); Abelev  _et al._ (2013a). In addition, we
note that for protons at very high collision energies we have $\left\langle
N_{p}\right\rangle\ll\left\langle M\right\rangle$ so that $\Delta_{2}$ is
suppressed parametrically. This behavior can indeed be illustrated with
simulations as presented in Sec. V. Since cumulants scale with the system
size, or in our case with the number of wounded nucleons, $\left\langle
N_{w}\right\rangle$, it is instructive to scale the (factorial) cumulants with
the mean number of particles

$\frac{\bar{\kappa}_{2}[N]}{\left\langle
N\right\rangle}=\frac{\kappa_{2}[N]}{\left\langle
N\right\rangle}-\frac{\left\langle N\right\rangle}{\left\langle
M\right\rangle}\left(\frac{C_{2}[M]}{\left\langle
N\right\rangle}-\frac{\bar{C}_{2}[M]}{\left\langle
N\right\rangle}\right)=\frac{\kappa_{2}^{corr}[N]}{\left\langle
N\right\rangle}+\frac{\Delta_{2}}{\left\langle N\right\rangle}$ (17)

The scaled bias is then given by

$\frac{\Delta_{2}}{\left\langle N\right\rangle}=\frac{\left\langle
N\right\rangle}{\left\langle M\right\rangle}\bar{c}_{2}[M].$ (18)

Typically, the scaled cumulants are of order unity,
$\kappa_{j}[N]/\left\langle N\right\rangle\sim\mathcal{O}(1)$. In addition,
the scaled factorial cumulants, $C_{k}[N]/\left\langle N\right\rangle$, are
expected to depend only weakly on the multiplicity. Therefore, the scaled bias
should be much smaller than one, $\Delta_{2}/\left\langle N\right\rangle\ll
1$, for the volume correction to be reliable.

Finally, one may express the bias term $\Delta_{2}$ also in terms of cumulants
by using the relation between cumulants and factorial cumulants (see Appendix
B), $C_{2}[M]=\kappa_{2}[M]-\left\langle M\right\rangle$ and and so forth.
This gives,

$\displaystyle\Delta_{2}$ $\displaystyle=\frac{\left\langle
N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\left(\bar{\kappa}_{2}[M]-\left\langle
M\right\rangle\right)$ (19)

A note of caution may be useful in this context. One might be inclined to
express the fluctuations of the wounded nucleon directly using the cumulants
of the multiplicity distribution, in which case one would get

$\left\langle N\right\rangle^{2}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}=\frac{\left\langle N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\left(\kappa_{2}[M]-\left\langle
N_{w}\right\rangle\kappa_{2}[m]\right).$

And since $\left\langle N_{W}\right\rangle\kappa_{2}[m]=\bar{\kappa}_{2}[M]$
is not directly accessible to experiment, one may further assign the bias to
be $\Delta_{2}=\frac{\left\langle N\right\rangle^{2}}{\left\langle
M\right\rangle^{2}}\bar{\kappa}_{2}[M]$. This, however, would considerably
overestimate its true value, Eq. 19, as cumulants always contain a “trivial”
component proportional to the number of particles, which in principle is
measurable.

##  IV Higher-Order Results

Let us now discuss the corrections for volume fluctuations up to fourth order.
Given the discussion in the previous section the strategy is straightforward.
First we express the fluctuations of the wounded nucleons in terms of
factorial cumulants of the multiplicity distribution. Then we identify the
parts which are experimentally accessible and those which are not. The latter
will be the bias while the former will be subtracted from the expression for
the cumulants in order to remove most of the effect of volume fluctuations.
The terms involving cumulants of the wounded-nucleon distribution as they
appear in the expressions for the cumulants as $\kappa_{j}[N_{w}]/\left\langle
N_{w}\right\rangle^{j}$, see Eqs. (6-8) are:

$\displaystyle\frac{\kappa_{2}[N_{w}]}{\left\langle N_{w}\right\rangle^{2}}$
$\displaystyle=\frac{C_{2}[M]-\bar{C}_{2}[M]}{\langle M\rangle^{2}}$ (20)
$\displaystyle\frac{\kappa_{3}[N_{w}]}{\left\langle N_{w}\right\rangle^{3}}$
$\displaystyle=-3\frac{\bar{C}_{2}[M]}{\langle
M\rangle^{2}}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+\frac{C_{3}[M]-\bar{C}_{3}[M]}{\langle M\rangle^{3}}$
(21) $\displaystyle\frac{\kappa_{4}[N_{w}]}{\left\langle
N_{w}\right\rangle^{4}}$ $\displaystyle=-6\frac{\bar{C}_{2}[M]}{\langle
M\rangle^{2}}\frac{\kappa_{3}[N_{w}]}{\left\langle
N_{w}\right\rangle^{3}}-\frac{4\bar{C}_{3}[M]\langle
M\rangle+3\bar{C}_{2}[M]{}^{2}}{\langle
M\rangle^{4}}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}+\frac{C_{4}[M]-\bar{C}_{4}[M]}{\langle M\rangle^{4}}$
(22)

We note that binomial efficiency corrections do not affect the results as
both, numerators and denominators of the right hand side of the above
expressions, scale with the same power of the efficiency.

Inserting these expressions into Eqs. (2-4) for the cumulants $\kappa_{j}[N]$,
we can solve for the cumulants of the system with fixed number of wounded
nucleons, namely the $\bar{\kappa}_{j}[N]$. The results are given in the
following general form

$\bar{\kappa}_{j}[N]=\kappa_{j}^{corr}[N]+\Delta_{j}[N]$ (23)

with $\bar{\kappa}_{j}[N]$ the cumulant of order $j$ for a system with fixed
$N_{w}$ nucleons, $\kappa_{j}^{corr}[N]$ the cumulant including the measurable
corrections for volume fluctuations, and $\Delta_{j}$ the corresponding bias
due to quantities that are not measurable. The second-order result we already
derived in Sec. A, Eqs. (15) and (14), namely

$\displaystyle\kappa_{2}^{corr}[N]$
$\displaystyle=\kappa_{2}[N]-\frac{\left\langle
N\right\rangle^{2}}{\left\langle M\right\rangle^{2}}C_{2}[M]$
$\displaystyle\Delta_{2}$ $\displaystyle=\frac{\left\langle
N\right\rangle^{2}}{\left\langle M\right\rangle^{2}}\bar{C}_{2}[M].$ (24)

For the third order we have

$\displaystyle\kappa_{3}^{corr}[N]$
$\displaystyle=\kappa_{3}[N]-\frac{3C_{2}[M]\kappa_{2}[N]\langle
N\rangle}{\langle M\rangle^{2}}+\frac{3C_{2}[M]{}^{2}\langle
N\rangle^{3}}{\langle M\rangle^{4}}-\frac{C_{3}[M]\langle
N\rangle^{3}}{\langle M\rangle^{3}}$ $\displaystyle\Delta_{3}$
$\displaystyle=\bar{C}_{2}[M]\left(\frac{3\kappa_{2}[N]\langle
N\rangle}{\langle M\rangle^{2}}-\frac{3C_{2}[M]\langle N\rangle^{3}}{\langle
M\rangle^{4}}\right)+\frac{\bar{C}_{3}[M]\langle N\rangle^{3}}{\langle
M\rangle^{3}}.$ (25)

And the fourth order result reads

$\displaystyle\kappa_{4}^{corr}[N]$
$\displaystyle=\kappa_{4}[N]-\left(\frac{6\kappa_{2}[N]\langle
N\rangle^{2}\left(C_{3}[M]\langle M\rangle-3C_{2}[M]{}^{2}\right)}{\langle
M\rangle^{4}}+\frac{4C_{2}[M]\kappa_{3}[N]\langle N\rangle}{\langle
M\rangle^{2}}\right.$
$\displaystyle\left.+\frac{3C_{2}[M]\kappa_{2}(N){}^{2}}{\langle
M\rangle^{2}}+\frac{\langle N\rangle^{4}\left(-10C_{3}[M]C_{2}[M]\langle
M\rangle+15C_{2}[M]{}^{3}\right)}{\langle M\rangle^{6}}+\frac{C_{4}[M]\langle
N\rangle^{4}}{\left\langle M\right\rangle^{4}}\right)$
$\displaystyle\Delta_{4}$
$\displaystyle=\bar{C}_{2}[M]\left(-\frac{18C_{2}[M]\kappa_{2}[N]\langle
N\rangle^{2}}{\langle M\rangle^{4}}+\frac{15C_{2}[M]{}^{2}\langle
N\rangle^{4}}{\langle M\rangle^{6}}-\frac{4C_{3}[M]\langle
N\rangle^{4}}{\langle M\rangle^{5}}+\frac{4\kappa_{3}[N]\langle
N\rangle}{\langle M\rangle^{2}}+\frac{3\kappa_{2}[N]{}^{2}}{\langle
M\rangle^{2}}\right)$
$\displaystyle+\bar{C}_{3}[M]\left(\frac{6\kappa_{2}[N]\langle
N\rangle^{2}}{\langle M\rangle^{3}}-\frac{6C_{2}[M]\langle
N\rangle^{4}}{\langle M\rangle^{5}}\right)+\frac{\bar{C}_{4}[M]\langle
N\rangle^{4}}{\langle M\rangle^{4}}.$ (26)

The corresponding relations for correction and bias terms of any order can be
obtained with the provided software package tes .

Equivalent expressions for the factorial cumulants, $C_{n}[N]$, and their
related biases, $\Delta_{n,F}$, may then be obtained by using the relation
between factorial cumulants and regular cumulants Eq. 56:

$\displaystyle C_{n}^{corr}[N]$
$\displaystyle=\sum_{j=1}^{n}B_{n,j}\left(1,-1,2,\ldots,(-1)^{j-1}(n-j+1)!\right)\,\kappa^{corr}_{j}[N]$
$\displaystyle\Delta_{n,F}$
$\displaystyle=\sum_{j=1}^{n}B_{n,j}\left(1,-1,2,\ldots,(-1)^{j-1}(n-j+1)!\right)\,\Delta_{j}[N]$

The results for the corrected factorial cumulants, $C_{k}^{corr}$ and
$c_{k}^{corr}$, and the associated biases, $\Delta_{k,F}$ and $\delta_{k,F}$,
are given in Appendix E.

##  V Simulations

Experimental data are usually analyzed in centrality percentiles, i.e. event
classes corresponding to the $n\%$ most central collisions, by introducing
selection criteria on e.g. the energy deposited in a forward detector system
covering typically the projectile (target) spectator region or the
multiplicity of charged particles emitted from the mid-rapidity region, with
an acceptance reaching close to the projectile (target) rapidity regions in
case of low beam energies Abelev  _et al._ (2013a); Adamczewski-Musch  _et
al._ (2018). For the latter, care must be taken to ensure that the evaluated
particles are not simultaneously used to determine the critical fluctuations
Adam  _et al._ (2021).

Figure 1: NBD adjusted to HADES (a), STAR (b) and ALICE (c) charged-particle
multiplicity distributions (blue line), shown together with a Poisson
distribution for the HADES and STAR data (red lines), the parameters used are
listed in Table I.

Figure 2: Distribution of wounded nucleons in Au+Au collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 2.4 GeV for four selected centrality classes, as
obtained from the Glauber Monte Carlo simulations.

The respective distributions, like e.g., the forward energy deposit or the
charged particle multiplicity, are commonly modelled using the Glauber Monte
Carlo Model Loizides  _et al._ (2015). The model provides event by event and
for a given impact parameter the number of projectile/target nucleons which
are “wounded” and responsible for the event activity (multiplicity), and
those, which proceed nearly undisturbed into the phase space region covered by
the forward detectors. To determine centrality using charged particle
multiplicity the respective distribution is generally modelled assuming that
particles are “produced” independently from distinct sources following a
negative binomial distribution (NBD). Its probability mass function is defined
as

$P(n;\mu,k)=\frac{\Gamma(n+k)}{\Gamma(n+1)\Gamma(k)}\left(\frac{\mu}{k}\right)^{n}\left(\frac{\mu}{k}+1\right)^{-(n+k)},$
(27)

where $\mu$ denotes the mean of the NBD, while the combination of $\mu$ and
$k$ determines its higher-order cumulants

$\kappa_{n}^{NBD}=\left.\frac{\partial^{n}\ln M(t)}{\partial
t^{n}}\right\rvert_{t=0},$ (28)

with

$M(t)=\sum_{n=0}^{\infty}e^{tn}P(n;\mu,k)=\left(\frac{k}{k+\mu(1-e^{t})}\right)^{k}$
(29)

being the moment-generating function of the NBD. The first four cumulants read

$\displaystyle\kappa_{1}^{\text{NBD}}=\mu,$ (30)
$\displaystyle\kappa_{2}^{\text{NBD}}=\frac{\mu(k+\mu)}{k},$ (31)
$\displaystyle\kappa_{3}^{\text{NBD}}=\frac{\mu(k+\mu)(k+2\mu)}{k^{2}},$ (32)
$\displaystyle\kappa_{4}^{\text{NBD}}=\frac{\mu(k+\mu)(k^{2}+6k\mu+6\mu^{2})}{k^{3}},$
(33)

The parameters of the NBD are fixed in each experiment by the fitting
procedure. In a first step the number of particle-emitting sources $n_{s}$ is
determined according to Abelev  _et al._ (2013a); Abdallah  _et al._ (2023)

$n_{\text{s}}=fN_{w}+(1-f)N_{coll},$ (34)

where $N_{W}$ and $N_{coll}$ are the numbers of wounded nucleons and binary
collisions, respectively. Sampling impact parameters according to
$d\sigma=b\,db$ a list of number of sourses,$n_{\text{s}}^{i}$, is generated,
with $i\in[1,\cdots,N_{\text{Event}}]$. Then, for each event $i$, NBD is
sampled $n_{\text{s}}^{i}$ times and the parameters of the NBD, $\mu$, $k$ and
$f$, are adjusted such that the obtained multiplicity distribution agrees with
the corresponding experimental one. The mixing parameter $f$ is introduced to
improve the description by accounting also particles produced in hard (prompt)
processes.

Figure 3: Simulation of particle production within the model of independent
sources. The circles indicate the individual sources sampled according to Eq.
34. $n_{1}$ and $n_{2}$ show particle species used to sample emission from a
single source. The distributions per single source can be chosen arbitrarily.
See models A and B discussed below.

Figure 4: Left panel: Mean number of simulated protons used in Model A as a
function of centrality. Right panel: Second-order cumulants of protons in
Model A including volume fluctuations (black stars), corrected with Eq. 24
(blue stars) and reconstructed with Eq. 23 (open red stars). The results are
normalized to $\bar{\kappa}_{2}[p]_{sim}$, corresponding to the second-order
cumulants of protons in the absence of volume fluctuations.

Figure 5: Left panel: Third order cumulants of protons in Model A including
volume fluctuations (black stars), corrected with Eq. 25 (blue stars) and
reconstructed with Eq. 23 (open red stars). The results are normalized to
$\bar{\kappa}_{3}[p]_{\text{sim}}$, corresponding to the second-order
cumulants of protons in the absence of volume fluctuations. Right panel:
Fourth order cumulants of protons in Model A including volume fluctuations
(black stars), corrected with Eq. 26 (blue stars) and reconstructed with Eq.
23 (open red stars). The results are normalized to
$\bar{\kappa}_{4}[p]_{\text{sim}}$, corresponding to the fourth order
cumulants of protons in the absence of volume fluctuations.  Figure 6: The
normalised bias terms in Model A.

Panel (a) of Fig. 1 represents the NBD distribution as observed by the HADES
experiment for Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=2.4$ GeV, with
parameters $\mu=0.24$, $k=20.34$, and $f=1$ taken from Adamczewski-Musch  _et
al._ (2018). For comparison, we also present a Poisson distribution with the
same mean, $\mu$. Similar distributions from the STAR Abdallah  _et al._
(2023) and ALICE Abelev  _et al._ (2013a) experiments are presented in panels
(b) and (c). Figure 1 shows that at the HADES energy the fitted NBDs are very
close to Poisson distributions. Quantitatively this can be seen by evaluating
the cumulants of the HADES NBD ($\mu=0.24$, $k=20.34$)

$\displaystyle\kappa_{1}^{\text{NBD}}(\text{HADES})=0.24$ (35)
$\displaystyle\kappa_{2}^{\text{NBD}}(\text{HADES})=0.2428,$ (36)
$\displaystyle\kappa_{3}^{\text{NBD}}(\text{HADES})=0.2486,$ (37)
$\displaystyle\kappa_{4}^{\text{NBD}}(\text{HADES})=0.2602.$ (38)

For a Poisson distribution all cumulants are equal to its mean and the HADES
data are indeed close to fulfilling this condition. The statement, to a lesser
extent, is also valid for the STAR Au+Au data at 3 GeV (see Fig. 1). The
corresponding ALICE distribution, however, is much wider compared to the
Poisson baseline, but the ALICE NBD is obtained for very different acceptance
than that used for the cumulant analysis. In Table 1 we also provide the NBD
parameters as obtained by the STAR and ALICE collaborations for Au+Au and
Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=3$ GeV and 2.76 TeV, respectively.

experiments | $\mu$ | $k$ | $f$
---|---|---|---
HADES | 0.24 | 20.34 | 1
STAR | 0.31 | 5.66 | 0.94
ALICE | 29.3 | 1.6 | 0.8
Table 1:  NBD parameters as extracted from Glauber fits to particle
distributions observed in different experiments. For simulations, the
distributions should be folded within the experimental acceptance in which the
cumulants are measured.

In the following we test the proposed method using two different simulations
referred to as Model A and Model B. While the sampling of wounded nucleons is
the same for both models, in Model A we sample different particle species
independently while in Model B we introduce correlations between pions and
protons via cluster production and decay. We will concentrate on the HADES
data. For both methods, we use the Glauber model to extract the distributions
of wounded nucleons corresponding to four different centrality classes in
Au+Au collisions at $\sqrt{s_{\mathrm{NN}}}=2.4$ GeV. They are presented in
Figure 2.

In both models particles are produced from independent sources (cf. Eq. 34).
For the HADES data the extracted number of binary collision is zero ($f=1$,
see Table 1 ), the number of sources per event are sampled exclusively from
the wounded nucleon distributions presented in Fig. 2. The simulation process
is schematically illustrated in Fig. 3.

Figure 7: Number of wounded nucleons (left panel) and binary collisions (right
panel) as generated with a Glauber Monte Carlo simulation using input from the
ALICE experiment Abelev  _et al._ (2013a). The selection corresponds to the 5%
most central Pb-Pb collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV.

##  VI Model A

In model A we first generate the charged-particle multiplicity for individual
events using the NDB distribution as extracted by experimental measurements.
In doing so we sample the NBD distribution $n_{s}$ times. Different particle
species are then taken as fractions of the total number of charged particles.
For example, from a randomly sampled NBD distribution a respective fraction is
assigned to protons. From the remaining charged particles another fraction is
assigned to positively charged pions and the rest is taken as negatively
charged pions. These fractions are chosen such that the overall probability of
having protons, positively and negatively charged pions correspond to 75$\%$,
9$\%$ and 16 $\%$ of all charged particles, respectively, based on the actual
HADES measurement in Au+Au collisions (see Harabasz  _et al._ (2023) and
references therein). In addition, we account for acceptance effects, because
the NDB parameters are obtained in different acceptance than that used for the
fluctuation analyses. We therefore fold the entire NDB distribution with a
binomial distribution such that the experimentally measured mean
multiplicities of particles in the acceptance used for fluctuation analysis
are reproduced. Volume fluctuations are naturally accounted for as for each
event the number of sources $n_{s}$ are randomly sampled from the
corresponding distributions.

The simulated mean numbers of protons are shown in the left panel of Fig. 4
for the four centrality classes. In the right panel of Fig. 4 the
reconstructed second-order cumulants of protons are presented, normalized to
the expected true cumulant, $\bar{\kappa}_{2}[p]_{sim}$. The black stars
represent the results which include volume fluctuations. The values
$\kappa_{2}^{corr}$ as calculated using Eq. 15 are shown with blue stars,
while the red stars correspond to $\bar{\kappa}_{2}[p]$ =
$\kappa_{2}^{corr}[p]$ \+ $\Delta_{2}[p]$. The results for the third and
fourth order cumulants are shown in Fig. 5. The corresponding normalized
biases $\Delta_{n}/\bar{\kappa}_{n}$ are presented in Fig. 6. We find the
normalized biases to be very small, of the order of a few present, so that the
corrected cumulants, $\kappa_{n}^{corr}$ are very close to their expected true
values, $\bar{\kappa}_{n}$. As already discussed, this is to be expected since
the multiplicity distribution per wounded nucleon in the present Model is
close to Poisson.

###  VI.1 High energy limit

In this section we apply the method to high energy collisions using ALICE
data. In doing so we first run Glauber Monte Carlo simulations for Pb-Pb
collisions at $\sqrt{s_{\mathrm{NN}}}=2.76$ TeV. The input parameters are
taken from Ref. Abelev  _et al._ (2013a). Following the ALICE experiment
Abelev  _et al._ (2013a), different centrality classes are chosen by
introducing sharp cuts on the charged-particle distributions. The
distributions of wounded nucleons and binary collisions corresponding to the
5$\%$ most central collisions are presented in Fig. 7 Braun-Munzinger  _et
al._ (2017). The reconstructed mean number of wounded nucleons and binary
collisions, corresponding to the 5$\%$ most central collisions are $\langle
N_{W}\rangle\approx 382$ and $\langle N_{coll}\rangle\approx 1685$,
respectively, consistent with the numbers given in Abelev  _et al._ (2013a).
With these numbers one can estimate a mean number of particle emitting
sources, yielding $\langle n_{s}\rangle\approx 642$ (see Eq. 34). The
corresponding mean number of charged particles can be estimated as $\langle
N_{ch}\rangle=\langle n_{s}\rangle\times\mu_{NBD}\approx$ 18811\. On the other
hand, the total number of charged particles measured inside the ALICE
acceptance is about 1601 Abelev  _et al._ (2013b). We thus folded the ALICE
NBD distribution with a binomial with the acceptance factor of
$\epsilon=1601/18811\approx 8.5\%$ to obtain the distribution within the
experimental acceptance.111We note that binomial folding of the NDB
distribution, with the acceptance factor $\epsilon$, changes only the
parameter $\mu$ of the original NDB distribution
($\mu\rightarrow\epsilon\mu$), while the parameter $k$ stays unchanged. The so
obtained NDB distribution from ALICE is presented in Fig. 8. We further note
that only the acceptance in rapidity is considered. Fluctuation analyses are
performed within a finite momentum range. Inclusion of the latter will further
reduce the discrepancy between NBD and the corresponding Poisson distribution
shown with the red histogram in Fig. 8.

Figure 8: NDB distributions from ALICE as measured by fitting the signal
amplitude in VZEROs (left panel) and after folding within the experimental
acceptance 8.5$\%$ (see the text for details). The red histogram corresponds
to a Poisson distribution with the same mean, $\mu$ as the folded NBD
distribution. The original NBD distribution from ALICE is presented in the
panel (c) of Fig. 1.

Finally using the NBD distribution presented in Fig. 8, and measured proton
number, $\langle N_{p}\rangle\approx$ 35 Abelev  _et al._ (2013b), we
estimated $\Delta_{2}[p]\approx$ 1.2. This corresponds to a bias of
$\Delta_{2}[p]/\bar{\kappa}_{2}[p]\approx 3.3\%$.

##  VII Model B

In model B we introduce correlations between charged particles, specifically
pions and protons by generating clusters (or rather resonances). Especially
for HADES energies most of the observed pions are believed to originate from
decays of Delta resonances. Therefore, the effect of such decay correlations,
while no treated quantitatively here, needs to be taken into account.
Specifically, this is done by generating clusters of particles (e.g.,
resonances) from each source and letting them decay into two different
particle species. Moreover, the clusters are generated from a Poisson
distribution. In addition we produce independent particles from each source as
well, sampled also from a Poisson distribution. Schematic illustration of the
model for a single source is given in Fig. 9. Fluctuations of sources are
introduced like in the model A. The input parameters for model B are given in
Table 2.

particles | mean numbers per source
---|---
clusters | 0.03
independent protons | 0.23
other particles | 0.21
Table 2:  Parameters for model B are mean numbers of different particles
species per source (see Harabasz  _et al._ (2023) and references therein). In
addition, each cluster decays into one proton and one pion. Numbers of
clusters, independent protons and other particles are sampled from independent
Poisson distributions.

Figure 9: Left panel: Schematic illustration of model B. From a single source
two clusters and three protons are produced. The clusters further decay to
protons and pions. In simulations both, clusters and independent protons are
sampled from Poisson distributions. The corresponding mean values of clusters
and independent protons are provided in Table 2. Right panel: Schematic
rapidity distributions for different particles in order to study acceptance
effects.

By construction the simulated protons, pions and clusters follow a Poisson
distribution. However the distribution of the total number of particles does
not, due to the correlation between pions and protons introduced via the
cluster decay (see Appendix F).

In experiments measurements are performed inside the finite acceptance by
imposing selection criteria in momentum space, e.g., on rapidity and/or
transverse momentum of particles. Moreover, such conditions typically lead to
different acceptances for different particle species. In order to study the
impact of the finite acceptance on the presented formalism, we introduce
arbitrary rapidity distributions for protons, pions and other particles as
illustrated in the rigth panel of Fig. 9. To this end we generate rapidity
values for protons, pions and other particles from the corresponding
distributions presented in Fig. 9.

In the left panel of Fig. 10 we present mean multiplicities of protons
produced via clusters (red circles) and independently (blue circles), while
the black circles correspond to the total number of protons. In addition, we
produce pions from clusters, and, by construction, their mean values are equal
to those of protons from clusters. The right panel of Fig. 10 shows the
second-order cumulants of protons divided by the expected value
$\bar{\kappa}_{2}[p]_{sim}$. The black stars represent those including
participant (volume) fluctuations, $\kappa_{2}[p]/\bar{\kappa}_{2}[p]_{sim}$.
The corrected cumulants $\kappa_{2}^{corr}[p]/\bar{\kappa}_{2}[p]_{sim}$ (see
Eq. 24) are shown with blue symbols, while the open red stars represent the
true reconstructed values of fluctuations of protons
$\bar{\kappa}_{2}[p]/\bar{\kappa}_{2}[p]_{sim}$ as calculated using Eq. 23.
Similar results for the third and fourth order cumulants of protons are
presented in Fig. 11 (see Eqs. 25, 26, 23). In Fig. 12 the normalized
cumulants as a function of cumulant order are presented for the 10$\%$ most
central collisions. The right panel of Fig. 12 corresponds to the full
acceptance, while in the right panel the results in the finite acceptance,
delimited as $|y|<1$, are presented. One clearly observes that in the finite
acceptance the precision of the method is significantly better. In Fig. 13 we
show the magnitude of the corresponding normalized biases,
$\Delta_{n}/\bar{\kappa}_{n}$, for the full acceptance (left panel) and for
$|y|<1$ (right panel). While the biases for the full acceptance may at first
sight appear rather large ($\sim 40\%$ ) one should realize that for the most
central events the uncorrected fourth order cumulant is more than a factor of
50 larger in magnitude than the true cumulants. In other words the proposed
corrections, while not perfect are a huge improvement of the measurement. The
situation gets better for the limited acceptance.

Figure 10: Left panel: Mean number of protons produced via clusters and
independently are presented with the red and blue circles respectively. The
black circles represent total multiplicity of protons. Pions are produced only
via clusters. Right panel: Reconstructed second-order cumulants of protons
including participant fluctuations (black stars). Corrected values for
cumulants $\kappa_{2}^{corr}[p]$, i.e., without the bias term $\Delta_{2}[p]$
are presented with blue stars, while red stars represent fully corrected,
against volume fluctuation. The results are normalized to the true second
order cumulant, $\bar{\kappa}_{2}[p]_{sim}$.

Figure 11: Left panel: Reconstructed third-order cumulants of protons
including participant fluctuations (black stars) for model B. Corrected values
for cumulants $\kappa_{3}^{corr}[p]$, i.e. without the bias term
$\Delta_{3}[p]$ are presented with blue stars, while red stars represent fully
corrected, against volume fluctuations, values
$\bar{\kappa}_{3}$[p]=$\kappa_{3}^{corr}[p]$ \+ $\Delta_{3}[p]$. Right panel:
Similar to the left panel but for the fourth-order cumulants. Note that the
values for the fourth-order cumulants with volume fluctuations (black stars)
need to be multiplied by 50 for the first three centrality classes. The
results are normalized to the true third or fourth order cumulant,
$\bar{\kappa}_{3}[p]_{sim}$ ,$\bar{\kappa}_{4}[p]_{sim}$, respectively

Figure 12: Left panel: Cumulants of protons in the full acceptance, presented
for the 10$\%$ most central collisions (cf. Figs. 10 and 11). Right panel:
Similar plot for the 10$\%$ most central collisions, but inside the finite
acceptance delimited with the $|y|<1$ criterion. The results are normalized to
the true cumulants, $\bar{\kappa}_{n}[p]_{sim}$.

Figure 13: The bias terms for model B, in the full (left panel) and finite
(right panel) acceptances. The acceptance, $|y|<1$ is introduced using
rapidity distributions of pions and protons as shown in Fig. 9.

##  VIII Software package

A Python package is provided to derive analytical formulas for both mixed and
pure cumulants of multiplicity distributions, including participant/volume
fluctuations. The correction formulas and their bias terms can be derived as
well. The dedicated graphical user interface is presented in Fig. 14 and can
be downloaded via Ref. tes .

Figure 14: The GUI for deriving analytic formulas.

##  IX Discussion and Summary

  * •

We have shown that using mixed events to determine the contributions of
wounded nucleon or volume fluctuations is equivalent to extracting the latter
from the track multiplicity distribution. However event mixing may offer an
advantage since it allows to generate an almost arbitrarily large event
ensemble with the same multiplicity distribution, and thus eliminate possible
constraints due to limited event statistics.
In either case, not all contributions can be accessed by a direct measurement.
The remaining terms lead to biases, $\Delta_{k}$, which depend on the
multiplicity distribution per wounded nucleon. These biases are, however,
parametrically suppressed by powers of $\left\langle
N\right\rangle/\left\langle M\right\rangle$ depending on the order of the
cumulants. The biases are also small if the multiplicity distribution per
wounded nucleon is close to Poisson. In addition, we suggest to constrain
these biases in experiment with fits to the observed multiplicity distribution
within the wounded-nucleon model.

  * •

We have worked here within the wounded nucleon model to formulate volume
fluctuations. Alternatively, one may introduce generic volume fluctuations as
done e.g. in Skokov  _et al._ (2013). It is easy to show (see Appendix D) that
the resulting expressions for the corrected cumulants, $\kappa_{j}^{corr}$,
and the biases, $\Delta_{j}$, are identical to those derived here, i.e. Eqs.
(24-26).

  * •

We note that one gets similar expressions for the fluctuations from the
wounded nucleons, Eqs. (20-22), in terms of the factorial cumulants of, for
example, pions instead of the total track multiplicity. This has the advantage
that the corrections do not involve the particles of interest, protons, in our
case. Of course, this approach requires that sufficiently many pions are
produced and thus may be limited to collisions at higher energies.

  * •

We have checked that the proposed method also works if the multiplicity
distribution is determined for a different acceptance than the particle
distribution of interest. In this case all quantities in the expression for
the corrected cumulants, Eqs.(24-26) involving the multiplicity should be
evaluated in the multiplicity acceptance while all quantities involving the
particles of interest, such as $\left\langle N\right\rangle$ or
$\kappa_{j}[N]$ should be determined in the particle acceptance.

  * •

We have verified that the proposed method is not affected by potentially
different rapidity distributions for different particle types.

  * •

We note that the correction term and bias for the second-order cumulant
depends on the properties of the multiplicity distribution only while those
for higher-order cumulants also involve the (uncorrected) cumulants of
interest (at a lower order), $\kappa_{n}[N]$ (see Eq. (25,26).

  * •

The proposed method is also applicable for mixed cumulants, such as the
covariance between two particle species. The relevant formulas for mixed
cumulants between species $A$ and $B$ up to $\kappa_{2,2}[A,B]$ are given in
Appendix G.

The corresponding relations for correction and bias terms for mixed cumulants
of any order can be obtained with the provided software package tes .

  * •

Here we have not explicitly discussed corrections for net proton cumulants.
However, those can be easily obtained using the corrections to the mixed
cumulants provided in Appendix G. For example, the second-order cumulant of
the net proton distribution is given by

$\displaystyle\kappa_{2}[N-\bar{N}]=\kappa_{2}[N]+\kappa_{2}[\bar{N}]-2\kappa_{1,1}[N,\bar{N}]$
(39)

The corrected cumulant is the (using Eqs. (24, 101)

$\displaystyle\kappa_{2}^{corr}[N-\bar{N}]=\kappa_{2}[N-\bar{N]}-\frac{\left(\left\langle
N\right\rangle-\left\langle\bar{N}\right\rangle\right)^{2}}{\left\langle
M\right\rangle^{2}}C_{2}[M]$ (40)

with the bias

$\displaystyle\Delta=\frac{\left(\left\langle
N\right\rangle-\left\langle\bar{N}\right\rangle\right)^{2}}{\left\langle
M\right\rangle^{2}}\bar{C}_{2}[M]$ (41)

For systems at vanishing baryon number chemical potential, such as those
created at very high collision energies, $\left\langle
N\right\rangle=\left\langle\bar{N}\right\rangle$ so that the corrected
cumulant is identical to the measured one (as discussed in Skokov  _et al._
(2013) ) and the bias vanishes.

In summary, we have presented a method to correct experimentally measured
particle number cumulants for the effect of participant or volume
fluctuations. The essential idea is to extract the contribution from the
volume fluctuations from the distribution of charged particles which, after
appropriate re-scaling, may be subtracted from the measured cumulants of
interest. Our proposed method is not exact as there remains a bias or remnant
which can not be accessed directly from experiment. However, we have shown by
model calculations that this bias is very small compared to the contribution
from participant fluctuations and we hence consider our method an important
step towards measuring the true dynamical fluctuations of the system.

###### Acknowledgements.

V.K would like to thank GSI and the Institute for Nuclear Theory at the
University of Washington for their kind hospitality and stimulating research
environment. V.K. has been supported by the U.S. Department of Energy, Office
of Science, Office of Nuclear Physics, under contract number DE-
AC02-05CH11231, by the INT’s U.S. Department of Energy grant No. DE-
FG02-00ER41132, and by the ExtreMe Matter Institute EMMI at the GSI
Helmholtzzentrum für Schwerionenforschung, Darmstadt, Germany.

##  Appendix A Wounded-nucleon model

Here we briefly discuss the wounded-nucleon model following the Appendix of
Ref. (Bzdak  _et al._ , 2017). The wounded-nucleon model assumes that
particles are produced by independent sources, called wounded nucleons or
participants. Therefore, the probability to find $A$ particles of type $A$ and
$B$ particles of type $B$ can be written as

$P\left(A,B\right)=\sum_{w}W\left(w\right)\sum_{a_{1},\cdots
a_{w}}\sum_{b_{1},\cdots b_{w}}p\left(a_{1},b_{1}\right)\cdots
p\left(a_{w},b_{w}\right)\,\delta_{A,\sum_{k=1}^{w}a_{k}}\delta_{B,\sum_{k=1}^{w}b_{k}}.$
(42)

Here $W(w)$ denotes the probability to have $w$ wounded nucleons, and $p(a,b)$
is the probability to have $a$ particles of type $A$ and $b$ particles of type
$B$ from one wounded nucleon. The moment-generating function,
$h\left(t_{A},t_{B}\right)$ is then

$\displaystyle H\left(t_{A},t_{B}\right)$
$\displaystyle=\sum_{A,B}e^{t_{A}A}e^{t_{b}B}P\left(A,B\right)$
$\displaystyle=\sum_{w}W\left(w\right)\sum_{a_{1},\cdots
a_{w}}\sum_{b_{1},\cdots b_{w}}p\left(a_{1},b_{1}\right)\cdots
p\left(a_{w},b_{w}\right)e^{t_{A}\sum_{k=1}^{w}a_{k}}e^{t_{B}\sum_{k=1}^{w}b_{k}}$
$\displaystyle=\sum_{w}W\left(w\right)\sum_{a_{1},b_{1}}p\left(a_{1},b_{1}\right)e^{t_{A}a_{1}+t_{b}b_{1}}\cdots\sum_{a_{w},b_{w}}p\left(a_{w},b_{w}\right)e^{t_{A}a_{w}+t_{B}b_{w}}$
$\displaystyle=\sum_{w}W\left(w\right)\left[\sum_{a,b}p\left(a,b\right)e^{t_{A}a+t_{B}b}\right]^{w}$
$\displaystyle=\sum_{w}W\left(w\right)\left[h_{w}\left(t_{A},t_{B}\right)\right]^{w}$
$\displaystyle=\sum_{w}W\left(w\right)e^{w\,g_{w}(t_{A},t_{B})}$ (43)

where
$h_{w}\left(t_{A},t_{B}\right)=\sum_{a,b}p\left(a,b\right)e^{t_{A}a+t_{B}b}$
is the moment generating function and
$g_{w}(t_{A},t_{B})=\ln\left[h_{w}\left(t_{A},t_{B}\right)\right]$ the
cumulant-generating function for one wounded nucleon, respectively. The
cumulant-generating function,
$G\left(t_{A},t_{B}\right)=\ln\left(h\left(t_{A},t_{B}\right)\right)$, is then
given by

$G\left(t_{A},t_{B}\right)=\ln\left[H\left(t_{A},t_{B}\right)\right]=\ln\left[\sum_{w}W\left(w\right)e^{w\,g_{w}(t_{1},t_{2})}\right]=G_{W}\left(g_{w}(t_{A},t_{B})\right)$
(44)

where $G_{W}(t)=\ln\left[\sum_{w}W\left(w\right)e^{w\,t}\right]$ in the
cumulant-generating function for the wounded-nucleon distribution, $W(w)$. We
note, that $g_{w}(0,0)=G_{W}(0)=0$ by construction. The cumulants are then
obtained as

$\kappa_{j,k}[A,B]=\left.\frac{\partial^{j}\partial^{k}}{\partial
t_{A}\partial t_{B}}g\left(t_{A},t_{B}\right)\right|_{t_{A}=t_{B}=0}.$ (45)

For example:

$\displaystyle\kappa_{1}[A]$ $\displaystyle=\frac{\partial}{\partial
t_{A}}\left.G(t_{A},0)\right|_{t_{A}=0}=\frac{dG_{w}}{dg_{w}}\left.\frac{dg_{w}(t_{A},0)}{dt_{A}}\right|_{t_{A}=0}=\left.\frac{dG_{w}}{dg_{w}}\right|_{g_{w}=0}\left.\frac{dg_{w}(t_{A},0)}{dt_{A}}\right|_{t_{A}=0}$
$\displaystyle=\kappa_{1}[N_{w}]\kappa_{1}[a]=\left\langle
N_{w}\right\rangle\left\langle a\right\rangle$ (46)

where $\kappa_{1}[n]=\left\langle a\right\rangle$ denotes the mean number of
particles of type $A$ per wounded nucleon and $\kappa_{1}[w]=\left\langle
N_{w}\right\rangle$ the mean number of wounded nucleons. The second-order
cumulant is

$\displaystyle\kappa_{2}[A]$ $\displaystyle=\frac{\partial^{2}}{\partial
t_{A}^{2}}\left.G(t_{A},0)\right|_{t_{A}=0}$
$\displaystyle=\frac{d^{2}G_{w}}{dg_{w}^{2}}\left.\left(\frac{dg_{w}(t_{A},0)}{dt_{A}}\right)^{2}\right|_{t_{A}=0}+\frac{dG_{w}}{dg_{w}}\left.\frac{d^{2}g_{w}(t_{A},0)}{dt_{A}^{2}}\right|_{t_{1}=0}$
$\displaystyle=\kappa_{2}[N_{w}]\kappa_{1}[a]^{2}+\kappa_{1}[N_{w}]\kappa_{2}[a]=\kappa_{2}[N_{w}]\left\langle
a\right\rangle^{2}+\left\langle N_{w}\right\rangle\kappa_{2}[a]$ (47)

The covariance is

$\displaystyle cov[A,B]$ $\displaystyle=\frac{\partial}{\partial t_{A}\partial
t_{B}}\left.G(t_{A},t_{B})\right|_{t_{A},t_{B}=0}=\frac{\partial}{\partial
t_{B}}\left.\left(\frac{dG_{w}}{dg_{w}}\frac{\partial
g_{w}(t_{A},t_{B})}{\partial t_{A}}\right)\right|_{t_{A},t_{B}=0}$
$\displaystyle=\left.\frac{d^{2}G_{w}}{dg_{w}^{2}}\frac{\partial
g_{w}(t_{A},t_{B})}{\partial t_{A}}\frac{\partial g_{w}(t_{A},t_{B})}{\partial
t_{B}}\right|_{t_{A},t_{B}=0}+\left.\left(\frac{dG_{w}}{dg_{w}}\frac{\partial^{2}g_{w}(t_{A},t_{B})}{\partial
t_{A}\partial t_{B}}\right)\right|_{t_{A},t_{B}=0}$
$\displaystyle=\kappa_{2}[N_{W}]\kappa_{1}[a]\kappa_{1}[a]+\kappa_{1}[N_{w}]cov[a,b]$
$\displaystyle=\kappa_{2}[N_{w}]\left\langle a\right\rangle\left\langle
b\right\rangle+\left\langle N_{w}\right\rangle cov[a,b]$ (48)

Using the relation between the cumulant and factorial cumulant generating
function, Eq. 53 , the factorial cumulant generating function is given by

$G_{F}\left(z_{\rm A},z_{\rm B}\right)=G\left(\ln(z_{\rm A}),\ln(z_{\rm
B})\right)=G_{W}\left(g_{w}\left(\ln(z_{\rm A}),\ln(z_{\rm
B})\right)\right)=G_{W}\left(g_{F,w}\left(z_{\rm A},z_{\rm B}\right)\right),$
(49)

with

$g_{F,w}\left(z_{\rm A},z_{\rm B}\right)=\sum_{a,b}p\left(a,b\right)z_{\rm
A}^{a}z_{\rm B}^{b}$

the factorial cumulant generating function for the distribution of one wounded
nucleon, $p\left(a,b\right)$. The structure is the same as for the cumulant
generating function, except that now the argument of the wounded-nucleon
cumulant-generating function is the factorial cumulant generating function,
$g_{F,w}$. Thus the factorial cumulants are easily obtained by simply
replacing all the cumulants of the particle distribution for a given wounded
nucleon, $\kappa_{i,j}[a,b]$ with the corresponding factorial cumulants, with
$C_{1}[X]=\kappa_{1}[X]=\left\langle X\right\rangle$

$\displaystyle C_{1}[A]$ $\displaystyle=\frac{\partial}{\partial z_{\rm
A}}\left.g_{F}(z_{\rm A},1)\right|_{z_{\rm
A}=1}=\kappa_{1}[N_{w}]C_{1}[a]=\left\langle N_{w}\right\rangle\left\langle
a\right\rangle$ $\displaystyle C_{2}[A]$
$\displaystyle=\kappa_{2}[N_{w}]\left\langle a\right\rangle^{2}+\left\langle
N_{w}\right\rangle C_{2}[a]$ $\displaystyle C_{1,1}[A,B]$
$\displaystyle=cov[A,B]$ (50)

##  Appendix B Cumulant and factorial cumulant generating functions

Given a multiplicity distribution for particles of type $A$ and $B$,
$P\left(A,B\right)$ the generating functions for cumulants, $g(t_{A},t_{B})$,
and factorial cumulants. $g_{F}(t_{A},t_{B})$ are given by

$\displaystyle g\left(t_{A},t_{B}\right)$
$\displaystyle=\ln[\sum_{A,B}P\left(A,B\right)e^{t_{A}A}e^{t_{B}B}$ (51)
$\displaystyle g_{F}\left(z_{\rm A},z_{\rm B}\right)$
$\displaystyle=\ln[\sum_{A,B}P\left(A,B\right)(z_{\rm A})^{A}(z_{\rm B})^{B}.$
(52)

By construction, $g\left(t_{A}=0,t_{B}=0\right)=0$ and $g_{F}\left(z_{\rm
A}=1,z_{\rm B}=1\right)=0$. Cumulants of order $(i,j)$, $\kappa_{j,k}[A,B],$
are then obtained through

$\kappa_{j,k}[A,B]=\left.\frac{\partial^{j}\partial^{k}}{\partial
t_{A}\partial t_{B}}g\left(t_{A},t_{B}\right)\right|_{t_{A}=t_{B}=0},$

while the factorial cumulants, $C_{j,k}[A,B]$, are given by

$C_{j,k}[A,B]=\left.\frac{\partial^{j}\partial^{k}}{\partial z_{\rm A}\partial
z_{\rm B}}g\left(z_{\rm A},zt\right)\right|_{z_{\rm A}=z_{\rm B}=1}.$

The generating functions are related through

$g_{F}\left(z_{\rm A},z_{\rm B}\right)=g\left[\ln(z_{\rm A}),\ln(z_{\rm
B})\right]$ (53)

or vice versa

$g\left(t_{A},t_{B}\right)=g_{F}\left(e^{t_{B}},e^{t_{B}}\right)$ (54)

These relations may also be used to convert cumulants into factorial cumulants
and vice versa. For example, for the diagonal cumulants, $\kappa_{n}[A]$ we
have

$\kappa_{n}[A]=\sum_{j=1}^{n}S\left(n,j\right)C_{j}[A]$ (55)

where $S(n,j)$ denotes the Sterling numbers of the second kind. The inverse
relation is

$C_{n}[A]=\sum_{j=1}^{n}B_{n,j}\left(1,-1,2,\ldots,(-1)^{j-1}(n-j+1)!\right)\,\kappa_{j}[A]$
(56)

with $B_{n,j}$ being Bell polynomials. For the first four orders this
evaluates to

$\displaystyle\kappa_{1}$ $\displaystyle=C_{1}$ $\displaystyle\kappa_{2}$
$\displaystyle=C_{1}+C_{2}$ $\displaystyle\kappa_{3}$
$\displaystyle=C_{1}+3C_{2}+C_{3}$ $\displaystyle\kappa_{4}$
$\displaystyle=C_{1}+7C_{2}+6C_{3}+C_{4}$ (57)

and

$\displaystyle C_{2}$ $\displaystyle=\kappa_{2}-\kappa_{1}$ $\displaystyle
C_{3}$ $\displaystyle=2\kappa_{1}-3\kappa_{2}+\kappa_{3}$ $\displaystyle
C_{4}$ $\displaystyle=-6\kappa_{1}+11\kappa_{2}-6\kappa_{3}+\kappa_{4}$ (58)

##  Appendix C Multiplicity Distribution

The multiplicity distribution, $P(M)$, is given by summing over all (charged)
particles,

$\displaystyle P_{M}\left(M\right)$
$\displaystyle=\sum_{A,B,X}P\left(A,B,X\right)\delta_{M,A+B+X},$ (59)

where

$P\left(A,B,X\right)=\sum_{w}W\left(w\right)\sum_{a_{1},\cdots
a_{w}}\sum_{b_{1},\cdots b_{w}}\sum_{x_{1}.\cdots
x_{w}}p\left(a_{1},b_{1},x_{1}\right)\cdots
p\left(a_{w},b_{w},x_{w}\right)\,\delta_{A,\sum_{k=1}^{w}a_{k}}\delta_{B,\sum_{k=1}^{w}a_{k}}\delta_{X,\sum_{k=1}^{w}x_{k}}$

is the distribution of particles of type $A$, $B$ and all others, denoted by
$X$. The distribution for particles $A$ and $B$ defined in Appendix A are then
given by $P\left(A,B\right)=\sum_{X=0}^{\infty}P\left(A,B,X\right)$ while that
for the particles per wounded nucleons are given by
$p\left(a,b\right)=\sum_{x=0}^{\infty}p\left(a,b,x\right)$

The moment generating function is then given by (proceeding analogously to Eq.
43:

$\displaystyle H_{M}(t)$
$\displaystyle=\sum_{M}P_{M}(M)\,e^{tM}=\sum_{M,A,B,X}P(A,B,X)\delta_{M,A+B+X}\,e^{tM}=\sum_{A,B,X}P(A,B,X)\,e^{t(A+B+X)}$
$\displaystyle=\sum_{w}W\left(w\right)\sum_{a_{1},\cdots
a_{w}}\sum_{b_{1},\cdots b_{w}}\sum_{x_{1}.\cdots
x_{w}}p\left(a_{1},b_{1},x_{1}\right)\cdots
p\left(a_{w},b_{w},x_{w}\right)e^{t\sum_{k=1}^{w}a_{k}}e^{t\sum_{k=1}^{w}b_{k}}e^{t\sum_{k=1}^{w}x_{k}}$
$\displaystyle=\sum_{w}W\left(w\right)\sum_{a_{1}}\sum_{b_{1}}\sum_{x_{1}}p\left(a_{1},b_{1},x_{1}\right)e^{t(a_{1}+b_{1}+c_{1})}\cdots\sum_{a_{w}}\sum_{b_{w}}\sum_{x_{w}}p\left(a_{w},b_{w},x_{w}\right)e^{t(a_{w}+b_{w}+c_{w})}$
$\displaystyle=\sum_{w}W\left(w\right)\left[\sum_{a}\sum_{b}\sum_{x}p\left(a,b,x\right)e^{t(a+b+c)}\right]^{w}$
$\displaystyle=\sum_{w}W\left(w\right)\left[h_{m,w}(t)\right]^{w}$
$\displaystyle=\sum_{w}W\left(w\right)e^{wg_{m,w}(t)}$

where, $h_{m,w}(t)$ the moment generating function and
$g_{m,w}(t)=\ln\left[h_{m,w}(t)\right]$ the cumulant generating function of
the multiplicity distribution for _one_ wounded nucleon,
$p(m)=\sum_{a,b,x}p(a,b,x)\delta_{m,a+b+c}$. The cumulant generating function,
$G_{M}(t)$, for the multiplicity distribution, $P\left(M\right)$, is then
given by

$G_{M}(t)=\ln\left[H_{M}(t)\right]=\ln\left[\sum_{w}W\left(w\right)e^{wg_{m,w}(t)}\right]=G_{W}\left(g_{m,w}(t)\right)$
(60)

The cumulants of the multiplicity distribution are given by (following the
analogous Eqs. 46 and 47)

$\displaystyle\kappa_{1}[M]$ $\displaystyle=\left\langle
N_{w}\right\rangle\left\langle m\right\rangle$ $\displaystyle\kappa_{2}[M]$
$\displaystyle=\kappa_{2}[N_{w}]\left\langle m\right\rangle^{2}+\left\langle
N_{w}\right\rangle\kappa_{2}[m]$ (61)

where $\kappa_{i}[m]$ denote the cumulants of the multiplicity distribution of
_one_ wounded nucleon and $\kappa_{i}\left[N_{W}\right]$ those of the wounded
nucleon distribution. Analogous to Eq. 49 the factorial cumulant generating
function is readily obtained

$G_{F,M}\left(z,\right)=G_{M}\left(\ln(z)\right)=G_{W}\left(g_{m,w}\left(\ln(z)\right)\right)=G_{W}\left(g_{F,m,w}\left(z\right)\right),$
(62)

with

$g_{F,m,w}(z)=\sum_{m}p(m)z^{m}$

the factorial cumulant generating function for the multiplicity distribution
of one nucleon, $p(m)$. Again, the factorial cumulants are obtained by
replacing the cumulants of the distribution $p(m)$, $\kappa_{i}[m]$ with the
corresponding factorial cumulants, $C_{i}[m]$, in Eq. 61 by the factorial
cumulants

$\displaystyle C_{1}[M]$ $\displaystyle=\left\langle
N_{w}\right\rangle\left\langle m\right\rangle$ $\displaystyle C_{2}[M]$
$\displaystyle=\kappa_{2}[N_{w}]\left\langle m\right\rangle^{2}+\left\langle
N_{w}\right\rangle C_{2}[m]$ (63)

##  Appendix D Wounded Nucleon vs Volume Fluctuations

Here we will discuss the relation between wounded nucleon fluctuations and so-
called volume fluctuations as they are discussed e.g. in Skokov  _et al._
(2013). Following Ref. Skokov  _et al._ (2013) the cumulant generating
function is given by

$\Phi(t)=\ln\text{$\left[\int dVP(V)e^{V\xi(t)}\right]$=$\chi^{V}(\xi(t))$}$
(64)

with $\chi^{V}(t)$ the cumulant generating function for the distribution of
volumes, $P(V)$, and

$\xi(t)=\frac{1}{V}\ln\left[\sum_{N}p(N;V)e^{Nt}\right]$ (65)

the generating function for scaled cumulants, $\kappa/V$, given for the
distribution of particles at fixed volume, $p(N;V)$. Then, for a fixed volume
$V$, the scaled cumulants are given by

$\frac{\kappa_{j}}{V}=\frac{\partial^{j}}{\partial\,t^{j}}\left.\xi(t)\right|_{t=0}.$
(66)

For the wounded nucleon model we have (see Appendix A)

$G(t)=\ln\left[\sum_{w}W\left(w\right)e^{w\,g_{w}(t)}\right]=G_{W}(g_{w}(t))$
(67)

with $G_{W}(t)$ the cumulant generating function for the wounded nucleon
distribution, $W(w),$ and

$g_{w}(t)=\ln\left[\sum_{n}p(n)e^{nt}\right]$ (68)

the generating function for the distribution of particles for one wounded
nucleon. The cumulants per wounded nucleons for a fixed number of wounded
nucleons, $N_{w},$ are given by

$\frac{\kappa_{j}[N]}{N_{w}}=\kappa_{j}[n]=\frac{\partial^{j}}{\partial\,t^{j}}\left.g_{w}(t)\right|_{t=0}.$
(69)

Comparing the above expressions, one finds that the cumulants for volume
fluctuations can be obtained from those for the wounded-nucleon number by the
following replacements

$\displaystyle\kappa_{j}[N_{w}]$ $\displaystyle\rightarrow\kappa_{j}[V]$
$\displaystyle\bar{\kappa}_{j}[N]=\left\langle
N_{w}\right\rangle\kappa_{j}[N]$ $\displaystyle\rightarrow\left\langle
V\right\rangle\frac{\kappa_{j}}{V}$

Indeed comparing the second-order cumulants for both scenarios we have

$\displaystyle\kappa_{2}[N]$ $\displaystyle=\left\langle
N_{w}\right\rangle\kappa_{2}[n]+\left\langle
n\right\rangle^{2}\kappa_{2}[N_{w}]=\bar{\kappa}_{2}[N]+\left\langle
N\right\rangle^{2}\frac{\kappa_{2}[N_{w}]}{\left\langle
N_{w}\right\rangle^{2}}$ $\displaystyle\kappa_{2}[N]$
$\displaystyle=\left\langle
V\right\rangle\frac{\kappa_{2}}{V}+\left(\frac{\kappa_{1}}{V}\right)^{2}\kappa_{2}[V]=\bar{\kappa}_{2}[N]+\left\langle
N\right\rangle^{2}\frac{\kappa_{2}[V]}{\left\langle V\right\rangle^{2}}$

where in the second line we used $\kappa_{1}=\left\langle N\right\rangle$ and
$\bar{\kappa}_{2}=\left\langle V\right\rangle\frac{\kappa_{2}}{V}$. Obviously,
analogous replacements also hold for the factorial cumulants

$\bar{C}_{j}[N]=\left\langle N_{w}\right\rangle
C_{j}[n]\rightarrow\left\langle V\right\rangle\frac{C_{j}}{V}$

with $C_{j}/V$ the volume scaled factorial cumulants.

##  Appendix E Results for factorial cumulants

Here we provide the formulas for the corrected factorial cumulants,
$C_{k}^{corr}$ and $c_{k}^{corr}$, and the the associated biases,
$\Delta_{k,F}$ and $\delta_{k,F}$. Both the factorial cumulants and the biases
are related to the corresponding cumulants via the linear relation Eqs. 57 and
58. The corrected factorial cumulants and the associated biases are

$\displaystyle C_{2}^{corr}$ $\displaystyle=C_{2}[N]-\frac{\langle
N\rangle^{2}}{\langle M\rangle^{2}}C_{2}[M]$ (70) $\displaystyle C_{3}^{corr}$
$\displaystyle=C_{3}[N]-\frac{3C_{2}[M]C_{2}[N]\langle N\rangle}{\langle
M\rangle^{2}}+\frac{3C_{2}[M]{}^{2}\langle N\rangle^{3}}{\langle
M\rangle^{4}}-\frac{C_{3}[M]\langle N\rangle^{3}}{\langle M\rangle^{3}}$ (71)
$\displaystyle C_{4}^{corr}$
$\displaystyle=C_{4}[N]-\left(\frac{6C_{2}[N]\langle
N\rangle^{2}\left(C_{3}[M]\langle M\rangle-3C_{2}[M]{}^{2}\right)}{\langle
M\rangle^{4}}+\frac{4C_{2}[M]C_{3}[N]\langle N\rangle}{\langle
M\rangle^{2}}\right.$ (72)
$\displaystyle\left.+\frac{3C_{2}[M]C_{2}[N]{}^{2}}{\langle
M\rangle^{2}}+\frac{\langle N\rangle^{4}\left(-10C_{3}[M]C_{2}[M]\langle
M\rangle+C_{4}[M]\langle M\rangle^{2}+15C_{2}[M]{}^{3}\right)}{\langle
M\rangle^{6}}\right)$ (73)

$\displaystyle\Delta_{2,F}$ $\displaystyle=\frac{\langle N\rangle^{2}}{\langle
M\rangle^{2}}\bar{C}_{2}[M]$ (74) $\displaystyle\Delta_{3,F}$
$\displaystyle=\bar{C}_{2}[M]\left(\frac{3C_{2}[N]\langle N\rangle}{\langle
M\rangle^{2}}-\frac{3C_{2}[M]\langle N\rangle^{3}}{\langle
M\rangle^{4}}\right)+\frac{\bar{C}_{3}[M]\langle N\rangle^{3}}{\langle
M\rangle^{3}}$ (75) $\displaystyle\Delta_{4,F}$
$\displaystyle=\bar{C}_{2}[M]\left(\frac{4C_{3}[N]\langle N\rangle}{\langle
M\rangle^{2}}+\frac{3C_{2}[N]{}^{2}}{\langle
M\rangle^{2}}-\frac{18C_{2}[M]C_{2}[N]\langle N\rangle^{2}}{\langle
M\rangle^{4}}-\frac{4C_{3}[M]\langle N\rangle^{4}}{\langle
M\rangle^{5}}+\frac{15C_{2}[M]{}^{2}\langle N\rangle^{4}}{\langle
M\rangle^{6}}\right)$ (76)
$\displaystyle+\bar{C}_{3}[M]\left(\frac{6C_{2}[N]\langle
N\rangle^{2}}{\langle M\rangle^{3}}-\frac{6C_{2}[M]\langle
N\rangle^{4}}{\langle M\rangle^{5}}\right)+\frac{\bar{C}_{4}[M]\langle
N\rangle^{4}}{\langle M\rangle^{4}}$ (77)

For the scaled factorial cumulants we have

$\displaystyle c_{2}^{corr}$ $\displaystyle=c_{2}[N]-\frac{\left\langle
N\right\rangle}{\left\langle M\right\rangle}c_{2}[M]$ (78) $\displaystyle
c_{3}^{corr}$ $\displaystyle=c_{3}[N]-3\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}c_{2}[M]c_{2}[N]+\left(\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}\right)^{2}\left(3c_{2}[M]^{2}-c_{3}[M]\right)$ (79)
$\displaystyle c_{4}^{corr}$ $\displaystyle=c_{4}[N]-\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}\left(3c_{2}[M]c_{2}[N]{}^{2}+4c_{2}[M]c_{3}[N]\right)$ (80)
$\displaystyle+\left(\frac{\left\langle N\right\rangle}{\left\langle
M\right\rangle}\right)^{2}\left(18c_{2}[M]^{2}c_{2}[N]-6c_{3}[M]c_{2}[N]\right)$
(81) $\displaystyle+\left(\frac{\left\langle N\right\rangle}{\left\langle
M\right\rangle}\right)^{3}\left(-15c_{2}[M]{}^{3}+10c_{2}[M]c_{3}[M]-c_{4}[M]\right)$
(82)

$\displaystyle\delta_{2,F}$ $\displaystyle=\frac{\left\langle
N\right\rangle}{\left\langle M\right\rangle}\bar{c}_{2}[M]$ (83)
$\displaystyle\delta_{3,F}$ $\displaystyle=3\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}c_{2}[N]\bar{c}_{2}[M]+\left(\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}\right)^{2}\left(\bar{c}_{3}[M]-3c_{2}[M]\bar{c}_{2}[M]\right)$
(84) $\displaystyle\delta_{4,F}$ $\displaystyle=\frac{\left\langle
N\right\rangle}{\left\langle
M\right\rangle}\left(3c_{2}[N]{}^{2}\bar{c}_{2}[M]+4c_{3}[N]\bar{c}_{2}[M]\right)$
(85) $\displaystyle+\left(\frac{\left\langle N\right\rangle}{\left\langle
M\right\rangle}\right)^{2}\left(6c_{2}[N]\bar{c}_{3}[M]-18c_{2}[M]c_{2}[N]\bar{c}_{2}[M]\right)$
(86) $\displaystyle+\left(\frac{\left\langle N\right\rangle}{\left\langle
M\right\rangle}\right)^{3}\left(15c_{2}[M]{}^{2}\bar{c}_{2}[M]-6c_{2}[M]\bar{c}_{3}[M]-4c_{3}[M]\bar{c}_{2}[M]+\bar{c}_{4}[M]\right)$
(87)

##  Appendix F Particle production through clusters

Let us assume that particles are produced via clusters and that each cluster
further decays into two particles. Moreover, clusters are generated from a
Poisson distribution. As each cluster decays into two particles the
probability of measuring $k$ particles is equivalent to measuring $k/2$
clusters and can be presented as:

$p\left(k;\langle N_{cl}\rangle\right)=e^{-\langle N_{cl}\rangle}\frac{\langle
N_{cl}\rangle^{k/2}}{(k/2)!}$ (88)

The corresponding moment generating function reads:

$M(t)=\sum_{k/2=0}^{\infty}e^{tk}e^{-\langle N_{cl}\rangle}\frac{\langle
N_{cl}\rangle^{k/2}}{(k/2)!}=e^{\langle N_{cl}\rangle\left(e^{2t}-1\right)},$
(89)

where $\langle N_{cl}\rangle$ denotes mean number of clusters produced.

The cumulants of total particle number $k$ can be computed as:

$\kappa_{n}[k]=\left.\frac{\partial ln(M(t))}{dt}\right\rvert_{t=0}$ (90)

For the first two cumulants one gets:

$\kappa_{1}[k]=2\langle N_{cl}\rangle$ (91)

$\kappa_{2}[k]=4\langle N_{cl}\rangle$ (92)

One clearly sees from Eqs. 91 and 92 that $\kappa_{1}(k)\neq\kappa_{2}(k)$,
i.e the total number of particles does not follow a Poisson distribution,
although the clusters do. Moreover, one observes that particle production
through clusters enhances fluctuations. In general, for clusters following a
Poisson distribution and decaying into $m$ particles, the cumulants of total
particle number can be written as:

$\kappa_{n}[k]=m^{n}\langle N_{cl}\rangle$ (93)

##  Appendix G Mixed Cumulants

Here we provide the relevant formulas for mixed cumulants. Given the
generating function, Eq. 44, the mixed cumulants for particles of type $A$ and
$B$ are given by (see Eq. 45)

$\kappa_{j,k}[A,B]=\left.\frac{\partial^{j}\partial^{k}}{\partial
t_{A}\partial t_{B}}g\left(t_{A},t_{B}\right)\right|_{t_{A}=t_{B}=0}$ (94)

The explicit formulas for the four lowest-order mixed cumulants are:

$\displaystyle\kappa_{1,1}[A,B]$
$\displaystyle=\bar{\kappa}_{1,1}[A,B]+\langle A\rangle\langle
B\rangle\frac{\kappa_{2}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{2}}$ (95) $\displaystyle\kappa_{2,1}[A,B]$
$\displaystyle=\bar{\kappa}_{2,1}[A,B]+\left(2\langle
A\rangle\bar{\kappa}_{1,1}[A,B]+\langle
B\rangle\bar{\kappa}_{2,0}[A,B]\right)\frac{\kappa_{2}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{2}}+\langle A\rangle^{2}\langle
B\rangle\frac{\kappa_{3}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{3}}$ (96) $\displaystyle\kappa_{1,2}[A,B]$
$\displaystyle=\bar{\kappa}_{1,2}[A,B]+\left(\langle
A\rangle\bar{\kappa}_{0,2}[A,B]+2\langle
B\rangle\bar{\kappa}_{1,1}[A,B]\right)\frac{\kappa_{2}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{2}}+\langle A\rangle\langle
B\rangle^{2}\frac{\kappa_{3}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{3}}$ (97) $\displaystyle\kappa_{2,2}[A,B]$
$\displaystyle=\bar{\kappa}_{2,2}[A,B]+\left(\langle
A\rangle^{2}\bar{\kappa}_{0,2}[A,B]+4\langle A\rangle\langle
B\rangle\bar{\kappa}_{1,1}[A,B]+\langle
B\rangle^{2}\bar{\kappa}_{2,0}[A,B]\right)\frac{\kappa_{3}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{3}}$ $\displaystyle+\left(2\langle
A\rangle\bar{\kappa}_{1,2}[A,B]+2\langle
B\rangle\bar{\kappa}_{2,1}[A,B]+2\bar{\kappa}_{1,1}[A,B]{}^{2}+\bar{\kappa}_{0,2}[A,B]\bar{\kappa}_{2,0}[A,B]\right)\frac{\kappa_{2}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{2}}$ $\displaystyle+\langle A\rangle^{2}\langle
B\rangle^{2}\frac{\kappa_{4}\left[N_{W}\right]}{\left\langle
N_{W}\right\rangle{}^{4}}$ (98)

where, analogous to the notation for the regular cumulants,
$\bar{\kappa}_{j,k}[A,B]$ denotes the mixed cumulant for constant number of
wounded nucleons $\left\langle N_{W}\right\rangle$. Note, that
$\kappa_{j,0}[A,B]=\kappa_{j}[A]$ and $\kappa_{0,j}[A,B]=\kappa_{j}[B]$
correspond to the regular cumulant for particles of type $A$ and $B$
respectively. The firs order mixed cumulant, $\kappa_{1,1}[A,B]=cov[A,B]$ is
also referred to as the co-variance between the distributions of particles $A$
and $B$. In order to obtain the corrected mixed cumulants we proceed in the
same fashion as for the regular cumulant. We express the terms involving
cumulants of the wounded nucleons, $\kappa_{i}\left[N_{W}\right]/\left\langle
N_{W}\right\rangle^{i}$ in terms of the factorial cumulants of the
multiplicity distribution (See Eqs. 20-22) and solve for the the mixed
cumulants with fixed number of wounded nucleons, $\bar{\kappa}_{j,k}[A,B]$.
Again, the results are given in the form

$\displaystyle\bar{\kappa}_{j,k}[A,B]=\kappa_{j,k}^{corr}[A,B]+\Delta_{j,k}$
(99)

where $\kappa_{j,k}^{corr}[A,B]$ are the cumulants including the measurable
corrections and $\Delta_{j,k}$ are the corresponding biases due to quantities
which are not directly measurable.

$\displaystyle\kappa_{1,1}^{corr}[A,B]$
$\displaystyle=\kappa_{1,1}[A,B]-\frac{\langle A\rangle\langle
B\rangle}{\langle M\rangle^{2}}C_{2}[M]$ (100) $\displaystyle\Delta_{1,1}$
$\displaystyle=\frac{\langle A\rangle\langle B\rangle}{\langle
M\rangle^{2}}\bar{C}_{2}[M]$ (101) $\displaystyle\kappa_{2,1}^{corr}[A,B]$
$\displaystyle=\kappa_{2,1}[A,B]-\frac{\langle B\rangle
C_{2}[M]\bar{\kappa}_{2,0}[A,B]}{\langle M\rangle^{2}}-\frac{2\langle A\rangle
C_{2}[M]\kappa_{1,1}[A,B]}{\langle M\rangle^{2}}+\frac{\langle
A\rangle^{2}\langle B\rangle\left(2C_{2}[M]{}^{2}-\langle M\rangle
C_{3}[M]\right)}{\langle M\rangle^{4}}$ (102) $\displaystyle\Delta_{2,1}$
$\displaystyle=\frac{1}{\langle M\rangle^{4}}\left[2\langle A\rangle\langle
M\rangle^{2}\bar{C}_{2}[M]\kappa_{1,1}[A,B]+\langle B\rangle\langle
M\rangle^{2}\bar{C}_{2}[M]\bar{\kappa}_{2,0}[A,B]\right.$
$\displaystyle\left.-\langle A\rangle^{2}\langle
B\rangle\left(\bar{C}_{2}[M]\left(\bar{C}_{2}[M]+C_{2}[M]\right)-\langle
M\rangle\bar{C}_{3}[M]\right)\right]$ (103)
$\displaystyle\kappa_{1,2}^{corr}[A,B]$
$\displaystyle=\kappa_{1,2}[A,B]-\frac{\langle A\rangle
C_{2}[M]\bar{\kappa}_{0,2}[A,B]}{\langle M\rangle^{2}}-\frac{2\langle B\rangle
C_{2}[M]\kappa_{1,1}[A,B]}{\langle M\rangle^{2}}+\frac{\langle
B\rangle^{2}\langle A\rangle\left(2C_{2}[M]{}^{2}-\langle M\rangle
C_{3}[M]\right)}{\langle M\rangle^{4}}$ (104) $\displaystyle\Delta_{1,2}$
$\displaystyle=\frac{1}{\langle M\rangle^{4}}\left[2\langle B\rangle\langle
M\rangle^{2}\bar{C}_{2}[M]\kappa_{1,1}[A,B]+\langle A\rangle\langle
M\rangle^{2}\bar{C}_{2}[M]\bar{\kappa}_{0,2}[A,B]\right.$
$\displaystyle\left.-\langle A\rangle\langle
B\rangle^{2}\left(\bar{C}_{2}[M]\left(\bar{C}_{2}[M]+C_{2}[M]\right)-\langle
M\rangle\bar{C}_{3}[M]\right)\right]$ (105)
$\displaystyle\kappa_{2,2}^{corr}[A,B]$
$\displaystyle=\kappa_{2,2}[A,B]+\frac{2\langle
A\rangle^{2}C_{2}[M]{}^{2}\bar{\kappa}_{0,2}[A,B]}{\langle
M\rangle^{4}}+\frac{2\langle
B\rangle^{2}C_{2}[M]{}^{2}\bar{\kappa}_{2,0}[A,B]}{\langle M\rangle^{4}}$
$\displaystyle-\frac{C_{2}[M]\bar{\kappa}_{0,2}[A,B]\bar{\kappa}_{2,0}[A,B]}{\langle
M\rangle^{2}}-\frac{\langle
A\rangle^{2}C_{3}[M]\bar{\kappa}_{0,2}[A,B]}{\langle
M\rangle^{3}}-\frac{\langle
B\rangle^{2}C_{3}[M]\bar{\kappa}_{2,0}[A,B]}{\langle M\rangle^{3}}$
$\displaystyle+\frac{12\langle A\rangle\langle B\rangle
C_{2}[M]{}^{2}\kappa_{1,1}[A,B]}{\langle
M\rangle^{4}}-\frac{2C_{2}[M]\kappa_{1,1}[A,B]{}^{2}}{\langle
M\rangle^{2}}-\frac{2\langle A\rangle C_{2}[M]\kappa_{1,2}[A,B]}{\langle
M\rangle^{2}}$ $\displaystyle-\frac{2\langle B\rangle
C_{2}[M]\kappa_{2,1}[A,B]}{\langle M\rangle^{2}}-\frac{4\langle
A\rangle\langle B\rangle C_{3}[M]\kappa_{1,1}[A,B]}{\langle
M\rangle^{3}}-\frac{10\langle A\rangle^{2}\langle
B\rangle^{2}C_{2}[M]{}^{3}}{\langle M\rangle^{6}}$
$\displaystyle+\frac{8\langle A\rangle^{2}\langle
B\rangle^{2}C_{3}[M]C_{2}[M]}{\langle M\rangle^{5}}-\frac{\langle
A\rangle^{2}\langle B\rangle^{2}C_{4}[M]}{\langle M\rangle^{4}}$ (106)
$\displaystyle\Delta_{2,2}$
$\displaystyle=\bar{C}_{2}[M]{}^{2}\left(-\frac{\langle
A\rangle^{2}\bar{\kappa}_{0,2}[A,B]}{\langle M\rangle^{4}}-\frac{\langle
B\rangle^{2}\bar{\kappa}_{2,0}[A,B]}{\langle
M\rangle^{4}}+\frac{3C_{2}[M]\langle A\rangle^{2}\langle B\rangle^{2}}{\langle
M\rangle^{6}}\right)$
$\displaystyle+\bar{C}_{2}[M]\left(-\frac{C_{2}[M]\langle
A\rangle^{2}\bar{\kappa}_{0,2}[A,B]}{\langle
M\rangle^{4}}-\frac{C_{2}[M]\langle
B\rangle^{2}\bar{\kappa}_{2,0}[A,B]}{\langle
M\rangle^{4}}+\frac{\bar{\kappa}_{0,2}[A,B]\bar{\kappa}_{2,0}[A,B]}{\langle
M\rangle^{2}}\right.$ $\displaystyle\;\;-\frac{2\bar{C}_{3}[M]\langle
A\rangle^{2}\langle B\rangle^{2}}{\langle
M\rangle^{5}}-\frac{12C_{2}[M]\langle A\rangle\langle
B\rangle\kappa_{1,1}[A,B]}{\langle M\rangle^{4}}+\frac{2\langle
A\rangle\kappa_{1,2}[A,B]}{\langle
M\rangle^{2}}+\frac{2\kappa_{1,1}[A,B]{}^{2}}{\langle M\rangle^{2}}$
$\displaystyle\left.\;+\frac{2\langle B\rangle\kappa_{2,1}[A,B]}{\langle
M\rangle^{2}}+\frac{6C_{2}[M]{}^{2}\langle A\rangle^{2}\langle
B\rangle^{2}}{\langle M\rangle^{6}}-\frac{2C_{3}[M]\langle A\rangle^{2}\langle
B\rangle^{2}}{\langle M\rangle^{5}}\right)$
$\displaystyle+\bar{C}_{3}[M]\left(\frac{\langle
A\rangle^{2}\bar{\kappa}_{0,2}[A,B]}{\langle M\rangle^{3}}+\frac{\langle
B\rangle^{2}\bar{\kappa}_{2,0}[A,B]}{\langle M\rangle^{3}}+\frac{4\langle
A\rangle\langle B\rangle\kappa_{1,1}[A,B]}{\langle
M\rangle^{3}}-\frac{4C_{2}[M]\langle A\rangle^{2}\langle B\rangle^{2}}{\langle
M\rangle^{5}}\right)$ $\displaystyle+\frac{\bar{C}_{2}[M]{}^{3}\langle
A\rangle^{2}\langle B\rangle^{2}}{\langle
M\rangle^{6}}+\frac{\bar{C}_{4}[M]\langle A\rangle^{2}\langle
B\rangle^{2}}{\langle M\rangle^{4}}$ (107)

## References

  * Rustamov  _et al._ (2023) A. Rustamov, J. Stroth, and R. Holzmann, Nucl. Phys. A 1034, 122641 (2023), arXiv:2211.14849 [nucl-th] .
  * Gross  _et al._ (2023) F. Gross _et al._ , Eur. Phys. J. C 83, 1125 (2023), arXiv:2212.11107 [hep-ph] .
  * Stephanov (2009) M. Stephanov, Phys.Rev.Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
  * Stephanov (2011) M. Stephanov, Phys.Rev.Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Bzdak  _et al._ (2020) A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rept. 853, 1 (2020), arXiv:1906.00936 [nucl-th] .
  * Asakawa  _et al._ (2009) M. Asakawa, S. Ejiri, and M. Kitazawa, Phys. Rev. Lett. 103, 262301 (2009), arXiv:0904.2089 [nucl-th] .
  * Bazavov  _et al._ (2017) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 96, 074510 (2017), arXiv:1708.04897 [hep-lat] .
  * Borsanyi  _et al._ (2018) S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, and C. Ratti, JHEP 10, 205 (2018), arXiv:1805.04445 [hep-lat] .
  * Friman  _et al._ (2011) B. Friman, F. Karsch, K. Redlich, and V. Skokov, Eur.Phys.J. C71, 1694 (2011), arXiv:1103.3511 [hep-ph] .
  * Abdallah  _et al._ (2021) M. Abdallah _et al._ (STAR), Phys. Rev. C 104, 024902 (2021), arXiv:2101.12413 [nucl-ex] .
  * Abdallah  _et al._ (2023) M. Abdallah _et al._ (STAR), Phys. Rev. C 107, 024908 (2023), arXiv:2209.11940 [nucl-ex] .
  * Adamczewski-Musch  _et al._ (2020) J. Adamczewski-Musch _et al._ (HADES), Phys. Rev. (2020), arXiv:2002.08701 [nucl-ex] .
  * Acharya  _et al._ (2020) S. Acharya _et al._ (ALICE), Phys. Lett. B 807, 135564 (2020), arXiv:1910.14396 [nucl-ex] .
  * Acharya  _et al._ (2023) S. Acharya _et al._ (ALICE), Phys. Lett. B 844, 137545 (2023), arXiv:2206.03343 [nucl-ex] .
  * Fu  _et al._ (2016) W.-j. Fu, J. M. Pawlowski, F. Rennecke, and B.-J. Schaefer, Phys. Rev. D 94, 116020 (2016), arXiv:1608.04302 [hep-ph] .
  * Bleicher  _et al._ (2000) M. Bleicher, S. Jeon, and V. Koch, Phys. Rev. C62, 061902 (2000), arXiv:hep-ph/0006201 [hep-ph] .
  * Bzdak  _et al._ (2013) A. Bzdak, V. Koch, and V. Skokov, Phys. Rev. C87, 014901 (2013), arXiv:1203.4529 [hep-ph] .
  * Braun-Munzinger  _et al._ (2019) P. Braun-Munzinger, A. Rustamov, and J. Stachel, Nucl. Phys. A 982, 307 (2019), arXiv:1807.08927 [nucl-th] .
  * Braun-Munzinger  _et al._ (2021) P. Braun-Munzinger, B. Friman, K. Redlich, A. Rustamov, and J. Stachel, Nucl. Phys. A 1008, 122141 (2021), arXiv:2007.02463 [nucl-th] .
  * Braun-Munzinger  _et al._ (2023) P. Braun-Munzinger, K. Redlich, A. Rustamov, and J. Stachel, (2023), arXiv:2312.15534 [nucl-th] .
  * Kitazawa and Asakawa (2012) M. Kitazawa and M. Asakawa, Phys. Rev. C85, 021901 (2012), arXiv:1107.2755 [nucl-th] .
  * Skokov  _et al._ (2013) V. Skokov, B. Friman, and K. Redlich, Phys. Rev. C88, 034911 (2013), arXiv:1205.4756 [hep-ph] .
  * Braun-Munzinger  _et al._ (2017) P. Braun-Munzinger, A. Rustamov, and J. Stachel, Nucl. Phys. A960, 114 (2017), arXiv:1612.00702 [nucl-th] .
  * Bialas  _et al._ (1976) A. Bialas, M. Bleszynski, and W. Czyz, Nucl. Phys. B 111, 461 (1976).
  * (25) https://github.com/anarrustamov/VolumeFluctuations.
  * Adamczewski-Musch  _et al._ (2018) J. Adamczewski-Musch _et al._ (HADES), Eur. Phys. J. A 54, 85 (2018), arXiv:1712.07993 [nucl-ex] .
  * Abelev  _et al._ (2013a) B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044909 (2013a), arXiv:1301.4361 [nucl-ex] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR), Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Loizides  _et al._ (2015) C. Loizides, J. Nagle, and P. Steinberg, SoftwareX 1-2, 13 (2015), arXiv:1408.2549 [nucl-ex] .
  * Harabasz  _et al._ (2023) S. Harabasz, J. Kołaś, R. Ryblewski, W. Florkowski, T. Galatyuk, M. Gumberidze, P. Salabura, J. Stroth, and H. P. Zbroszczyk, Phys. Rev. C 107, 034917 (2023), arXiv:2210.07694 [nucl-th] .
  * Abelev  _et al._ (2013b) B. Abelev _et al._ (ALICE), Phys. Rev. C 88, 044910 (2013b), arXiv:1303.0737 [hep-ex] .
  * Bzdak  _et al._ (2017) A. Bzdak, V. Koch, and V. Skokov, Eur. Phys. J. C77, 288 (2017), arXiv:1612.05128 [nucl-th] .


# Quarkophobic W’ for LHC searches

A. Gurrola1, J. D. Ruiz-Álvarez2 1Department of Physics and Astronomy,
Vanderbilt University, Nashville, TN, 37235, USA
2Instituto de Física, Universidad de Antioquia, A.A. 1226 Medellín, Colombia.

###### Abstract

We consider a simplified model where a W’ boson is added to the standard model
with negligible couplings to quarks, but generic couplings to leptons and
electroweak bosons. We study the implications of such a model for LHC
searches. Consequently, we propose an LHC search through the vector boson
fusion topology which would have sensitivity for such a new particle with the
current proton-proton collisions’s energy and available luminosity.

##  I Introduction

In recent years, there have been several measurements that have reported
deviations form the predictions of the standard model (SM) of particle
physics. Experiments with different setups and experimental configurations
have measured some B mesons decays with branching fractions that challenge SM
theoretical predictions BaBar:2006tnv ; BaBar:2015wkg ; Belle:2009zue ;
Belle:2016xuo ; Belle:2019oag ; LHCb:2019hip ; LHCb:2015wdu ; LHCb:2015wdu ;
LHCb:2014cxe . However these measurement results have been less than 5 sigma,
they have awaken the interest of the particle physics community. A plethora of
models beyond the standard model have been proposed in order to explain such
anomalies.

It has been proved recently that these low energy anomalies should have
connections with high energy physics as for proton-proton collisions at the
Large Hadron Collider Greljo:2018tzh ; Fuentes-Martin:2020lea . The most
traditional model building for b-anomalies consider a new physics mediator
that couples to quarks and leptons of the SM. Fig. 1 displays the Feynman
diagram for a W’ that couples to both leptons and quarks from the SM.

Figure 1: Feynman diagram for a W’ mediating a process involved in b-anomalies
where the W’ couples to SM quarks.

However, one can consider a W’ with negligible or no couplings to the SM
quarks (quarkophobic). This W’ can couple to the electroweak bosons through
triple gauge couplings (TGC) which could also mediate the b-anomalies
processes. The corresponding Feynman diagram for such a model is depicted in
Fig. 2

Figure 2: Feynman diagram for a quarkophobic W’ mediating a process involved
in b-anomalies through TGC.

A quarkophobic W’ produces naturally a vector boson fusion (VBF) topology at
the LHC proton-proton collisions. This fact opens the door for searching for
this new particle at high energies, and additionally motivates a search for a
single high energy lepton, from the W’ decay, missing transverse momentum and
two VBF jets. In Fig. 3 is shown the VBF production of a quarkophobic W’.

Figure 3: Feynman diagram for the production of a quarkophobic W’ from proton-
proton collisions

We report, in a very summarized manner, a proposal for a search at the LHC of
a quarkophobic W’ produced through a VBF topology in proton-proton collisions
at the LHC.

##  II A minimal model

We consider to include to the SM content a new W’ particle with the following
minimal couplings to SM particles, as shown in Eq. 1\- 3.

$\displaystyle\mathcal{L}_{VWW^{\prime}}^{1}$ $\displaystyle=$ $\displaystyle
g_{1}^{V}V^{\mu}(W^{-}_{\mu\nu}W^{\prime+\nu}-W^{+}_{\mu\nu}W^{\prime-\nu}$
(1)
$\displaystyle+W^{\prime-}_{\mu\nu}W^{+\nu}-W^{\prime+}_{\mu\nu}W^{-\nu}),$

$\mathcal{L}_{VWW^{\prime}}^{2}=g_{2}^{V}(W^{+}_{\mu}W^{\prime-}_{\nu}V^{\mu\nu}+W^{\prime+}_{\mu}W^{-}_{\nu}V^{\mu\nu}),$
(2)

where $V_{\mu\nu}=\partial_{\mu}V_{\nu}-\partial_{\nu}V_{\mu}$ and
$V=Z\;\text{or}\;\gamma$. The coefficients $g_{1,2}^{V}$ represent the
coupling strength between the $W^{\prime}$ and SM weak bosons, and thus govern
the VBF production cross section. The interactions between the $W^{\prime}$
and the SM leptons are described by

$\mathcal{L}_{l}=\sum_{l}\bar{\nu}_{l}\gamma_{\mu}(g^{R}_{l}(1+\gamma^{5})+g^{L}_{l}(1-\gamma^{5}))W^{\prime\mu}l,$
(3)

where the coefficients $g_{l}^{L,R}$ are the left-handed and right-handed
couplings that govern the $W^{\prime}\to l\nu$ decays. For the purpose of the
studies shown in this paper, we set the values of the couplings to
$\\{g_{l}^{L},g_{l}^{R},g^{V}\\}=\\{1,1,1\\}$.

To have a complete model with a W’ with similar couplings as the ones
described, a serious amount of theoretical work must be done. However, for our
purpose we only need a model that allow us to simulate the interactions of
interest to have a W’ produced via a VBF topology at the LHC with subsequent
decays to SM leptons. In this sense, the model we propose is a simplified W’
model, as the dark matter community has been doing, for LHC searches with the
specificity of having TGC and negligible couplings of the W’ to quarks.

##  III Samples and Strategy

The signal model was implemented with FeynRules Alloul:2013bka  and used in
its UFO Degrande:2011ua  implementation for the production of montecarlo
samples. The partonic processes for signal and background events have been
produced with MadGraph5_aMC (v2.8.2) Alwall:2014hca . The showering and
hadronization have been performed by PYTHIA8 Bierlich:2022pfr , and the
detector simulation has been done by Delphes (v3.4.2) deFavereau:2013fsa .

The strategy for the search is based on the characteristics of the jets coming
from the VBF topology. These are, having a high pseudorapidity separation,
high invariant mass and being produced in separate hemispheres of a detector
like CMS or ATLAS. Additionally we have studied the kinematics of the leptons
from the decay of the W’ in order to further suppress SM backgrounds.

For the optimization of the selection criteria we use a signal sample with a
W’ with a mass of 1 TeV. As our signature is one lepton, missing transverse
momentum and two jets we have considered as our main backgrounds W+jets,
$t\bar{t}$ and diboson SM production. We use the significance $S/\sqrt{S+B}$
as the figure to maximize for the optimization of the selection criteria. $S$
stands for the number of signal events while $B$ represents the sum of events
from all the background processes. For the proposed optimization, a dummy
cross section of 0.1 pb has been used for the W’. Real cross sections can be
easily obtained from MadGraph accordingly to specific values of the couplings
and W’ mass.

##  IV Results

We describe in the following the results for the case where the W’ decays into
a muon. We have also studied the electron and tau channels, and these results
will be put together in a complete separate paper.

In Fig. 4 are shown the characteristics of the muon and the missing transverse
momentum that were found to be useful to discriminate signal from backgrounds.
Additionally, the transverse mass distribution for the muon is shown. For a
complete search this last distribution could be considered as the variable of
interest.

Figure 4: Signal and backgrounds normalized to unit distributions for (top-
left) $p_{T}$ of the muon, (top-right) $p_{T}^{miss}$, (bottom-left) the
absolute value of the $\Delta\phi$ betwwen the muon and the missing transverse
momentum and (bottom-right) the transverse mass of the muon and the missing
transverse momentum

The selection criteria developed is:

  1. 1.

A baseline selection consisting of requiring at least two jets with $p_{T}>60$
GeV, $N(b)=0$, $N(\mu)=1$, and veto leptons with $p_{T}(l)>25$ GeV

  2. 2.

$\eta(j_{1})\times\eta(j_{2})<0$

  3. 3.

$m_{jj}>1000$ GeV

  4. 4.

$|\Delta(\eta(j_{1}),\eta(j_{2}))|>4.0$

  5. 5.

$p_{T}(\mu)>200$ GeV

  6. 6.

$p_{T}^{miss}>200$ GeV

  7. 7.

$|\Delta(\phi(\mu),p_{T}^{miss})|>1.0$

In Table 1 are shown the event yields for the signal with $m(W^{\prime})=1$
TeV and the SM backgrounds at each step of the selection. Also is shown the
significance as previously defined for each selection criteria. The number of
events are calculated for proton-proton collisions at the LHC at a centrer of
mass energy of 13 TeV and a luminosity of 150 $\text{fb}^{-1}$.

Table 1: Event yields for muon channel signal and backgrounds at each selection step. | Signal | W+jets | $t\bar{t}$ | VV | $S/\sqrt{S+B}$
---|---|---|---|---|---
Initial | 1259.4 | 152156.8 | 48246.6 | 1870.3 | 2.8
Cut 1 | 930.3 | 68480.3 | 14364.6 | 788.2 | 3.2
Cut 2 | 593.7 | 13461.8 | 1796.9 | 147.8 | 4.7
Cut 3 | 527.1 | 2771.8 | 365.0 | 27.9 | 8.7
Cut 4 | 474.0 | 118.4 | 17.6 | 2.1 | 19.2
Cut 5 | 459.9 | 91.7 | 11.6 | 1.9 | 19.3
Cut 6 | 459.6 | 0.0 | 0.0 | 0.9 | 21.4

##  V Conclusions

We have proposed a model for searches at the LHC looking for a W’ produced
through the VBF topology. We have also developed a search that is capable of
efficiently looking for such a W’ in the context of the data collected at the
LHC for the run 2. We have shown the details for this search when the W’
decays into a muon a neutrino but we will be producing a complete paper were
also the tau and electron channels are included. Our finding strongly
motivates for searching for a heavy quarkophobic W’ with the current data at
the LHC.

###### Acknowledgements.

A.G. gratefully acknowledge the constant support of the Physics & Astronomy
department at Vanderbilt University and the US National Science Foundation.
This work is supported in part by NSF Award PHY-1945366. J. D. R.-Á.
gratefully acknowledges the support of Universidad de Antioquia, the Colombian
Science Ministry and Sostenibilidad-UdeA.

## References

  * (1) B. Aubert et al. [BaBar], Phys. Rev. D 73, 092001 (2006) doi:10.1103/PhysRevD.73.092001 [arXiv:hep-ex/0604007 [hep-ex]].
  * (2) J. P. Lees et al. [BaBar], Phys. Rev. D 93, no.5, 052015 (2016) doi:10.1103/PhysRevD.93.052015 [arXiv:1508.07960 [hep-ex]].
  * (3) J. T. Wei et al. [Belle], Phys. Rev. Lett. 103, 171801 (2009) doi:10.1103/PhysRevLett.103.171801 [arXiv:0904.0770 [hep-ex]].
  * (4) A. Abdesselam et al. [Belle], [arXiv:1604.04042 [hep-ex]].
  * (5) A. Abdesselam et al. [Belle], Phys. Rev. Lett. 126, no.16, 161801 (2021) doi:10.1103/PhysRevLett.126.161801 [arXiv:1904.02440 [hep-ex]].
  * (6) R. Aaij et al. [LHCb], Phys. Rev. Lett. 122, no.19, 191801 (2019) doi:10.1103/PhysRevLett.122.191801 [arXiv:1903.09252 [hep-ex]].
  * (7) R. Aaij et al. [LHCb], JHEP 09, 179 (2015) doi:10.1007/JHEP09(2015)179 [arXiv:1506.08777 [hep-ex]].
  * (8) R. Aaij et al. [LHCb], JHEP 09, 179 (2015) doi:10.1007/JHEP09(2015)179 [arXiv:1506.08777 [hep-ex]].
  * (9) R. Aaij et al. [LHCb], JHEP 06, 133 (2014) doi:10.1007/JHEP06(2014)133 [arXiv:1403.8044 [hep-ex]].
  * (10) A. Greljo, J. Martin Camalich and J. D. Ruiz-Álvarez, Phys. Rev. Lett. 122, no.13, 131803 (2019) doi:10.1103/PhysRevLett.122.131803 [arXiv:1811.07920 [hep-ph]].
  * (11) J. Fuentes-Martin, A. Greljo, J. Martin Camalich and J. D. Ruiz-Alvarez, JHEP 11, 080 (2020) doi:10.1007/JHEP11(2020)080 [arXiv:2003.12421 [hep-ph]].
  * (12) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. 185, 2250-2300 (2014) doi:10.1016/j.cpc.2014.04.012 [arXiv:1310.1921 [hep-ph]].
  * (13) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, Comput. Phys. Commun. 183, 1201-1214 (2012) doi:10.1016/j.cpc.2012.01.022 [arXiv:1108.2040 [hep-ph]].
  * (14) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * (15) C. Bierlich, S. Chakraborty, N. Desai, L. Gellersen, I. Helenius, P. Ilten, L. Lönnblad, S. Mrenna, S. Prestel and C. T. Preuss, et al. [arXiv:2203.11601 [hep-ph]].
  * (16) J. de Favereau et al. [DELPHES 3], JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].