Centrality fluctuations and decorrelations in heavy-ion collisions

The centrality or the number of initial-state sources $V$ of the system produced in heavy ion collision is a concept that is not uniquely defined and subject to significant theoretical and experimental uncertainties. We argue that a more robust connection between the initial-state sources with final-state multiplicity could be established from the event-by-event multiplicity correlation between two subevents separated in pseudorapidity, $N_a$ vs $N_b$. This correlation is sensitive to two main types of centrality fluctuations (CF): 1) particle production for each source $p(n)$ which smears the relation between $V$ and $N_a$ used for experimental centrality, and 2) decorrelations between the sources in the two subevents $V_b$ and $V_a$. The CF is analyzed in terms of cumulants of $V_b$ and $N_b$ as a function of $N_a$, i.e. experimental centrality is defined with $N_a$. We found that the mean values $\langle V_b\rangle_{N_a}$ and $\langle N_b\rangle_{N_a}$ increase linearly with $N_a$ in mid-central collisions, but flatten out in ultra-central collisions. Such non-linear behavior is sensitive to the centrality resolution of $N_a$. In the presence of centrality decorrelations, the scaled variances $\langle(\delta V_b)^2\rangle/\langle V_b\rangle$ and $\langle(\delta N_b)^2\rangle/\langle N_b\rangle$ are found to decrease linearly with $N_a$ in mid-central collisions, while the $p(n)$ leads to another sharp decrease in the ultra-central region. The higher-order cumulants of $V_b$ and $N_b$ show interesting but rather complex behaviors which deserve further studies. Our results suggest that one can use the cumulants of the two-dimensional multiplicity correlation, especially the mean and variance, to constrain the particle production mechanism as well as the longitudinal fluctuations of the initial-state sources.


I. INTRODUCTION
In heavy-ion collisions, centrality is an important concept for characterizing the size of the produced fireball. The centrality of an event is often represented by the number of particle production sources V (for example participating nucleons or constituent quarks) in the initial state. These sources are then used to define other important geometry quantities, such as the nuclear overlap function [1,2] for jet quenching studies [3,4] and eccentricities [5,6] for collective flow studies. [7,8] Since the sources in an event are not directly measurable, a Glauber model is used to calculate V and relate it to the final observed particle multiplicity N , with N = ∑ V i=1 n i calculated as a sum of the multiplicity for each source n i sampled independently from a common distribution p(n) [9][10][11]. Due to fluctuations in the particle production, events selected with the same V can have different values of N and vice versa. The fluctuation of sources at fixed multiplicity, often referred to as the "volume fluctuations", is an irreducible "centrality fluctuations" (CF) or centrality resolution [12][13][14][15]. The main issue concerning centrality is to understand how the final-state p(N ) is related to the source distribution p(V ) and particle production for each source p(n).
The centrality or volume of an event is often assumed to be a global concept, long-range in pseudorapidity (η): a central event should be a central event independent of η. This assumption can be checked by studying the the correlation of multiplicities between two different η windows, also referred to as forward-backward (FB) multiplicity correlation [16][17][18][19][20]. In experimental analysis, the observed centrality is usually defined as the particle multiplicity N a in a forward pseudorapidity window A, and the measurement is performed using particles with total multiplicity N b in a different pseudorapidity window B, usually around mid-rapidity. If the system is boost invariant and same sources are responsible for particles in both subevents (see Fig 1(a)), the particle multiplicity distribution in subevent B for fixed N a can be expressed as, The distribution p(V ) Na , reflecting the extent of the CF, is controlled by the smearing from p(n) in subevent A. Narrow p(V ) Na distribution would imply poor centrality resolution in subevent A and vice versa. Eq. (1) shows that the CF arising from subevent A directly influences the multiplicity fluctuations in subevent B, and this is the picture assumed by most previous analyses [9][10][11].
Recently, several studies also consider the possibility that the number of sources and their transverse distributions may fluctuate along η in a single event. This is because the number of forward-going and backward-going participating nucleons, i.e. N F part and N B part , are not the same in a given event [8,20], and they contribute differently in the forward FIG. 1: Schematic illustration of the relation between sources (large circles) and final-state particles (small circles) produced in subevent A for centrality determination and subevent B for measurements for two cases: 1) same sources, wounded nucleons Npart, for both subevents (panel a), and 2) different sources, forward-going and backward-going wound nucleons N F part and N B part , for the two subevents (panel b). The second scenario is particularly relavent for fixed-target experiments, which often use the forward multiplicity for centrality and mid-or backward rapidity for multiplicity measurements.
(backward) rapidity. This forward-backward asymmetry is realized in one of the two ways in dynamical models: 1) In string picture [21][22][23], each string originates from a participant nucleon and extends to the opposite direction. Therefore, the number of strings, their lengths and end-points fluctuate in rapidity. 2) In a gluon saturation picture, the sub-nucleonic degree-of-freedoms evolve with rapidity [24]. In the forward rapidity, the projectile nucleons are dominated by a few large-x partons, while the target nucleons are expected to contribute mainly low-x partons. In this case, the Eq. (1) can be revised as, The fluctuation between the number of sources in the two subevents, denoted by p(V b ) Va , can weaken the centrality correlation between different rapidities (see Fig 1(b)). Such "centrality decorrelations" effects is expected to have a strong influence on the correlation of multiplicities between two different η windows. Eqs. (1) and (2) show that information about centrality fluctuations and particle productions is encoded in a set of one-dimensional distributions, which are customarily analyzed via multiplicity cumulants of p(x) (x = N , n or V ): meanx ≡ ⟨x⟩, variance ⟨(x −x) 2 ⟩, skewness ⟨(x −x) 3 ⟩,... In this paper, we use the reduced form to scale out the overall system size: Previous studies mainly focused on the scaling behavior of the mean multiplicity per source, N ⟨V ⟩ [11,[25][26][27][28][29]. Our paper extends these studies to higher-order moments/cumulants of p(N ) and p(V ), and their dependences on particle production p(n). These cumulants reveal detailed structures of the two-dimensional correlation between N b and N a , and provide constraints on p(n), p(V a ) and p(V b ) Va . A good understanding of CF is also important for other areas of heavy-ion physics: for any physics observable x that varies with V , its fluctuation p(x) should be sensitive to p(V ) [15]. The CF is one of the main background sources in the ongoing search for the critical end point in the QCD phase diagram based on the fluctuations of conserved quantities [30,31]. The CF also strongly affects the fluctuations of harmonic flow [15], and is responsible for the observed sign-change behavior in ultra-central A+A collisions (UCC) [32]. Understanding the longitudinal fluctuations of sources and the resulting centrality decorrelation effects can also help to describe previous measurements of harmonic flow decorrelations in η [22,[33][34][35].
In this paper, we study the influence of centrality smearing from particle production (Fig 1a)) and centrality decorrelations from the initial sources (Fig 1b)) on the event-by-event forward-backward multiplicity fluctuations. The study is based on an independent source Glauber model with particle production tuned to reproduce the ATLAS data [32], and the initial number of sources is chosen to be the V = N part for the scenario shown in Fig 1a) and V a ≡ N F part and V b ≡ N B part for the scenario shown in Fig 1b). We found that the fluctuations in the UCC region are very sensitive to the source distribution p(V ) and particle production p(n), due to the steeply falling distributions (upper boundary effect) [13,36,37]. We also found that the multiplicity dependence of the mean and variance of p(V ) have a very simple interpretation in terms of the two types of centrality fluctuations, and they can constrain the longitudinal structure of initial sources and particle production mechanism via comparison between data and model. The structure of the paper is the following. Section II and III discuss the expected general behavior of the CF in the presence of multiplicity smearing and centrality decorrelations, respectively. After a brief description of the Glauber model setup in Section IV, the results for the two types of CF are given in Section V and VI, respectively.
The summary of the results are given in Section VII. The Appendix A gives a derivation of the formulas used in this paper, and some additional results are given in Appendix B.

II. CENTRALITY FLUCTUATIONS DUE TO MULTIPLICITY SMEARING
The CF arising from multiplicity smearing effects (Fig. 1a)) is considered in an independent source model framework, by assuming that the total particle multiplicity is a sum of particles from each source N ≡ ∑ V i=1 n i sampled from p(n). In this case, one can show that cumulants of p(N b ) V in Eq. (1) is related to p(n) by a factor V [13], i.e.
From this, Skokov et al. [13] find the following relations connecting the cumulants for the three distributions in Eq. (1): where we have used a simpler notation, The behavior of centrality cumulants k n,V Na for p(V ) Na in Eq. (1) can be understood qualitatively from the correlation between V and N a in Fig. 2. In the absence of smearing from particle production, i.e. by setting p(n a ) = δ(n a −n a ), V and N a are related by a linear function N a = Vn a (thick solid line). In the presence of smearing, p(n a ) broadens the correlation along the x-axis (shaded area). As a result, the ⟨V ⟩ Na deviates from the linear relation in the UCC region (thick dashed line). This deviation arises because the distribution p(V ) Na is bounded above by V max (twice of the atomic number for N part ), leading to an asymmetric shape in the UCC region as illustrated by the middle and right insert panels. For even larger N a values, the p(V ) Na becomes narrower and more asymmetric, and eventually converges to V max . For this reason, the cumulants of p(V ) as a function of N a are strongly modified in the UCC collisions: in addition to the departure of the mean from the linear relation, the variance of p(V ) Na is expected to approach zero, the skewness of p(V ) Na becomes negative and then approaches zero, and the kurtosis of p(V ) Na oscillates from negative to positive then approaches zero. Outside the UCC region, where p(V ) Na is not constrained by the upper boundary effect (left insert panel), the ⟨V ⟩ Na is expected to overlap with the solid line, and the higher-order cumulants of p(V ) Na are determined mainly by p(n a ). In this paper, we focus on the first four cumulants of CF, as their importance can be clearly identified from the shape of p(V ) Na .
The origin of CF and the behavior of centrality cumulants k n,V Na for n = 2-4 were studied in detail in Ref. [15]. In the UCC region, the moments of CF (and therefore k n,V Na ) are driven mainly by the shape of p(V ), where we assume p(N a ) V ≈ is the relative width for p(n a ). In the mid-central region where p(V ) is slowly changing, k n,V Na are found to be approximately constant, and the scaled variance of the multiplicity fluctuation in subevent B is sensitive to both p(n a ) and p(n b ), This relation implies that for mid-central collisions and within the independence source picture, the variance of the measured multiplicity fluctuations has the same functional dependence on p(n a ) and p(n b ) in the two subevents.

III. CENTRALITY DECORRELATIONS DUE TO LONGITUDINAL FLUCTUATIONS
Next we consider the CF due to centrality decorrelations as illustrated by Fig. 1b). Equation (2) can be decomposed into two equations The cumulants for the first equation are given by Eq. (5). The second equation contains the effects of centrality decorrelations from p(V b ) Va , which, unlike the p(N b ) V b , is not described by the independent source picture (i.e V b is not an independent sum of V a number of sources via some common distribution). Indeed, we find that the cumulants of p(V b ) Va are not strictly proportional to V a . However, by assuming the departure from such proportionality is small, we have derived the relations including the first-order correction (see Appendix A): The k n,ab are the reduced cumulants for p(V b ) Va , while k ′ n,ab also include the leading-order correction, k ′ n,ab = ∂(k 1,ab k n,ab V a ) k 1,ab ∂V a = k n,ab + V a ∂(k 1,ab k n,ab ) k 1,ab ∂V a , n ≥ 2 The total multiplicity cumulants in subevent B, K n,b , can be obtained by replacing the k n,V Na in Eq. (5) by Eq. (10). If the centrality resolution of subevent A is very good, k n,Va Na ≈ 0 and k n,V b ≈ k n,ab for n > 1, and we recovers the Eq. (5). In this case, the behavior of K n,b is dictated by the FB fluctuations.

IV. GLAUBER MODEL SETUP
For quantitative study of the effects of multiplicity smearing and the longitudinal centrality decorrelations, we follow the implementation of our previous work [15], which is described briefly here. The number of sources in subevent A and B are chosen to be N F part and N B part calculated for each event in a standard Glauber model framework [10]. The nucleons are assumed to have a hard-core of 0.3 fm in radii, their transverse positions are generated according to the Woods-Saxon distribution as provided by Ref. [38]. A nucleon-nucleon cross-section of σ = 68 mb is used to simulate the collisions at √ s NN = 5.02 TeV. Since the p(V ) distribution as well as the correlation between N F part and N B part change with system size, we studied several nuclei covering a broad range of the V : Xe+Xe, Cu+Cu, S+S and O+O collisions, with total number of nucleons 2A = 258, 126, 64, 32 respectively. The particle productions for each source are chosen to follow the negative binomial distribution (NBD): wheren is the average number of particles in a given subevent. One important quantity is the relative widthσ: σ 2 ≡ ⟨(δn) 2 ⟩ n 2 , which controls the strength of fluctuations for each source. Table I lists the three NBD parameter sets used to produce particles in subevents A and B, taken directly from Ref. [15]. The Par0 and Par1 are adjusted to approximately describe the shapes of the experimental distributions of N rec ch ( η < 2.5) and ΣE T (3.2 < η < 4.9) from the ATLAS Collaboration [32], while the Par2 corresponds to a case with much larger fluctuations.  (12)) used for modeling the particle production in the wounded nucleon model from Ref. [15]. The corresponding values of reduced cumulants k2, k3 and k4 are also listed. The calculation of cumulants follows the standard procedure. Each A+A event has two subevents: subevent A for centrality selection and subevent B for the calculation of multiplicity cumulants. The particle multiplicities in these two subevents, N a and N b , are generated independently from 1) the same sources V a = V b = N part for the study of multiplicity smearing effect or 2) V a = N F part and V b = N B part for the study of longitudinal centrality decorrelations. The generated events are divided into classes according to N a . The cumulants are first calculated from V and N b distributions for events with the same N a , which are then combined into broader N a ranges to reduce the statistical uncertainty. In the following, we first discuss the behavior of cumulants by assuming V a = V b = N part , we then show results for cumulants calculated by assuming V a = N F part and V b = N B part .

V. RESULTS ON CENTRALITY FLUCTUATIONS FROM PARTICLE PRODUCTION
We first consider the case when participant nucleons are used as common sources for subevents A and B, Figure 3 illustrates the behavior of the centrality fluctuations by selecting events with fixed N a in Pb+Pb collisions. The top three panels show the correlation between V and N a generated with Par0, Par1 and Par2 in Table I. The centrality cumulants k n,V for distribution p(V ) Na are calculated and shown in the bottom panels.
The behavior of these cumulants follows the naive expectation of Fig. 2. The k 1,V = ⟨V ⟩ is proportional to N a , except in UCC region where it turns over. For largerσ value, this turn-over starts earlier and extends to larger N a range, as expected from a poorer centrality resolution. The insert panel shows the ratio N a (n ⟨V ⟩), to quantify the deviation of particle per source fromn. This ratio should be unity in the absence of centrality smearing effects. The ratio exhibits a suppression in the peripheral region and an enhancement in the central region, as expected from the boundary effects. We notice that in the central region, the ratios for different parameter sets are nearly parallel to each other. We found that the slopes of the ratios are mainly controlled by the shape of the p(N part ) distributions for  Table I. Bottom row: the corresponding cumulants of V , k n,V , calculated for the three parameter sets, with n = 1, 2, 3, and 4 from left to right panels. The insert panel shows the ratio of Na (n ⟨V ⟩), to quantify the deviation of particle per source fromn (or from the linear relation).
the UCC events. The largerσ value only shifts the turn-over point to a smaller N a value and extends the suppression region to a larger N a value, but has little impact on the slope. The behavior of higher-order centrality cumulants follows what was observed before in Ref. [15]. The k n,V are flat in mid-central collisions but show strong variations in central region: the k 2,V decreases to 0, the k 3,V decreases to negative minimum values and then approaches 0 from below, while the k 4,V decreases to a negative minimum, increases to reach a positive maximum and then decreases to approach 0. The larger smearing associated with Par2 further enhances the magnitudes of k n,V and their oscillating behaviors in the ultra-central region.
The source distributions for events with fixed N a , p(V ) Na , is then used to produce the distribution of N b according to a common distribution p(n b ) for each source. The top panels of Fig. 4 shows the correlation between N b and N a , obtained by smearing the 2D distributions in Fig. 3 along the y-axis using the Par1 for p(n b ). The cumulants of p(N b ) Na , K n,b , are calculated and shown in the bottom panels. According to Eq. (5), they are expressed as a linear combination of k n, where the coefficients are determined from the k n for Par1 in Table I. Eq. (13) shows that the ⟨N b ⟩ is related to ⟨V ⟩ in Fig. 3 by a constant scale factor. The K 2,b is related to k 2,V by a constant offset and a constant scale factor, etc. In general, one could isolate the k n,V iteratively order-by-order: k 2,V can be extracted from K 2,b by identifying and subtracting a constant, which then can be subtracted from K 3,b to isolate the k 3,V , etc. Note that this is only true in the independent source model framework and considering only multiplicity smearing effects. The centrality fluctuations are sensitive to the shape of p(V ≡ N part ), which changes with the size of the collision system. The top-left panel of Fig. 5 shows the distribution p(V ) from five collision systems in central collision region. The same distribution is replotted in the bottom-left panel, but the x-axes have been rescaled by V max = 2A.
The smaller system shows a broader tail in the V distributions. For system-size dependence studies, the centrality cumulants k n,V are calculated with Par1 from Table I. The results are presented both as a function of N a and as a function of N a N knee , where the knee is defined as the average multiplicity for maximum V = 2A nucleons, N knee = 2An. The 2nd column of Fig. 5 show the k 1,V = ⟨V ⟩ as a function of N a and N a N knee . The ⟨V ⟩ increases almost linearly with N a , and flattens out close to N knee . To quantify the non-linear behavior in the UCC region, we calculate the ratios N a (n a ⟨V ⟩) and present them in the third column. In the very small N a region, the ratios are   Table I. below unity and agree with each other, reflecting the fact that the p(V ) distributions have very similar shape in the small V region. Towards the large N a region, the ratios separate from each other and increase to above unity with the smaller system showing a larger deviation from unity. The third bottom panel shows the ratios as a function of N a N knee , where the increase in UCC is linear for all systems. For a smaller system, the rate of increase is smaller but over a broader range in N a N knee , consistent with the smoother fall-off for smaller system in the central region  Table I. shown of the bottom-left panel. Previous studies on the particle production often focused on the scaling behavior of the particle production per source, i.e. N ⟨V ⟩, with V defined as N part or the number of quark participants. They found that the sources based on the quark participant number has a better scaling behavior than that based on the nucleon participant number from pp, p+A to A+A systems [11,[25][26][27][28]. However, a detailed study by ALICE Collaboration showed that both types sources show a sharp increase of N ⟨V ⟩ vs ⟨V ⟩ in the UCC region [28]. To perform a similar check, the right column of Fig. 5 shows N a (n a ⟨V ⟩) as a function of ⟨V ⟩ and ⟨V ⟩ 2A, where ⟨V ⟩ is obtained by mapping from N a using the data in the top panel in the second column. The linear increase in the third column now appears as a sharp increase in the UCC region, similar to the data. The increase is stronger and affects most of the centrality range in a small system, consistent with the experimental observation.
On the other hand, the experimental observable is defined as N b ⟨V ⟩ as appose to N a ⟨V ⟩ in our study. In the independent source model framework, N b ⟨V ⟩ would always equal ton by construction. Therefore, one natural explanation would be that there are strong correlations between the particle production in subevents A and B, i.e. n a and n b are strongly correlated. In an extreme case of n a = n b , the N b = N a would be valid for each event, and N b ⟨V ⟩ would also show an increase in the UCC region, we leave this point to a future investigation. Figure 6 shows higher-order k n,V as a function of N a (top row) or N a N knee (bottom row). The k n,V have the same constant values in mid-central collisions, but deviates from each other towards more peripheral or more central region where the lower and upper boundary effects from the p(V ) distribution are important. Such deviation appears over a larger fraction of the N a range for smaller collision systems. In the central region, the magnitudes of the k 3,V and k 4,V are smaller for smaller collision systems. At the same time, the widths of the maxima region are larger in terms of N a N knee , implying that the centrality fluctuations influences a larger fraction of the centrality region. This can be understood from the bottom-left panel in Fig. 5 , which shows that the decrease in the central region is smoother in smaller collision systems.

VI. RESULTS ON FORWARD-BACKWARD CENTRALITY DECORRELATIONS
To study the effects of centrality decorrelations, we consider V a = N F part and V b = N B part as the sources for particles in subevent A and B, respectively. According to Eq. (10), the cumulants of p(V b ) are expressed in terms of cumulants k n,ab and k ′ n,ab describing p(V b ) Va , as well as k n,Va describing the fluctuation of V a . In the following, we first discuss the behaviors of k n,Va , k n,ab and k ′ n,ab and then show quantitatively how well k n,V b are described by Eq. (10). Figure 7 shows the centrality cumulants k n,Va for Pb+Pb and O+O collisions. The magnitudes of k n,Va depend on the centrality resolution of subevent A, and are calculated with Par0-Par2 as a function of N a . The results for Pb+Pb are similar to those presented in the bottom row of Fig. 3. The only difference is that the N F part instead of N part is used as sources in this figure, therefore the maximum range of N a reaches only about half of that in Fig. 3. The results for O+O are much more affected by the boundary effects, such that no plateau is observed for the secondand higher-order cumulants. Next we discuss the behavior of k n,ab and k ′ n,ab describing p(V b ) Va . The latter is rather spread out as shown in the left panels of Fig. 8 for Pb+Pb and O+O collisions. In principle, p(V b ) Va is not necessarily described by independent source picture, i.e. the value of V b can not be treated as a sum of independent contributions from V a number of sources. In practice, the deviation from the independent source picture is small and the cumulants of p(V b ) Va are approximately proportional to V a . Nevertheless, we need to use Eq. (11) to calculate k ′ n,ab , which includes the firstorder correction to account for the residual dependence of k n,ab on V a . Such corrections can be quite important if the centrality resolution of subevent A is poor (i.e. k n,Va are large).
The right panels of Fig. 8 show k n,ab and k ′ n,ab as a function of V a . The first-order cumulants k 1,ab and k ′ 1,ab are close to unity in mid-central collisions, where ⟨V b ⟩ Va ≈ V a . The rather sharp decrease of k ′ n,ab in central region is due to V a ∂k 1,ab ∂Va in Eq. (11), which are quite large in the central region. Overall the influence of k 1,ab and k ′ 1,ab to the higher-order k n,V b should be small except in the very central and peripheral collisions.
The third column of Fig. 8 shows that the scaled variance k 2,ab increases sharply and then decreases gradually for the remaining V a range. This behavior can be understood from the leaf-like shape of the 2D correlation between V b and V a : although the variance of V b is small towards small and large V a region, the scaled variance of V b (i.e. normalized by ⟨V b ⟩) is large in small V a region and decreases to nearly zero at largest V a . The higher-order cumulants k 3,ab and k 4,ab have more complex behaviors. In the small V a region, k n,ab are positive since V b > 0. At large V a , the fluctuations of V b is bounded by V b ≤ A, leading to a negative k n,ab . Note that the sign-change of cumulants occurs in mid-central region, around V a ≈ 110 for k 3,ab and V a ≈ 60 for k 4,ab for Pb+Pb. This behavior is very different from the influence of pure multiplicity fluctuations in Section V, where the sign-change happens only in the UCC region. For higher-order cumulants, we note that the value of k ′ n,ab are very different from k n,ab over the full range of V a due to V a ∂(k 1,ab k n,ab ) k 1,ab ∂Va in Eq. (11). This implies that the deviation from the independent source assumption has a stronger impact for the third and higher-order cumulants.
In the smaller O+O system, even the k 1,ab and k ′ 1,ab already show a quite sizable deviation from unity. In general, the higher-order cumulants k n,ab and k ′ n,ab have smoother variations but over a much broader V a range. From k n,Va in Fig. 7 and k n,ab and k ′ n,ab in Fig. 8, we calculate the k n,V b via Eq. (10) and compare with the results from full calculation. This comparison is shown in Fig. 9 for Pb+Pb and O+O collisions, together with the breakdown of the contributions from its individual components. The behavior of cumulants and the comparison can be summarized as follows: 2. The results of k 2,V b are shown in the second column. The N a dependence of k 2,V b in Pb+Pb collisions is described nearly perfectly by Eq. (10), while a small deviation is observed in O+O collisions around the maximum of k 2,V b . In Pb+Pb collisions, the decrease of k 2,V b as a function of N a is mainly driven by k 2,ab , while k 2,Va mainly contributes an offset in mid-central region and a sharp decrease in the UCC region. In O+O collisions, these two components have more similar shapes, although the k 2,ab term dominates in the central region and the most peripheral region.
3. The results of k 3,V b are shown in the third column. In Pb+Pb collisions, the description by Eq. (10) is also very good, except that it slightly underestimates the magnitude in the most central collisions. In O+O collisions, Eq. (10) underestimates the magnitude of k 3,V b in the central region, and around its maximum in the peripheral region. Looking at the three individual components of Eq. (10), the genuine FB decorrelation term k 3,ab dominates only in the peripheral region, but its contribution is not as important as the second term associated with k 2,Va Na . The last term associated with k 3,Va Na is nearly a constant, and it has little influence on the shape of k 3,V b , except in central collisions. These behaviors also imply that k 3,ab can only be reliably extracted if k 3,Va Na are small; otherwise, one has to first measure k ′ 2,ab and k ′ 1,ab using the information in the first two columns.
4. The behavior of k 4,V b shown in the right column are much more complex. In Pb+Pb collisions, the description of k 4,V b by Eq. (10) is only good in peripheral and the most central collisions. In O+O collisions, Eq. (10) has poor description over most of the region for the Par1 parameter set. When the Par0 parameter set is used, the agreement is much better in Pb+Pb collisions, but is still relatively poor in O+O collisions. Looking at the contributions of individual terms, the intrinsic kurtosis term k 4,ab from FB fluctuations is important only in the very peripheral and central regions. The mixing terms between lower-order k n,ab and k n,Va dominate in the mid-central collisions and are also important in the other centrality range. Figure 10 summarizes the results in the five collision systems for Par1 parameter set. The results of k n,V b in the bottom row are related to k n,Va in the top row and k n,ab and k ′ n,ab in the middle row via Eq. (10). For the first-order cumulants shown in the left column, the behavior of ⟨V b ⟩ Va largely resembles those of ⟨V a ⟩, except for very small systems where the k 1,ab show a significant deviation from one. For all higher-order cumulants, the values of k n,Va are much smaller than k n,ab . In fact, in the absence of the centrality resolution effect in subevent A, k n,Va vanish and k n,V b approach k n,ab . We find that if Par2 is chosen for multiplicity smearing in subevent A, the k n,Va are much larger and dominate the behavior of k n,V b . The shapes of k n,ab and k ′ n,ab are similar between different collision systems, but the magnitudes are smaller for smaller collision systems. This leads to a similar system-size ordering for the values of k n,V b .
Given the fact that k 1,ab and k ′ 1,ab are close to one, the features of k 2,V b can be used to separate k 2,ab and k 2,Va via Eq. (10): the increase and decrease from peripheral to mid-central collisions mainly reflects the contribution of k 2,ab , while the sharp decrease and flattening out behavior in central region are dominated by the k 2,Va . In contrast, we find that k 3,V b and k 4,V b are much more sensitive to mixing terms between k n,ab , k ′ n,ab and k n,Va , especially in mid-central and central regions, similar to Fig. 9. Separate these different components would be a rather challenging task.
In the final step, the centrality cumulants k n,V b shown in Figs. 9 and 10 are used to obtain the full results of multiplicity cumulants including the smearing effects in subevent B via Eq. (5) or Eq. (13) for Par1. Examples can be found in Sec. B.

VII. SUMMARY AND DISCUSSION
In heavy ion collisions, the centrality or the volume of the fireball, characterized by the number of sources V in the initial state geometry, is not fixed but fluctuates for events with the same final-state particle multiplicity N . This so-called centrality or volume fluctuations (CF) can arise from either 1) fluctuations in the particle production process which smears the mapping between N and V or 2) longitudinal fluctuations of V within the same event which decorrelates the N between different η ranges. In this paper, we propose to study CF using the correlation of multiplicities N b and N a in two subevents separated in pseudorapidity. This two-dimensional correlation are analyzed via cumulants of the conditional probability distribution p(N b ) Na as a function of N a . The contributions from the two types of CFs can be identified from features in the N a dependence of these cumulants.
For a quantitative study of the CF, a standard Glauber model based on an independent source picture is used: the Glauber model is used to produce the V for each event, and the final-state particle multiplicity in each event is then calculated as a sum of particles from each source N = ∑ V i=1 n i with n i sampled independently from a common p(n). For the study on the effects of multiplicity smearings, the sources for subevents A and B are chosen as V a = V b = N part . For the study on the effects of centrality decorrelations, they are chosen as V a = N F part and V b = N B part . The centrality selection is based on N a for both cases, i.e. cumulants are calculated upto fourth-order: mean ⟨V ⟩, scaled variance The ⟨V ⟩ is observed to be linearly proportional to N a , except in ultra-central collisions (UCC) where the increase slows down as N a approaches the upper bound. This means that N a ⟨V ⟩ exhibits a sharp increase in the UCC region similar to what is observed in the experimental data [28,32]. We also verified that (not shown) this linearity is not affected much by the centrality decorrelations considered in this paper, except in smaller systems where the correlation between V a and V b has significant deviations from diagonal. Such behavior in the UCC is very sensitive to the properties of p(V ) and p(n).
When including only effects of multiplicity smearings, the scaled variance k 2,V is nearly constant in mid-central collisions, but decreases to zero in UCC events as N a approaches the upper bound. In the presence of centrality decorrelations, the k 2,V also exhibits a linear decrease in the mid-central collisions. The rate of decrease should be a robust measure of the extent of longitudinal fluctuations in the particle production sources.
The N a dependence of higher-order centrality cumulants are more complex. In the presence of centrality decorrelations, they receive significant contributions from several mixing terms related to the lower-order cumulants, which are also much larger in smaller systems. Furthermore, additional correction terms need to be included to account for the fact that the cumulants for p(V b ) Va is not proportional to V a except for the first order. Due to these reasons, the ⟨V ⟩ and k 2,V are considered most valuable observables to probe the nature of p(V ) and p(n).
From the centrality cumulants k n,V , we calculated the cumulants for final-state multiplicity N b in subevent B, K n,b via Eq. (5). We find that the flattening behavior of ⟨N b ⟩ as a function of N a in UCC remains. Experimentally, one can calculate both ⟨N b ⟩ (N a ) and ⟨N a ⟩ (N b ) from p(N b , N a ) to quantify the relative centrality resolution of the two subevents. If ⟨N b ⟩ (N a ) is more diagonal than ⟨N a ⟩ (N b ), it would imply that N a has a better centrality resolution on V than N b (a similar conclusion was reached from ATLAS data [32]). For the scaled variance K 2,b , we find that most of the features of k 2,V has been preserved, therefore experimentally measured scaled variance in subevent B for events selected in subevent A can also be used to extract information about the centrality decorrelations.
The centrality decorrelations considered in this paper is rather simplistic. This can be improved by considering the partonic degree of freedom in the nucleons and assuming the quark participant number n qp to be a number that vary with η. In addition, one should consider situation used by most collider experiments, where the subevent A for centrality is defined in both forward and backward rapidity, and subevent B is defined at mid-rapidity. In this case, although V a and V b are largely correlated with N part , they should still be decorrelated due to difference in n qp . We leave this to future investigations.
In summary, we have studied the centrality fluctuations in terms of the number of initial sources V using the correlation of the final-state multiplicity between two subevents N a and N b . The sources and final multiplicity are generated with a Glauber model within an independent source picture. This correlation is analyzed in terms of cumulants for multiplicity distribution in one subevent (N b ) for events selected in another subevent with a fixed multiplicity (N a ). The N a dependence of the cumulants are affected by both pure multiplicity smearing in the particle production as well as the decorrelation between the sources in the two subevents. We find that the behavior of cumulants in ultra-central collisions are very sensitive to the particle production for each source p(n) and p(V ) due to the steeply falling p(V ) distribution. The FB centrality decorrelations have very little impact on the mean multiplicity but leads to an decrease of scaled variance ⟨(δV b ) 2 ⟩ ⟨V b ⟩ and ⟨(δN b ) 2 ⟩ ⟨N b ⟩ as a function of N a . These features can be used to constrain the particle production mechanism as well as the longitudinal fluctuations of the initial-state sources. Collect terms in power of t on both sides, one obtains the general formula of Ref. [13], which expresses the cumulants of p(B) in terms of those for p(b) and p(V ), i.e.
where the B n,i are Bell polynomials. This paper considers deviations from the independent source picture, i.e. when χ B V is not exactly proportional to V , where we has Taylor-expanded χ b in ln V aroundV , i.e., Dropping terms that are suppressed by the system size 1 V n , and plugging the first two terms back into Eq. (A7) gives Following the remaining steps of Eq. (A7), we obtain: This result can also be expressed with the normalized cumulant notation, i.e., For large collision systems, this leading-order approximation works very well. For small collision systems, where the variance of CF may not be small in comparison toV , higher-order correction terms in Eq. (A9) should also be considered. The expressions of the corrections for c n,B are quite lengthy, so we only show the result for first and second cumulants up to all order, and the third-order cumulants to 1 V , 1 , with the short-hand notation D (k + 1)!. As a concrete application of Eqs. (A12)-A13, it is useful to consider the CF due to the finite centrality bin-width effect [14]. Assuming centrality is defined in a finite multiplicity range in subevent A, whose distribution is described by p(A), the corresponding CF is given by: In general, p(V ) A is not described by the independent source picture, i.e. the number of V per particle is not independent, even if p(A) V is. But the correlation between V and A is expected to be close to linear, and thus the approximation described by Eq. (A12) or A13 can be directly used.
When one needs to consider both the centrality fluctuation in subevent B and the centrality bin width effect in subevent A, the relations can be expressed as, We consider the probability distribution p(v) of V per particle, whose cumulant generating function χ v (t) is related to those for p(V ) A as Combining Eqs.A8 and A12, we obtain B =Āvb,