Partial correlation analysis method in ultra-relativistic heavy-ion collisions

TL;DR

This paper introduces a partial covariance analysis method for two-particle correlations in ultra-relativistic heavy-ion collisions, effectively removing centrality fluctuation effects and providing clearer insights into initial-state physics.

Contribution

The paper presents a novel application of partial covariance to analyze correlations, allowing centrality constraints on sources and reducing trivial fluctuations in heavy-ion collision data.

Findings

Partial covariance effectively removes centrality fluctuation effects.

The method simplifies analysis of initial-state correlations.

Short-range correlations from resonance decays are largely reduced.

Abstract

We argue that statistical data analysis of two-particle longitudinal correlations in ultra-relativistic heavy-ion collisions may be efficiently carried out with the technique of partial covariance. In this method, the spurious event-by-event fluctuations due to imprecise centrality determination are eliminated via projecting out the component of the covariance influenced by the centrality fluctuations. We bring up the relationship of the partial covariance to the conditional covariance. Importantly, in the superposition approach, where hadrons are produced independently from a collection of sources, the framework allows us to impose centrality constraints on the number of sources rather than hadrons, that way unfolding of the trivial fluctuations from statistical hadronization and focusing better on the initial-state physics. We show, using simulated data from hydrodynamics followed with statistical hadronization, that the technique is practical and very simple to use, giving insight into the correlations generated in the initial stage. We also discuss the issues related to separation of the short- and long-range components of the correlation functions, and show that in our example the short-range component from the resonance decays is largely reduced by considering pions of the same sign. We demonstrate the method explicitly on the cases where centrality is determined with a single central control bin, or with two peripheral control bins.

Equations

$$\mathrm{c}(A,B)$$

$$\mathrm{v}(A)=⟨A^2⟩-⟨A{⟩}^2=\mathrm{c}(A,A).$$

$$\mathrm{c}(X,Y\bullet Z)=\mathrm{c}(X,Y)-\frac{\mathrm{c}(X,Z)\mathrm{c}(Z,Y)}{\mathrm{v}(Z)}.$$

$$\mathrm{v}(A\bullet Z)=\mathrm{v}(A)-\frac{\mathrm{c}(A,Z{)}^2}{\mathrm{v}(Z)}=\mathrm{c}(A,A\bullet Z),A=X,Y.$$

$$C(X,Y)=\frac{\mathrm{c}(X,Y)}{⟨X⟩⟨Y⟩},$$

$$C(X,Y\bullet Z)=C(X,Y)-\frac{C(X,Z)C(Z,Y)}{V(Z)}.$$

$$\rho (X,Y\bullet Z)=\frac{\mathrm{c}(X,Y\bullet Z)}{\sqrt{\mathrm{v}(X\bullet Z)\mathrm{v}(Y\bullet Z)}}$$

$$=\frac{\rho (X,Y)-\rho (X,Z)\rho (Z,Y)}{\sqrt{1-\rho (X,Z{)}^2}\sqrt{1-\rho (Z,Y{)}^2}},$$

$$\rho (A,B)=\mathrm{c}(A,B)/\sqrt{\mathrm{v}(A)\mathrm{v}(B)}.$$

$$\mathrm{c}(X,Y\bullet Z)\simeq \mathrm{c}(X,Y|Z),$$

$${\eta }_{\parallel }=\frac{1}{2}\log \left(\frac{t+z}{t-z}\right)$$

$$\eta \simeq k{\eta }_{\parallel },$$

$$N_A=\sum _{i=1}^{S_A}{m}_i,A=F,B$$

$$⟨S_A⟩⟨m⟩$$

$$\mathrm{c}(S_A,S_{A^{\prime }})⟨m{⟩}^2$$

$$\mathrm{v}(S_A)⟨m{⟩}^2$$

$$C(S_A,S_{A^{\prime }})=C(N_A,N_{A^{\prime }})-{\delta }^{A{A}^{\prime }}\frac{\omega (m)}{⟨N_A⟩}\equiv \overline{C}(N_A,N_{A^{\prime }}).$$

$$C(S_F,S_B\bullet S_C)$$

$$\mathrm{c}(S_A,N_{A^{\prime }})$$

$$C(S_F,S_B\bullet N_C)$$

$$f_{A,B}({\eta }_{\parallel })=\frac{y_b\pm {\eta }_{\parallel }}{y_{\mathrm{b}}}h({\eta }_{\parallel }),\mathrm{for}|{\eta }_{\parallel }|<y_{\mathrm{b}},$$

$$S({\eta }_{\parallel })=\left(Q_{+}+Q_{-}\frac{{\eta }_{\parallel }}{y_b}\right)h({\eta }_{\parallel }),$$

$$C(S_F,S_B)=\frac{\mathrm{v}(Q_{+})}{⟨Q_{+}{⟩}^2}+\frac{\mathrm{v}(Q_{-})}{⟨Q_{+}{⟩}^2}{u}_1{u}_2,$$

$$u_{1,2}=\frac{{\eta }_{\parallel 1,2}}{y_b}=\frac{\eta }{k{y}_b}$$

$$C(S_F,S_B\bullet S_C)=\frac{\mathrm{v}(Q_{-})}{⟨Q_{+}{⟩}^2}{u}_1{u}_2.$$

$$\rho (S_F,S_B\bullet S_C)=\mathrm{sgn}(u_1{u}_2),$$

$$C(S_F,S_B\bullet S_L+S_R)=\frac{\mathrm{v}(Q_{-})}{⟨Q_{+}{⟩}^2}{u}_1{u}_2,$$

$$C(S_F,S_B\bullet S_L,S_R)=0,$$

$$C(S_F,S_B\bullet S_C)={\overline{C}}_s(N_F,N_B)+{\overline{C}}_l(N_F,N_B)-$$

$$\frac{[{\overline{C}}_s(N_F,N_C)+{\overline{C}}_l(N_F,N_C)][{\overline{C}}_s(N_B,N_C)+{\overline{C}}_l(N_B,N_C)]}{{\overline{\mathrm{v}}}_s(N_C)+{\overline{\mathrm{v}}}_l(N_C)}.$$

$${\Sigma }_{XX\bullet Z}={\Sigma }_{XX}-{\Sigma }_{XZ}{\Sigma }_{ZZ}^{-1}{\Sigma }_{ZX},$$

$${\left({\Sigma }_{AB}\right)}_{ij}$$

Full text

relativistic heavy-ion collisions, partial correlation analysis, superposition model

I Introduction

From the outset of the studies of correlations in ultra-relativistic heavy-ion collisions, it has been known that event-by-event fluctuations due to the choice of the centrality class of the sample lead to spurious effects that should be separated from the physical correlations. A simple example involves the forward-backward (FB) fluctuations of the multiplicity of produced hadrons, with the centrality defined from multiplicity in a suitable reference bin, or via some other quantity obtained from detector response, correlated to multiplicity. As the measurements in the physical $F$ and $B$ bins are correlated to the reference bin, the FB correlations determined via a naive definition depend strongly on the fluctuations in the reference bin, thus the character of the “true” FB correlations is obscured.

Several remedies have been proposed to cure the problem of centrality fluctuations. Firstly, if the size of the data sample allows, one may use sufficiently narrow centrality bins, such that centrality fluctuations are negligible. Then, to improve statistics, one may average the obtained covariance matrices over several narrow centrality bins within a broader class. This method, essentially based on the concept of conditional correlation (see, e.g., Cramer (1946)), was successfully used in the analysis by the STAR Collaboration Abelev et al. (2009), with further proposals presented in De et al. (2013).

In this paper we explore the so-called partial correlation analysis (see, e.g., Cramer (1946); Krzanowski (2000)) in application to ultra-relativistic heavy-ion collisions. The method is closely related to the conditional correlations (see App. B), however, it offers an appealing simplicity as well as immediate insight into conditional independence of the studied variables. It has been widely used in other domains of statistical applications, ranging from physics (see, e.g., an interesting example from the X-ray spectroscopy, where partial correlations are used to remove the spurious effects of beam intensity fluctuation Frasinski (2016)) to medicine and psychology. The basic goal of the approach is to assess the correlation (or independence) of some “physical” variables, where the sample is determined by certain fluctuating control (or external or nuisance) variables, whose effects needs to be removed to accomplish the understanding of the relations between the physical variables.

We bring up the relationship of the partial covariance to the conditional covariance Lawrance (1976); Baba et al. (2004), which holds under conditions which are well satisfied in ultra-relativistic heavy-ion collisions. Then, the partial covariance with a control bin may be understood as conditional covariance with the fixed hadron multiplicity in the control bin. An important point, however, is that one should impose centrality constraints at the level of initial state rather than finally produced hadrons; that way, we limit the external fluctuations concerning the initial production, which are of our principal interest. One can accomplish this goal in the framework of the superposition model Olszewski and Broniowski (2013) of particle production, where hadrons are emitted from independent sources. We extend the partial correlation analysis to this physically interesting case. The result is a simple modification of the partial-correlation formulas, where variances have the autocorrelation terms removed. That way one can impose the constraints at the level of the initial sources, which is a non-trivial outcome of the partial covariance method in the superposition approach.

In our study, we use simulated events for ultra-relativistic heavy-ion collisions generated with event-by-event 3+1D viscous hydrodynamics Bożek (2012) and THERMINATOR Kisiel et al. (2006); Chojnacki et al. (2012) to argue that the partial covariance technique is a practical tool to analyze two-particle correlations. We apply it to examples where centrality is determined with the multiplicity in a single mid-rapidity control bin, as well as in two peripheral-rapidity control bins.

We verify in our examples that with the removed effects of resonance decays, the partial FB multiplicity correlations obtained from the simulated data reflect very closely the initial correlations in spatial rapidity, here implemented in the very simple Bzdak-Teaney (BT) Bzdak and Teaney (2013) form for the wounded quark model Bożek et al. (2016). This feature shows that the partial correlation technique may be used in practice to access information on the initial-state correlations. However, as in other approaches, this possibility relies on the separation of the short-range correlations (such as from the resonance decays, jets, or femtoscopic correlations) from the long-range correlations, generated in the initial state.

We recall that another approach separating centrality fluctuations, somewhat different in spirit and derived from the superposition approach, defines the so-called strongly-intensive measures Gaździcki and Mrówczyński (1992); Gorenstein and Gaździcki (2011), which by construction do not depend on fluctuations of the number of sources (or the volume) which emit the observed hadrons. Yet another technique which can be used to accomplish the goal invokes the Principal Components Analysis (PCA) Bhalerao et al. (2015). We discuss the relation of the partial covariance technique to these methods in App. C.

The outline of our paper is as follows: In Sec. II we provide the basic definitions of the partial covariance and partial correlation, and discuss their meaning linked to imposing conditions with control variables. Then in Sec. III we pass to the discussion of the FB multiplicity correlations, followed in Sec. IV by the description of the superposition approach to the multi-stage particle production process. Our key formulas, allowing to impose centrality fixing conditions at the level of sources rather than hadrons, are derived there. In Sec. V we obtain partial correlations in the initial state for the BT model, which are confronted with the corresponding results obtained with simulated data in Sec. VI. In Sec. VII we recapitulate our results, summarizing in practical terms the method based on partial correlations, which can be used in experimental analyses. The appendices contain the discussion of the relation of the partial covariance analysis to other methods, as well as some technical details.

II Partial correlation

In this section we establish the notation and provide the standard definitions, with more details presented in Appendix A.

The simplest case of partial correlations Cramer (1946); Krzanowski (2000) involves two physical random variables, $X$ and $Y$, and a single control random variable $Z$. One defines the element of the covariance matrix in the standard way as

[TABLE]

where $\langle.\rangle$ denotes the averaging over the sample of $n$ events. The variance is, of course, the diagonal term,

[TABLE]

The partial covariance between $X$ and $Y$ with the control variable $Z$ is defined by the formula (see Appendix A for interpretation)

[TABLE]

The second term in Eq. (2) removes the piece of the covariance between $X$ and $Y$ which is due to their correlation to $Z$, by means of projecting out the components of $X$ and $Y$ which are parallel to $Z$ (in the $n$-dimensional space, with $n$ denoting the number of events). The partial variance is, correspondingly,

[TABLE]

In experimental studies one often uses the correlation function defined as the covariance scaled with the multiplicities, i.e.,

[TABLE]

and, correspondingly, $V(A)={\rm v}(A)/\langle A\rangle^{2}$. Then the partial $C$-correlation following from Eq. (2) is

[TABLE]

Finally, one defines the partial analog of Pearson’s $\rho$-correlation Cramer (1946); Krzanowski (2000),

[TABLE]

where

[TABLE]

As discussed in a greater detail in Appendix B, the partial covariance, under quite general assumptions Lawrance (1976); Baba et al. (2004) which are typically fulfilled in ultra-relativistic heavy-ion collisions, is related to the conditional covariance, namely

[TABLE]

where the meaning of the condition in ${\rm c}(X,Y|Z)$ is that first $Z$ is fixed to a very narrow class (for instance, if it describes the discrete multiplicity of hadrons in the reference bin, it can be fixed to a natural number equal to this multiplicity), then the covariance between $X$ and $Y$ is evaluated within this subsample, and finally averaging of thus obtained covariances over various values of $Z$ within the sample is performed.

As a matter of fact, such a conditional procedure was used in the STAR experiment Abelev et al. (2009) to analyze the FB multiplicity correlations in Au+Au and $p$+$p$ collisions at $\sqrt{s_{NN}}=200$ GeV. In this context, relation (8) was derived by Lappi-McLerran Lappi and McLerran (2010) under the assumption of normal distributions, and by Bzdak Bzdak (2012) with the condition (38), similarly as in Lawrance (1976); Baba et al. (2004).

The practical significance of Eq. (2) or (5) is that the imposition of external constraints (such as fixing centrality) may be, under general assumptions, accomplished via the partial covariance technique.

III Forward-backward multiplicity correlations

In our study, $X$, $Y$, and $Z$ random variables are the multiplicities of produced hadrons in, correspondingly, a forward pseudorapidity bin $F$, a backward pseudorapidity bin $B$, and a reference bin. We will explore two cases: reference bin $C$ located at a central pseudorapidity bin, and the sum of bins $L$ and $R$, located symmetrically at peripheral rapidities $\eta_{L}$ and $\eta_{R}=-\eta_{L}$. Multiplicity in the reference bin determines the centrality of the event. The location of the reference bin is fixed, whereas $F$, centered around $\eta_{1}$, and $B$, centered around $\eta_{2}$, vary in the covered pseudorapidity range (in this more general arrangement we do not request that $F$ is forward and $B$ backward, but they assume any location in pseudorapidity within the acceptance range). We will also consider the case where the multiplicities in $L$ and $R$ bins are fixed independently, according to the case with two constraints described in App. A.

Experimental studies of the longitudinal multiplicity correlations have a long history. Early investigations of $pp$ and $p\bar{p}$ collisions Uhlig et al. (1978); Alpgard et al. (1983); Alner et al. (1987); Aivazian et al. (1989); Ansorge et al. (1988); Derado et al. (1988); Alexopoulos et al. (1995) and nuclear collisions Bächler et al. (1992); Akiba et al. (1997) were followed by the studies of ultra-relativistic heavy-ion and $pp$ reactions at RHIC Back et al. (2006); Abelev et al. (2009) and at the LHC Aad et al. (2012); Jia et al. (2016); Adam et al. (2015); Jia (2016); Aaboud et al. (2016); Sputowska (2016); Adam et al. (2016). Numerous physical models and theoretical methods have been invented in attempts to understand the mechanisms behind the generation of long-range correlations Capella and Krzywicki (1978); Kaidalov and Ter-Martirosian (1982); Chou and Yang (1984); Capella et al. (1994); Amelin et al. (1994); Braun et al. (2000); Giovannini and Ugoccioni (2002); Braun et al. (2004); Haussler et al. (2007); Brogueira and Dias de Deus (2007); Armesto et al. (2007); Vechernin and Kolevatov (2007); Braun (2008); Konchakovski et al. (2009); Bzdak and Woźniak (2010); Lappi and McLerran (2010); Bzdak (2009); Bożek et al. (2011); Yan et al. (2010); Dias de Deus and Pajares (2011); Białas and Zalewski (2011); Bzdak (2012); Bzdak and Teaney (2013); Vechernin (2015a); De et al. (2013); Ma and Bzdak (2014); Bzdak and Bożek (2016); Vechernin (2015b); Monnai and Schenke (2016); Broniowski and Bożek (2015); He et al. (2016a, b); Olszewski and Broniowski (2017). The importance of these investigations lies in the well-known fact that the correlations over a large rapidity separation can originate only from the earliest stages of the collision, thus they may reveal fingerprints of the early dynamics of the system.

In our study we demonstrate the methodology based on partial correlations with simulated events for Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, described in detail in Sec. VI.

IV Superposition approach

One should bear in mind that the particle production in ultra-relativistic heavy-ion collisions is an effect of a multi-stage process (see, e.g., Olszewski and Broniowski (2013, 2015)). First, we have an early production of entropy, resulting from partonic physics. It is modeled in various approaches, such as string formation Amelin et al. (1994); Brogueira and Dias de Deus (2007), the Color Glass Condensate theory Kovner et al. (1995); Iancu et al. (2001), or the Glauber approach Glauber (1959); Czyż and Maximon (1969); Białas et al. (1976); Kharzeev and Nardi (2001); Back et al. (2002); Broniowski et al. (2007), used in this work. The initial partonic entropy is distributed in space in a correlated way.

This distribution, treated event-by-event, serves as an initial condition for the intermediate hydrodynamic evolution (for reviews see, e.g., Heinz and Snellings (2013); Gale et al. (2013) and references therein) or for transport modeling Lin et al. (2005). In our work the applied hydrodynamics is deterministic, hence it does not introduce extra fluctuations, which should arise in a viscous system Kapusta et al. (2012), but which according to estimates are not very significant Yan and Gronqvist (2016); Gavin et al. (2016).

The intermediate evolution continues until freeze-out, where the Cooper-Frye formalism Cooper and Frye (1974) is applied at the freeze-out hyper-surface defined with a constant temperature (we take $T_{f}=150$ MeV). It generates primordial hadrons (stable and resonances) according to a thermal distribution. Subsequently, resonances undergo decays, possibly in cascades, into stable particles. Due to a statistical nature of the production process, the distribution of a hadron of a given species is Poissonian, hence (trivial) fluctuations are generated due to the sampling with a finite number of particles.

The goal of the data analyses is to unfold of the trivial fluctuations Gaździcki and Mrówczyński (1992); Fu and Liu (2003); Białas (2004); Gorenstein and Gaździcki (2011); Ling and Stephanov (2016); Nonaka et al. (2017), such as those from statistical hadronization, and acquire information on correlations generated in the earlier evolution phases. Below we describe how this is accomplished in the superposition approach.

Let us first bring up an important approximation underlying this approach, which may be termed as no bin mixing. The initial distribution of entropy may be divided into cells labeled with their space-time rapidity

[TABLE]

(and also by the transverse coordinates $x$ and $y$). The bins are “carried over” with hydrodynamics or transport to the freeze-out hyper-surface, with the assumption of no mixing between the bins, such that the final pseudorapidity of the fluid element, $\eta$, is a function of $\eta_{\parallel}$. The hydrodynamic push in the longitudinal direction is not very strong, and we estimate

[TABLE]

with $k=1.2$ for the model we apply.

At freeze-out, there is some thermal dispersion of the momenta of hadrons, as particles originating from the same hydrodynamic fluid cell acquire rapidities spread with $\Delta\eta\sim 1$. This causes some bin mixing of the two-particle correlation function in pseudorapidity. In addition, resonance decays generate extra short-range correlations with the width of $\Delta\eta\sim 1$ (we will return to the issue of separating the short-range components when describing our results in Sec. VI).

We remark that inclusion of the detector acceptance into the framework does not lead to a new element. If the acceptance is expressed with a Bernoulli trial of success rate $p$, then its folding with a Poisson distribution from the thermal motion leads to a Poisson distribution with a mean enhanced by the factor $p$. Sums of various particle species with Poisson distributions also lead to a Poisson distribution with the mean expressed as a sum means of the added distributions.

Our basic methodology is as follows: we will use simulated data (with removed short-range component coming from resonance decays) to obtain the two-particle partial correlation function with the trivial fluctuations unfolded with the help our formalism. The result will be directly compared to the partial correlation function of the initial state used in the simulations (cf. Sec. V), with arguments shifted according to Eq. 10.

We now come to the derivation of relevant formulas. One may consider the fluid cells at the freeze-out hyper-surface as sources, which emit hadrons. Each source, by definition, emits particles independently of other sources, but with the same (for the sake of simplicity) distribution. The number of fluid cells (sources) at the forward $F$ and backward $B$ pseudorapidity is denoted as $S_{F}$ and $S_{B}$, respectively (recall that the fluid cells at a given rapidity are located at various transverse positions and their number fluctuates). Then, the number of produced hadrons in each bin is

[TABLE]

where $m_{i}$ denotes the number of particles produced from source $i$. We assume that each sources produces particles with the same distribution, hence $\langle m_{i}\rangle=\langle m\rangle$ and ${\rm v}(m_{i})={\rm v}(m)$. The assumption of the independent production leads to simple formulas which connect moments of the number of sources with the moments of the number of the produced particles Olszewski and Broniowski (2013, 2015):

[TABLE]

where $\omega(X)={\rm v}(X)/\langle X\rangle^{2}$ is the scaled variance. We notice a subtraction of $\omega(m)\langle N_{A}\rangle$ appearing for the diagonal term $A=A^{\prime}$, i.e., for the variance:

[TABLE]

The origin of this term is the presence of variance (or autocorrelation) of the particles produced from the same sources, ${\rm v}(m_{i})={\rm c}(m_{i},m_{i})>0$, whereas the covariance of the overlaid distributions from different sources vanishes by the assumption of independent production, ${\rm c}(m_{i},m_{j})=0$ for $i\neq j$.

For the special case of the Poisson distribution of the overlaid variable $m$ (which is the case of our numerical study presented in Sec. VI)), we have $\omega(m)=1$ and $\overline{{\rm v}}(N_{A})=\langle N_{A}(N_{A}-1)\rangle-\langle N_{A}\rangle^{2}$, which corresponds to the subtraction of autocorrelations; the average number of pairs appears in the formula. However, Eq. (13) is more general, holding, under the assumptions of the superposition model, for any distribution of the overlaid variable $m$. In this paper, we refer to the subtraction of $\omega(m)\langle N_{A}\rangle$ from the variance as to “the removal of autocorrelations”, and indicate it with a bar.

In an experimental study, where in principle one does not know the distribution of $m$, one may find $\omega(m)$ in a numerical way with the following procedure: One examines the covariance matrix ${\rm c}(N_{A},N_{A^{\prime}})$ as a two-dimensional matrix in $A$ and $A^{\prime}$ indices. The diagonal term at $A=A^{\prime}$ forms a sharp discontinuous ridge, sticking up from a smooth ”background” function. One then adjusts the value of $\omega(m)$ to remove the sharp ridge from the function ${\rm c}(N_{A},N_{A^{\prime}})-\delta^{AA^{\prime}}\omega(m)\langle N_{A}\rangle$. This prescription conforms to Eq. (12), as the covariance of the sources ${\rm c}(S_{A},S_{A^{\prime}})$ is a smooth function of $A$ and $A^{\prime}$.

The meaning of Eq. (13) is that the covariance of the number sources is proportional to the covariance of the observed hadron multiplicities, but with the autocorrelations removed. Passing to the scaled covariance (4), we have

[TABLE]

We are now ready to build the partial covariance for the superposition approach. Some introductory discussion is in place. As stated in Sec. II and further discussed in Appendix B, the meaning of the partial covariance is, essentially, an imposition of a condition on the value of the control variables. In the case of multiplicity correlations in ultra-relativistic heavy-ion collisions a first instinct is to constrain the centrality fluctuations, i.e., the fluctuations of number of hadrons in a reference bin. However, a more desired constraint concerns the number of sources corresponding to the reference bin. Such a constraint is more directly related to a physical situation in the initial state. Suppose a reference bin, defined as $C$, has in a given event $S_{C}$ sources which “determine the physics”. On the other hand, the number of detected hadrons $N_{C}$ is sensitive to the fluctuation in the production from sources. Therefore a fixed value of $N_{C}$ corresponds to event-by event fluctuating values of $S_{C}$; constraining $N_{C}$ does not completely constrain $S_{C}$ and the physics of the initial condition remains washed out.

According to our formalism, the constraint imposed at the level of initial sources is realized with the equation

[TABLE]

This is the key formula used in our analysis of the simulated data in the following sections.

If one wishes, however, to impose the constraint at the level of the produced hadrons, then using the formula

[TABLE]

one arrives at the expression

[TABLE]

Note that the subtle but important difference between Eq. (15) and (16) is in the denominator of the subtracted term, where we find the variance of the multiplicity in the reference bin with autocorrelations subtracted, or present. In our sample, the autocorrelations increase the variance by $\sim 100$%, hence the effect is very important.

We stress that despite its simplicity, the meaning of Eq. (15) is non-trivial, as it allows to impose a strict centrality constraint at the level of sources and infer partial correlation of sources, whereas the evaluation is based solely on measured multiplicities of the produced hadrons.

V Modeling initial correlations

For our illustrative purposes we use the wounded quark model Białas et al. (1977); Anisovich et al. (1978) for the initial state. In this model, the initial sources are the wounded quarks, moving forward or backward, according to the motion of their parent nucleons from nuclei $A$ and $B$. An advantage of using the wounded quarks compared to wounded nucleons Białas et al. (1976) is that one obtains proper scaling Eremin and Voloshin (2003); Kumar Netrakanti and Mohanty (2004); Adler et al. (2014); Adare et al. (2016); Lacey et al. (2016); Bożek et al. (2016) of the multiplicities on the number of participants with no need for the binary-collision component Kharzeev and Nardi (2001). The event-by-event distribution of wounded quarks in the transverse plane is obtained from the Glauber simulations with GLISSANDO Broniowski et al. (2009); Rybczyński et al. (2014), corresponding the transverse location of the wounded quarks. The longitudinal profile in spatial rapidity $\eta_{\parallel}$ is taken according to the model of “triangles” Białas and Czyż (2005); Adil et al. (2006); Bożek and Wyskiel (2010), where each source has the entropy distributed preferentially in the direction of its motion, according to a simple formula

[TABLE]

where $A$ and $B$ denote the sources belonging to, respectively, the left- and right- moving nuclei, and $y_{\rm b}$ is the rapidity of the beam ($y_{\rm b}\simeq 8$ for Pb+Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV), and $h(\eta_{\parallel})$ is an additional profile, typically of a flatten Gaussian form Bożek and Wyskiel (2010). As $h(\eta_{\parallel})$ cancels from the formulas for symmetric $A$+$B$ collisions, we do not need to specify it explicitly.

Let us introduce the notation $Q_{A}$ and $Q_{B}$ for the number of wounded quarks belonging to the $A$ and $B$ nuclei, and $Q_{\pm}=Q_{A}\pm Q_{B}$. Then, according to Eq. (17), the number of sources (combined from nucleus $A$ and $B$) at location $\eta_{\parallel}$ is

[TABLE]

For symmetric collisions the average over events yields $\langle Q_{-}\rangle=0$, hence $\langle S(\eta_{\parallel})\rangle=\langle Q_{+}\rangle h(\eta_{\parallel})$.

Bzdak and Teaney (BT) Bzdak and Teaney (2013) computed the correlation function in the model of triangles (in the variant with wounded nucleons). It yields a very simple result

[TABLE]

where

[TABLE]

(cf. Eq. (10)), and indices $1$ and $2$ correspond to labels $F$ and $B$, respectively. Note that, as announced, the overall rapidity profile $h(\eta_{\parallel})$ cancels out. The moments of $Q_{\pm}$ are read out from GLISSANDO simulations via averaging over some chosen class of events.

Next, we use Eq. (5) to derive the partial covariance function for the BT model with the control bin $C$ placed at $\eta_{\parallel}=0$. A short calculation yields

[TABLE]

We note that Eq. (21) differs from Eq. 19 by not carrying the term with ${\rm v}(Q_{+})$. This is clear from the point of view of the conditional correlation, as in the present case $C(S_{F},S_{B}|S_{C})$ corresponds to fixing the multiplicity at $\eta_{\parallel}=0$. From Eq. (18) we get $S(0)=Q_{+}h(0)$, hence this is equivalent to fixing $Q_{+}$, thus ${\rm v}(Q_{+})=0$ in the calculation of the conditional covariance. The obtained consistency verifies in an obvious way the relation $C(S_{F},S_{B}\bullet S_{C})=C(S_{F},S_{B}|S_{C})$.

For the $\rho$-correlation of Eq. (6) we find the very simple formula

[TABLE]

where ${\rm sgn}$ denotes the sign function. This means, that the partial $\rho$-correlation of multiplicities in bins located at rapidities of the same (opposite) sign is $+1$ ($-1$), indicating maximum correlation (anti-correlation).

In experiments, it frequently happens that multiplicities in peripheral forward ($R$) and distant backward ($L$) bins are available and can be used for centrality determination. Below we consider two cases: 1) where the sum of the multiplicities in $L$ and $R$ is taken and Eq. (5) is used, and 2) the multiplicities in $L$ and $R$ are taken as separate constraints according to Eq. (36). We consider the case where the peripheral bins are symmetrically arranged, with $R$ and $L$ around pseudorapidities $\eta_{R}$ and $\eta_{L}=-\eta_{R}$.

In case 1) the BT model yields exactly the same result as Eq. (21),

[TABLE]

which follows from the fact that according to Eq. (18) $S(\eta_{L})+S(\eta_{R})\sim{Q_{+}}$, and the condition with symmetrically arranged peripheral bins fixes $Q_{+}$, as in the case of the central bin.

In case 2) Eq. (36) gives for the BT model a vanishing result,

[TABLE]

which is compatible with the simultaneous constraints $S(\eta_{L})=0$ and $S(\eta_{R})=0$, which fixes both $Q_{+}$ and $Q_{-}$, hence ${\rm v}(Q+)={\rm v}(Q_{-})=0$, and $C(S_{F},S_{B}|S_{L},S_{R})=0$. We note that this is a specific feature of the BT model, whereas models which include extra fluctuations in the initial state, e.g., the fluctuating-string model of Broniowski and Bożek (2015), would yield a non-zero result.

VI Results for the simulated events

Our final goal is to test to what extent the formulas obtained in the previous section are reproduced with simulated data and the master formula (15) and its equivalent for the case of left and right peripheral bins. Because of departures from assumptions of the superposition approach (bin mixing, resonance decays) this is not an academic exercise, but a verification if the method may be practical in actual data analyses.

The simulated data are obtained as follows: for the initial condition, we apply the wounded quark model for Pb+Pb collisions at $\sqrt{s_{NN}}=2.76\textrm{ TeV}$, as described in Sec. V. We take a rather broad sample corresponding to centrality (as determined with $Q_{+}$) of 30-40%. Then we use the results of event-by-event 3+1D viscous hydrodynamics Bożek (2010) obtained with the wounded-quark initial conditions. The statistical hadronization at freeze-out is carried out with THERMINATOR Kisiel et al. (2006); Chojnacki et al. (2012), which incorporates all hadrons from the Particle Data Tables and implements resonance decays. We label the results obtained with the products of resonance decays as “all charged”, which consists of $\pi^{\pm}$, $K^{\pm}$, $p$, and $\bar{p}$. To reduce the correlations induced by the resonance decays, we also present the results obtained with the “primordial” particles generated at freeze-out, i.e., before the resonance decays.

As we wish to have a view as broad as possible, we take a large continuous acceptance window in our numerical study, with $|\eta_{1,2}|\leq 5.1$ divided into 51 physical bins of width of 0.2. Our coverage is thus (on purpose) much larger than accessible in the LHC experiments, to better illustrate the approach. The central reference bin, where the number of hadrons in $N_{C}$, is taken as $\eta_{C}\leq 0.5$, thus it has a width of five physical bins (other values of width of the reference bin could be taken, with similar results up to statistical uncertainties).

Figure 1 presents the basic output from the simulated data, namely the $\overline{C}$-correlation (with autocorrelations removed), obtained with (a) all charged particles (i.e., after resonance decays), and (b) for the primordial particles. In case (a) we note a hallmark ridge, wide by about one unit of pseudorapidity, extending along the diagonal. It is due to the resonance decays, which is the basic difference between cases (a) and (b). We also note other non-trivial structures in the $\overline{C}$-correlation arising from the applied hydrodynamic model, such as its rise at the boundaries, however understanding these intricate details is not the goal of this work. Rather, we take the simulated data (which, as we see, are not trivial) as they are, and then carry out the partial correlation analysis outlined in the previous sections.

In Fig. 2 we compare the partial $C$-correlation function for the source multiplicities with the central reference bin, obtained for the BT model with Eq. (21), shown in panel (a), to the same quantity obtained from the simulated data with primordial hadrons with Eq. (15). We note that the overall agreement of (a) and (b) is quite remarkable. Hence, if we had the possibility of separating the resonance decays (and in real data, also other sources of short-range correlations from later stages of the evolution) we could infer the initial state correlations with the presented methodology. Note that the application of Eq. (15) leads to large cancellations when passing from Fig. 1(b) to Fig. 2(b). Also, the non-uniformities are “miraculously” smoothed out, hence they had to originate from centrality fluctuations.

Separating other sources of correlations is of course far from simple. One method (as done by the ATLAS collaboration in Aad et al. (2015a)), is to do a numerical fit to a function describing the short-range corrections, and then simply subtract it. In our illustration, it would correspond to “skimming” the ridge in Fig. 1(a), and then carrying out the calculation, which would lead to a result similar to Fig. 2(b).

We note that our model does not account for the correlations reflecting the conservation laws, which is not essential for testing the partial covariance approach.

When the short-range effects from the resonance decays are kept, i.e., we are using the data from Fig. 1(a), then the resulting partial $C$-correlation with the central reference bin has the shape shown in Fig. 3(a). We note the ridge along the $\eta_{1}=\eta_{2}$ diagonal, but also several other features. First, we can see a depletion along the lines $\eta_{1}=0$ and $\eta_{2}=0$. This is a simple artifact of the central reference bin placed at $\eta=0$, since the definition (5) has the feature, that when one of the measurement bins is equal to the reference bin, e.g., $Y=Z$, then $C(X,Z\bullet Z)=0$ identically. Suppose we decompose the measured correlation function into the short- ($s$) and long-range ($l$) components, $\overline{C}=\overline{C}_{s}+\overline{C}_{l}$. Then Eq. (15) becomes

[TABLE]

We see from here all the features coming out in Fig. 3(a). When $F\simeq B$, we get the ridge from $\overline{C}_{s}(N_{F},N_{B})$, and when $F\simeq C$ or $B\simeq C$, the function drops to zero for the reason discussed above Eq. (25). Note, however, that even when $F$, $B$, and $C$ are all sufficiently separated, we get still an artifact from the presence of ${\rm v}_{s}$ in the denominator of the second term in Eq. (25). It leads to a reduction of the subtraction, and, consequently, larger values of $C(S_{F},S_{B}\bullet S_{C})$ in Fig. 3(a) (in the region where $F$, $B$, and $C$ are separated) than in the BT model shown in Fig. 2(a).

Similar distortions are seen in the case of the $L+R$ reference bin, displayed in Fig. 3(b). We note that $C(S_{F},S_{B}\bullet S_{L}+S_{R})$ is “pulled down” at the peripheries, and too large in the region where $F$, $B$, $L$ and $R$ are well separated.

The conclusion of the above discussion is that the short-range components must be separated at the level of $\overline{C}(N_{F},N_{B})$ for the presented analysis to make practical sense.

A well-known method of reducing correlations from the resonance decays is to use particles of the same charge. We apply our procedures for $\pi^{+}$, as there are no resonances that decay into $\pi^{+}\pi^{+}$ pairs. Only some remnant correlation is expected from resonance decays proceeding in cascades. The result is shown in Fig. 4, with the central reference bin in panel (a), and sum of the peripheral bins $L+R$ in panel (b). We note a very weak correlation from cascade decays, visible along the $y_{1}=y_{2}$ diagonal. In the region away from the diagonal we see a very close agreement with Fig. 2. Therefore the use of hadrons of the same sign is an efficient way of getting rid of the short-range correlations from the resonance decays.

In Fig. 5 we show the partial $\rho$-correlation of Eq. 6, with the central reference bin. We recall that in the BT model it is given by Eq. (22), visualized in panel (a). As expected, this shape is reproduced within statistical noise by the primordial particles of the simulated data. We remark that the $\rho$-correlation, being the ratio of the covariance and variance, carries less information than the $C$-correlation. From the $C$-correlation we can read off both the scaled covariance and variance (along the diagonal).

Finally, we test the partial correlation approach for the case of the two peripheral reference bins discussed in App. A, with $-6.1<\eta_{L}<-5.1$ and $5.1<\eta_{R}<6.1$. We recall that in the BT model the corresponding partial correlation vanishes, cf. Eq. (24). This is also the case, within numerical uncertainties, for the simulated data, as visualized in Fig. 6, for the case of primordial particles (a) and all positively charged pions (b).

VII Conclusions

We have presented a simple method capable of providing information on the initial two-particle multiplicity correlations, which is insensitive to centrality fluctuations. The basic formalism relies on the concept of partial covariance, where a reference bin (or a few reference bins) are used to impose constraints on the sample. Application of the method to a superposition model, where particle production occurs in subsequent stages, allows to unfold the trivial statistical fluctuations from the hadronization at freeze-out and gain insight into the correlations in the initial stage of the reaction.

We have demonstrated the feasibility of the method by carrying out an illustrative analysis on simulated data obtained with hydrodynamics, run event-by-event on initial conditions provided by the wounded quark model, and followed by statistical hadronization. We have shown that performing the calculations for the partial correlations of hadrons which do not carry correlations from resonance decays (primordial hadrons, or pions of the same charge) reproduces efficiently the initial partial correlations. Thus the method can be used as a practical tool in experimental data analysis of two-particle correlations.

A nontrivial aspect of our approach is a simple way of unfolding trivial statistical fluctuations generated at statistical hadronization. In the superposition model, it amounts to removal of autocorrelations from the building blocks of the partial correlation function. Moreover, that way we are able to impose constraints on the number of sources in the reference bin, rather than the number of hadrons, which is desirable from the point of view of studying the initial state.

The method is directly extendable to imposition of more constraints, related to a possible involvement of more detectors. This allows for getting more information of the correlations generated in the initial state.

A crucial element of the correlation analyses is the separation of the short-range component, expected to be generated in later stages of the collision, and the long-range component, generated in the initial phase. We have demonstrated on the simulated data that the use of same-charge pions largely reduces the correlations due to resonances. Other sources of short-range correlations (jets, femtoscopy) should be removed with a suitable method.

To summarize, the procedure of obtaining the partial correlations in the initial state is as follows:

Obtain the two-particle correlation function in (pseudo)rapidity from the data. 2. 2.

Remove autocorrelations and the short-distance component. 3. 3.

Apply Eq. (15), or its generalizations in the case of more control bins. 4. 4.

Up to corrections from bin mixing and spatial rapidity – pseudorapidity mapping, the result represents the correlations in the initial state of the collision with centrality fluctuations removed at the level of sources.

Finally, we note that the technique of partial correlations is applicable to other observables which correlate to the “centrality” determination, for instance various charges or the transverse momentum.

Acknowledgements.

We are grateful to Piotr Bożek for providing us a sample of the hydrodynamic simulation results for the wounded quark model, used for illustration of the data analysis methods discussed in this paper. This research was supported by the Polish National Science Centre grant 2015/19/B/ST2/00937.

Appendix A Partial covariance with more control variables

In this appendix we list some basic definitions and properties referring to the partial correlations. In a general case we have $n$ physical random variables ${X}=\left(X_{1},\ldots,X_{n}\right)$ and $m$ control random variables ${Z}=\left(Z_{1},\ldots,Z_{m}\right)$. The quantities $X_{i}$ and $Z_{j}$ are vectors in the $N_{\rm ev}$-dimensional space, where $N_{\rm ev}$ (the number of events) is the number of the data points. Averaging over events for a quantity $A$ is defined as $\langle A\rangle=1/N_{\rm ev}\sum_{k=1}^{N_{\rm ev}}A_{i}$. One defines the partial covariance matrix as

[TABLE]

where $\Sigma_{{AB}}$ is the standard covariance matrix defined as

[TABLE]

where $i$ and $j$ label the variable types. Mathematically, Eq. (26) corresponds to projecting out from the vectors ${X}$ the components belonging to the space spanned by the vectors ${Z}$ (shifted to their central values). Indeed, introducing the projected physical vectors

[TABLE]

which by construction gives the orthogonality condition

[TABLE]

we straightforwardly verify Eq. (26).

The above formulas simplify when the control vectors $Z_{j}-\langle Z_{j}\rangle$ are orthonormal (which we can always achieve via the Gramm-Schmidt procedure or by diagonalization). Denoting these orthonormal vectors as $U_{j}\equiv Z_{j}-\langle Z_{j}\rangle$, $j=1,\dots,m$, we have

[TABLE]

and

[TABLE]

From construction, the diagonal terms are non-negative, $\left(\Sigma_{{XX}\bullet{Z}}\right)_{ii}\geq 0$, where the equality occurs when $X_{i}-\langle X_{i}\rangle$ is contained in the space spanned by the vectors $U_{k}$.

The interpretation of Eq. (31) is straightforward: we subtract from the covariance of the physical variables the covariance proceeding via a correlation, one-by-one, to the control variables $U_{j}$. The meaning of the general case (26) is the same, with the complication arising from the non-orthonormality of the control variables. Some further mathematical facts, in particular the connection to the linear regression analysis, may be found in wiki:partial.

Throughout the paper we use the short-hand notation

[TABLE]

The partial covariance matrix scaled with the multiplicities is defined as

[TABLE]

The Pearson-like partial correlation coefficient between $X_{i}$ and $X_{j}$ is defined as

[TABLE]

which makes sense as long as $c(X_{i},X_{j}\bullet{Z})>0$. From the Schwartz inequality $-1\leq\rho(X_{i},X_{j}\bullet{Z})\leq 1$, and for the diagonal terms $\rho(X_{i},X_{i}\bullet{Z})=1$, as in the case of the standard Pearson’s correlation $\rho(X_{i},X_{i})$.

For the simplest case of a single control variable, the above formulas reduce to Eqs. (2) and (6). For the special case of two physical and two control variables, explored in this paper, we have

[TABLE]

For the case of the scaled covariance Eq. (35) becomes

[TABLE]

where $V(Z_{j})={\rm v}(Z_{j})/\langle Z_{j}\rangle^{2}=C(Z_{j},Z_{j})/\langle Z_{j}\rangle^{2}$.

For the partial correlation coefficient the explicit formula reads

[TABLE]

Appendix B Relation of the partial covariance to the conditional covariance

Lawrance Lawrance (1976) has shown that if a sample satisfies the affine condition

[TABLE]

with $\alpha$ a constant and ${B}$ a constant matrix, then the equality of the partial covariance and the conditional covariance follows,

[TABLE]

The converse was shown by Baba, Shibata, and Sibuya Baba et al. (2004), hence conditions (38) and (39) are equivalent. The question if (38) holds, can be tested on the actual data. In the context of ultra-relativistic heavy-ion collisions, it should hold for sufficiently narrow centrality classes Bzdak (2012). Also, condition (38) holds for normal distributions Baba et al. (2004).

The important practical implication of Eq. (39) is that the partial covariance method is a simple way of imposing constraints in the correlation analysis.

Appendix C Relation to other methods

One may relate the partial covariance method to the Principal Components Analysis (PCA) Bhalerao et al. (2015). Suppose all bins are the measurement bins $X$, obtained, for instance, by dividing the full pseudorapidity acceptance of the detector info narrower bins. Let $U_{j}$ denote the eigenmodes of the covariance matrix $\Sigma_{XX}$. Then one may project out the eigenmodes (with the highest eigenvalues) from the covariance matrix, according to Eq. (31). Thus, PCA may be viewed as a special case of the Partial Covariance method, where the constraints have the form of the eigenmodes of the covariance matrix.

Whereas the algebra (projection) in the two methods is the same, the accents are different. In PCA one does not separate the measurement and reference variables, treating all the bins “democratically” and focuses on a possible hierarchy in the eigenvalues (large gaps in the spectrum are linked to strong correlations, or collectivity, in the fluctuations). In the partial covariance method the measured variables are grouped from the outset into the physical bins $X$ and the reference bins $Y$. This is more natural when the reference data come from different detectors (e.g., in ultra-relativistic heavy-ion collisions, from the peripheral detectors, transverse energy colorimeters, etc.) than the detector collecting the physical data (central TPC). In this case the constraints come entirely from the reference variables, and the constraint vectors $U_{j}$ are built from the vectors $Y$ (see App. A).

In PCA applied to pseudorapidity fluctuations, the multiplicity eigenmode with a highest eigenvalue, denoted in Ref. Bhalerao et al. (2015)as $v_{0}^{(1)}$, corresponds to an $\eta$-independent fluctuation, i.e., where multiplicities in all the bins vary by an equal amount. This may be interpreted as centrality fluctuation, hence the removal of the highest eigenvalue mode in PCA is equivalent to getting rid of the centrality fluctuations, with centrality defined via the total multiplicity from all the bins.

The issue of separating centrality fluctuations is also focal the methods aiming at strongly-intensive fluctuation measures Gaździcki and Mrówczyński (1992); Gorenstein and Gaździcki (2011); Sangaline:2015bma; Broniowski:2017tjq. The framework applied there is also based on the superposition model, but there is one common source (with its multiplicity traditionally termed as the “volume”) and two types of particles, $A$ and $B$, emitted from the source with multiplicities $m_{A}$ and $m_{B}$, respectively. One may then use the superposition approach to relate the statistical moments of $m_{A}$ and $m_{B}$ to the corresponding moments of the observed multiplicities of $A$ and $B$, $N_{A}$ and $N_{B}$, in such a way the fluctuation of the sources cancels out. Thus the object of the study is the emission process from the source, and not the fluctuation of the number of sources (volume), considered trivial.

With the method described in this paper we study the case where there are multiple types of sources (each in a given rapidity bin) and (for simplicity) just one type of hadrons (an extension to more types of hadrons is possible). Our objective is the fluctuation of the number of sources (at various rapidities), whereas the fluctuations of the overlaid variable $m$ are considered not interesting. Thus, the partial covariance method used with the superposition approach in this work is, in a sense, complementary to the studies based on the strongly-intensive fluctuation measures.

Appendix D The $a_{nm}$ coefficients

In this appendix we discuss the partial covariance for the case, where the $\overline{C}$-correlation function for hadrons (with autocorrelations removed) is expressed via the expansion Bzdak and Teaney (2013)

[TABLE]

where $T_{n}(x)$ is a set of orthonormal functions and $[-Y,Y]$ is the pseudorapidity domain. The choice of Aad et al. (2015b, a); Jia et al. (2016) is

[TABLE]

where $P_{n}(x)$ denote the Legendre polynomials. The orthonormality condition $\int_{-1}^{1}dxT_{n}(x)T_{m}(x)=\delta_{nm}$ is satisfied. The $a_{nm}$ coefficients are

[TABLE]

The transformation to the partial $C$-correlation function of Eq. (15) leads, in general, to mixing of the $a_{nm}$ coefficients, i.e., the coefficients $a^{C}_{nm}$for $C(S_{F},S_{B}\bullet S_{C})$ become complicated functions of the $a_{nm}$ coefficients for $\overline{C}(N_{F},N_{B})$.

We can, however, derive a simple formula connecting these coefficients, introducing the expansion (for a fixed $\eta_{C}$ of the reference bin)

[TABLE]

Then

[TABLE]

(up to the rescaling effects of Eq. (10), not included explicitly in the present discussion).

We remark that for the special case of the BT model with the central reference bin, $a^{C}_{11}=a_{11}$.