Multipion Bose-Einstein correlations in pp, p-Pb, and Pb-Pb collisions at the LHC

Three- and four-pion Bose-Einstein correlations are presented in pp, p-Pb, and Pb-Pb collisions at the LHC. We compare our measured four-pion correlations to the expectation derived from two- and three-pion measurements. Such a comparison provides a method to search for coherent pion emission. We also present mixed-charge correlations in order to demonstrate the effectiveness of several analysis procedures such as Coulomb corrections. Same-charge four-pion correlations in pp and p-Pb appear consistent with the expectations from three-pion measurements. However, the presence of non-negligible background correlations in both systems prevent a conclusive statement. In Pb-Pb collisions, we observe a significant suppression of three- and four-pion Bose-Einstein correlations compared to expectations from two-pion measurements. There appears to be no centrality dependence of the suppression within the 0-50% centrality interval. The origin of the suppression is not clear. However, by postulating either coherent pion emission or large multibody Coulomb effects, the suppression may be explained.

Coherent emission is known to suppress Bose-Einstein correlations below the expectation from a fully chaotic particle emitting source. Some of the earliest attempts to search for coherence relied solely on fits to two-pion correlation functions [17]. The intercepts of the fits at zero relative momentum were found to be highly suppressed. However, it was quickly realized that Coulomb repulsion and long-lived emitters (e.g. long-lived resonance decays) also suppress the correlation function significantly. Furthermore, the precise shape of the freeze-out spacetime distribution is unknown. As a consequence, the corresponding functional form of the correlation function in momentum space is also unknown. Being such, there is no reliable way to extrapolate the measured correlation function to the unmeasured intercept.
Multipion Bose-Einstein correlations could provide an increased sensitivity to coherence as the expected suppression increases with the order of the correlation function [5,18,19]. However, the analysis of multipion Bose-Einstein correlations comes at the expense of increased complexity. Some of the earliest attempts to measure three-pion Bose-Einstein correlations relied on a different methodology and gave rather ambiguous results [20][21][22][23]. Recently the methodology of isolating three-and four-pion Bose-Einstein correlations has been considerably improved [19]particularly in regards to the treatment of long-lived pion emitters. Our previous measurements of three-pion correlations revealed a suppression which may arise from a coherent fraction (G) of 23% ± 8% at low p T at kinetic freeze-out [24].
We present three-and four-pion QS correlations in pp, p-Pb, and Pb-Pb collisions at the LHC measured with ALICE using the methodology presented in Ref. [19]. The QS correlations are extracted from the measured multipion distributions. The extraction of QS correlations relies on the treatment of long-lived pion emitters and final-state interactions (FSI), e.g. Coulomb correlations. QS correlations between pions separated by large distances (>∼ 100 fm) are only observable at very low relative momentum, where track merging effects and finite momentum resolution prevent reliable measurements. The effect of long-lived emitters at measurable relative momentum is to simply dilute the correlation functions. The presented correlation functions are corrected for this dilution as well as FSI and therefore should represent the pure QS correlations from short-lived pion emitters, i.e. the core of particle production. We also present the mixed-charge four-pion correlations, which are used to demonstrate the effectiveness of all corrections in the analysis procedure.
The measured multipion QS correlations require a reference in order to quantify a possible suppression. Lower order QS correlation functions form the reference in this analysis. Two-pion QS correlations, in particular, provide a direct measurement of the pair-exchange magnitudes, which may be used as a building block to form an expectation for higher order correlation functions. These "expected" multibody correlations were termed "built" in Ref. [19]. This article is organized into 7 sections. We explain the detector setup and data selection in Sec. 2. In Sec. 3, we describe the analysis methodology. The results are presented in Sec. 4. In Sec. 5, we discuss all of the systematic uncertainties investigated. We discuss several possible origins of the suppression in Sec. 6. Finally, in Sec. 7 we summarize our findings.

Experimental setup and data selection
Data from pp, p-Pb and Pb-Pb collisions at the LHC recorded with ALICE [25] are analyzed. The data for pp collisions at √ s = 7 TeV were taken during 2010, during 2013 for p-Pb collisions at √ s NN = 5.02 TeV, and during 2011 for Pb-Pb at √ s NN = 2.76 TeV.
The trigger conditions are slightly different for each of the three collision systems. For pp collisions, at least one hit in the Silicon Pixel Detector (SPD), at central rapidity, or either of the V0 detectors [26], at forward rapidity, is required. For Pb-Pb and p-Pb collisions, the trigger is formed by requiring hits simultaneously in each V0.
The Inner Tracking System (ITS) and Time Projection Chamber (TPC) located at mid-rapidity are used for particle tracking [27]. There are 6 layers of silicon detectors in the ITS: two silicon pixel, two silicon drift, and two silicon strip detectors. The ITS provides high spatial resolution for the position of the primary vertex. The TPC alone is used for momentum and charge determination of particles through the radius of curvature of the particles traversing a 0.5 T longitudinal magnetic field. The TPC additionally provides particle identification through the specific ionization energy loss (dE/dx). To ensure uniform tracking, the z-coordinate (along the beam-axis) of the primary vertex is required to be within a distance of 10 cm from the detector center.
Tracks with a transverse momentum of 0.16 < p T < 1.0 GeV/c and a pseudorapidity of |η| < 0.8 are retained in this analysis. To ensure good momentum resolution a minimum of 70 tracking points in the TPC are required. The measured energy loss (dE/dx) of particles traversing the TPC and the corresponding uncertainty (σ ) are used to select charged pions [28]. Charged tracks observed in the TPC are identified as pions if their dE/dx is within 2σ of the Bethe-Bloch expectation for pions while being more than 2σ away from the Bethe-Bloch expectation for kaons and protons. The pion purity in our sample is studied with the HIJING generator [29], folded with the ALICE acceptance. In the sample selected with the procedure described above, about 96% of the particles are expected to be pions.
The effects of track merging and splitting are minimized by rejecting track pairs whose spatial separation in the TPC is smaller than a threshold value [24]. For three-pion and four-pion correlations, each same-charge pair in the triplet and quadruplet is required to satisfy this condition. Oppositely charged pairs are not required to satisfy this cut as they curve in opposite directions in the solenoidal magnetic field and are therefore easily distinguished.
The low multiplicity events produced in pp and p-Pb collisions contain a non-negligible nonfemtoscopic background arising from mini-jets [30][31][32]. We reduce this background by retaining only high multiplicity events in pp and p-Pb. For pp and p-Pb collisions, we retain events with at least 10 and 15 reconstructed charged pions, respectively. The choice of these boundaries are chosen to provide sufficient statistics while reducing non-femtoscopic background correlations. The multiplicity cut selects events from the top 46% and 42% of the cross-sections, respectively. In Pb-Pb collisions, all non-femtoscopic backgrounds are negligible. We ana-lyze Pb-Pb data from the top 50% collision centrality in ten equally divided intervals. The collision centrality in Pb-Pb is determined using the charged-particle multiplicity in the V0 detectors [26]. Approximately 13, 52, and 34 million events are used for pp, p-Pb, and Pb-Pb collisions, respectively.

Analysis technique
We follow the techniques outlined in Ref. [19] for the extraction of multipion QS correlations and a possible coherent fraction. Several types of multipion correlation functions are presented: QS 4 , and c QS 4 . The full three-pion correlation is given by C QS 3 and the cumulant correlation is given by c QS 3 . Four types of four-pion correlations are defined: the full correlation, C QS 4 ; two types of partial cumulant correlations, a QS 4 and b QS 4 ; and the cumulant correlation, c QS 4 . The full three-pion same-charge correlation function contains both pair and triplet symmetrization sequences while the cumulant contains only the triplet symmetrization sequence. The full four-pion same-charge correlation function contains four sequences of symmetrizations: singlepair, double-pair, triplet, and quadruplet symmetrizations. Partial cumulants, denoted by a QS 4 (b QS 4 ), have single-pair (single-and double-pair) symmetrizations explicitly removed. The cumulant correlation, denoted by c QS 4 , represents an isolation of the quadruplet symmetrization sequence.
Two-pion correlations are extracted from two types of pair momentum distributions, N 1 (p 1 )N 1 (p 2 ) and N 2 (p 1 , p 2 ), where p i is the momentum of particle i. N 1 (p 1 )N 1 (p 2 ) is measured by sampling two pions from different events with similar characteristic multiplicity and z-coordinate collision vertex class. N 2 (p 1 , p 2 ) is measured by sampling both pions from the same event. Three-pion QS correlations are extracted from three types of triplet distributions Four-pion QS correlations are extracted from the following quadruplet distributions The distributions in Eqs. 1-8 are formed by sampling the appropriate number of particles from the same event and the rest from different events. The subscript for N represents the number of pions taken from the same event. We normalize the distributions in Eqs. 1-2 to the distribution in Eq. 3 at a suitably large invariant relative momentum, q i j = −(p i − p j ) µ (p i − p j ) µ . Likewise, the distributions in Eqs. 4-7 are normalized to the distribution in Eq. 8. The q i j interval is chosen to be far away from the region of significant QS and FSI correlations. The normalization interval is 0.15 < q i j < 0.2 GeV/c in Pb-Pb while being 0.9 < q i j < 1.2 GeV/c in pp and p-Pb due to the wider QS correlations in smaller collision systems. The distributions are all corrected for finite momentum resolution and muon contamination [24].
The two-, three-, and four-pion distributions (N QS n ) are extracted from the measured distributions (N n ) with the appropriate coefficients according to the "core-halo" prescription [33] of short-and long-lived emitters [34]. In the core-halo model, a fraction of particles ( f c ) originate within a small radius component of particle production (the core). The rest, 1 − f c , originate within a much larger halo radius. The fraction of pairs, triplets, and quadruplets from the core is then given by f 2 c , f 3 c , and f 4 c , respectively. The other possibilities of mixed core-halo compositions are treated as well in this analysis. Pairs of particles from the core of particle production are separated by sufficiently short distances such that their QS and FSI correlations are experimentally observable. Pairs with one or both particles from the halo effectively dilute the correlation functions as no significant QS and FSI correlations are expected. The coefficients that isolate the multipion QS distributions are determined from the f c parameter [19].
The f c parameter is often associated with √ λ , where λ parametrizes the correlation strength, which is usually determined from fits to two-particle Bose-Einstein correlations. However, due to the unknown functional form of two-pion correlation functions, the λ parameter, determined this way, is convoluted with the arbitrary choice of fitting functions (e.g. Gaussian fits to non-Gaussian correlation functions). A more accurate extraction of f c is done by fitting mixedcharge two-pion correlations instead [24]. The correlation between π + and π − is dominated by Coulomb and strong FSI for which the wave functions are well known [35]. Owing to the large pion Bohr radius, π + π − correlations are less sensitive to the detailed structure of the source and can be fit less ambiguously wrt π + π + correlations. As part of the long-lived emitters correspond to weak decays (secondaries), f c is also sensitive to the specific tracking algorithm's ability to discriminate primary from secondary tracks. The value, f c = 0.84 ± 0.03, was used in Ref. [24] as well as in this analysis.
The distinction between core and halo may depend on the characteristic sizes and the dynamics of the system. Pions from decays of mid-lived emitters, such as the K * , Σ * , ω, and η ′ constitute a special case where the effect of QS correlations with other pions can be smaller than that of Coulomb correlations. Therefore, one might expect a slightly smaller core fraction for QS compared to Coulomb interactions. The magnitude of the difference should mainly relate to the fraction of pions produced from decays of mid-lived resonances. The resulting difference, which we assume to be small, is addressed by varying f c as discussed in the section on systematic uncertainties.
The treatment of multibody FSI (Coulomb and strong) is done according to the generalized Riverside approximation [19,21,23,24,36] where the n body FSI correlation is treated as the product of each pair FSI correlation, The two-pion FSI factor of pair (i, j) is given by K 2 (q i j ) and is calculated by averaging the modulus square of the Coulomb and strong wave function over an assumed freeze-out distribution.
We use the THERMINATOR model of particle production as an estimate for the freeze-out distribution [37,38]. The pair product approach to three-pion FSI correlations was shown to be a good approximation to the full asymptotic wave function calculation [19,24]. In this article we present QS correlation functions which are corrected for FSI and for the dilution of long-lived emitters according to Eqs. 33 and 39 in Ref. [19].
All distributions and correlation functions are projected onto the 1D the Lorentz invariant relative momentum. For three-and four-pion correlations, the sum quadrature of pair invariant relative momenta is used: The p T dependence of the correlation functions is studied by further projecting onto the average transverse momenta for two-, three-, and four-pion correlations, respectively. We form two intervals of K T3 defined by 0. 16 GeV/c and the RMS of the p T distribution is 0.03 GeV/c. At high K T3 , p T is 0.34 GeV/c and the RMS is 0.03 GeV/c. The same values also closely describe the low and high K T4 interval at low Q 4 (0.045 < Q 4 < 0.06 GeV/c). We further note that the p T is very similar for each q interval in this analysis. For 0.16 < K T2 < 0.3 GeV/c, p T increases linearly by about 0.015 GeV/c in the interval 0.005 < q < 0.2 GeV/c.

Extracting the pair-exchange magnitudes
The building blocks of Bose-Einstein correlations are the pair-exchange magnitudes (T i j ) and the coherent fraction (G) in the absence of multipion phases [5,18,19,39]. Multipion phases are expected when the space-time point of maximum pion emission is momentum dependent. However, the relative momentum dependence of the effect was shown to be rather weak [39]. Assuming a value of G, the pair-exchange magnitudes can be used to build all higher orders of correlation functions. We define the expected or built correlation functions, E n (i), which represent the expectation of higher order (n) QS correlations using lower order (i < n) experimental measurements as an input. The equations to build E n are given in appendix A. We define two types of expected correlation functions: The pair-exchange magnitudes can be extracted directly from two-pion correlation functions, which forms our primary expectation in Pb-Pb collisions. The twopion correlations are tabulated in four dimensions during the first pass over the data in the longitudinally co-moving system (q out , q side , q long , K T2 ). The interval width of each relative momentum dimension is 5 MeV/c, while it is 50 MeV/c in the K T2 dimension. In the second pass over the data, the previously tabulated two-pion correlations are interpolated for each pion pair from mixed events. We interpolate between relative momentum bins with a cubic interpolator. A linear interpolation is used in between K T2 bins, where a more linear dependence of correlation strength is observed.

E 4 (3) and e 4 (3):
We also extract the pair-exchange magnitudes from fits to C QS 3 (E 4 (3)) and c QS 3 (e 4 (3)) in 3D (q 12 , q 13 , q 23 ). The fit is performed according to an Edgeworth parametrization [40] as shown in equation 20 of Ref. [19]. This 2 nd approach is more limited as the pair-exchange magnitudes are extracted from a 3D projection of a 9D function. Similar to the 1 st type of expected correlations, the pair-exchange magnitudes are obtained from the first pass over the data and input into the second pass.
For the case of partial coherence, we assume that the pair-exchange magnitude of the coherent source is identical to the chaotic one (e.g. same radii) which might be expected for DCC radiation [7]. The value of G may then be extracted by minimizing the χ 2 difference between measured and expected correlations for each Q 3 or Q 4 bin. One may extract G from either of the six same-charge channels: The primary channel of extraction is C QS 4 for reasons of statistical precision and sensitivity to coherent emission. We also extracted G with several other multipion correlations and is shown in a separate note [41]. In pp and p-Pb collisions, where non-negligible non-femtoscopic backgrounds exist, we only use the 2 nd build technique as three-pion correlations have a larger signal to background ratio [42].
Both build techniques were tested using data generated by the THERMINATOR model, including a known coherent fraction [19]. The E 4 (2) correlations were typically 3% smaller than the "measured correlations" in THERMINATOR. The bias is attributed to the finite 4D projection of the true 6D two-pion correlation function. We correct for this potential bias in a data-driven approach. The interpolated two-pion correlation function from the 4D projection is compared to the true two-pion correlation function for each q interval. The ratio of the two correlation functions (subtracting unity from each), forms our correction factor.

Results
We now present the results of three-and four-pion QS correlations in pp, p-Pb, and Pb-Pb collisions. All correlations are corrected for FSI and for the dilution of pions from long-lived emitters. Mixed-charge correlations are first presented to demonstrate the effectiveness of all corrections in the analysis. Fits to same-charge three-pion correlations, which allow us to construct E 4 (3) and e 4 (3), are then presented. The comparison of measured to expected samecharge correlations assuming the null hypothesis (G = 0) is then presented. Finally we present the same comparison with non-zero values of G. charge quadruplets, the residue seen with the cumulant characterizes the effectiveness of several procedures. The baseline of the cumulant in pp collisions is offset from unity by about 10% and is due to statistical fluctuations in the high q normalization region of our data sample. It is included in the systematic uncertainty. The mixed-charge cumulant residues seen in pp and p-Pb collisions are similar in magnitude as seen in Pb-Pb collisions. Note that the FSI correlations are larger in pp and p-Pb with respect to Pb-Pb collisions. Isolation of the cumulant correlation function, c QS 4 , is done by subtracting several distributions as shown in Eqs. 4-7 after correcting for FSI. By default, we also utilize the distributions of two interacting opposite charge pions, N 2 (−, +)N 1 (−)N 1 (−) and N 2 (−, +)N 1 (−)N 1 (+) for π − π − π − π + and π − π − π + π + , respectively. After correcting for finite momentum resolution, muon contamination, and FSI corrections, such distributions should be identical to N 4 1 in the absence of additional correlations. A small difference in c QS 4 is observed without the subtraction of such terms [41].

Fits to three-pion correlation functions
The 2 nd build technique relies on the extraction of the pair-exchange magnitudes from fits to three-pion correlations. We separately fit both the cumulant (c QS 3 ) and full (C QS 3 ) correlations with an Edgeworth parametrization in 3D (q 12 , q 13 , q 23 ). The three-pion correlations and fits are projected onto Q 3 for pp, p-Pb, and Pb-Pb collisions in Figs  Four-pion measured correlations are compared to the E 4 (2) expectations in Pb-Pb in Figs. 6(a)-6(b) for low and high K T4 . Similar to the three-pion case, we observe a Q 4 dependent suppression of measured compared to the expected correlations.

Extracting a possible coherent fraction
We now investigate the expected correlations with non-zero values of the coherent fraction, G, and compare them to the measured correlations in Pb-Pb. We use the expected correlations of the 1 st type to extract the coherent fraction from four-pion correlations. Owing mostly to limitations of the three-pion fitting procedure, we do not extract the coherent fraction with the 2 nd type. The isospin effect relevant for charged-particle coherent states is neglected in this analysis [4,7,43,44]. Figure 7 presents same-charge four-pion correlations in Pb-Pb versus Q 4 at low K T4 . We observe that the suppression can be partially explained assuming G = 32% which minimizes the χ 2 of the difference of the ratio from unity for Q 4 < 0.105 GeV/c. The χ 2 /DOF of the minimum is quite low, 0.34, and is due to the inclusion of high Q 4 data in the calculation and the rapidly decreasing QS correlation with Q 4 . In Fig. 8 we present same-charge three-pion correlations in Pb-Pb versus Q 3 at low K T3 . In contrast to the four-pion case, the value of G = 32% does not   Fig. 4(a). The average of the charge conjugated correlation functions is shown. satisfactorily explain the suppression.
We also studied the centrality dependence of the suppression in Pb-Pb. Figures 9(a) and 9(b) show the centrality dependence of the extracted coherent fraction for low and high K T4 . Within statistical and systematic uncertainties, the coherent fractions are consistent for each centrality interval. We also parametrized the coherent component as a point source as opposed to the equal radii assumption used by default. The point source approximation may be expected to be more appropriate for gluon or pion condensate formation. The extracted coherent fractions with the point source approximation are shown in a separate note [41].
Previously [24], the coherent fractions were extracted from the r 3 observable which is intended to isolate the phase of three-pion correlations [39,45]. In contrast to the previous analysis, we estimate G by averaging the suppression in several Q 3  the unmeasured intercept. This approach was chosen due to the largely flat relative momentum dependence of previous r 3 measurements [19,24]. The values of G are obtained by averaging the bin-by-bin values within 0.03 < Q 4 < 0.105 GeV/c. Furthermore, our past analysis did not employ interpolation corrections which are relevant for the expected correlations. Correcting for the interpolation biases is expected to lower r 3 [19].
We extracted coherent fractions in Pb-Pb using the expected correlations of the 1 st type. The expected correlations of the 2 nd type were shown in all three collision systems but are expected to be less accurate due to more limited dimensionality and the fitting procedure of three-pion correlations. Being such, we could not reliably extract a value of G with the 2 nd build technique. The 2 nd type is, however, preferred in low multiplicity events, where non-negligible background correlations exist.
One of the most commonly cited sources of coherent pion emission is the DCC [8,10], which may occur as a consequence of chiral symmetry restoration. The most common prediction of the DCC is the fluctuation of charged to neutral pion production at low p T . If a single DCC domain is created within each event, we may expect a surplus of coherent charged pions in one event, while in another event, only coherent neutral pions are present. We investigated this possibility by first isolating a narrow multiplicity class at higher p T , 0.35 < p T < 0.5 GeV/c, within the 0-5% centrality class determined with the V0 detectors. From the multiplicity distribution of charged pions at the higher p T interval, we retain events which were within 1 standard deviation from the mean of the distribution. We then analyzed the multiplicity distribution of charged pions at low p T , 0.16 < p T < 0.25 GeV/c. Events with low p T multiplicities below the mean of the distribution were stored separately from those events above the mean. We do not observe a significant change of the suppression for events below or above the mean. The finding disfavors single-domain DCCs but does not rule out multidomain DCCs, for which independently coherent charged and neutral pions may be found in a single event [8,10].

Systematic uncertainties
We consider several sources of systematic uncertainty pertaining to the methodology and finite detector resolution. Below we describe each systematic uncertainty studied in order of decreasing magnitude. Some systematic uncertainties apply to only measured or expected correlations while others apply to both. The given values of the uncertainties apply to four-pion correlations. The values for three-pion correlations are generally smaller.
1. f c scale. The fraction of pion tracks from short-lived emitters for which QS and FSI correlations are experimentally observable is quantified with the f c parameter. From previous studies in ALICE using fits to π + π − FSI correlations, we estimate that f c = 0.84 ± 0.03 [24]. We vary f c within its uncertainties from the previous analysis. The uncertainty derived from varying f c applies to both measured and expected correlations and is about 6% at low 2. FSI variation. The default two-pion FSI correlation K 2 , together with the default value f c = 0.84, gives a satisfactory description of π + π − correlations [24]. We find that increasing the FSI correlation strength, |K 2 − 1|, by 5% while decreasing f c to 0.806 also provides a satisfactory description of π + π − correlations. The analysis was redone with such modifications, and the ratio of measured to expected four-pion correlations changed by less than 0.5%.
3. T i j extraction at high q. The 1 st type of expected correlations use the pair-exchange magnitudes (T i j ) extracted from two-pion correlations. The extraction of T i j becomes problematic at large q, where the measured two-pion QS correlations fluctuate beneath the baseline due to finite statistics. For such bins we set T i j = 0. We also constructed a separate expected correlation where the entire triplet or quadruplet was skipped if any pair T i j was negative. Half of the difference between these two builds was assigned as an uncertainty which is about 4% at high Q 4 and less than 0.1% at low Q 4 .

4.
Interpolation. We apply a data-driven approach to correct for interpolation biases, as already mentioned. From studies with different interpolation schemes, we find a 1% systematic uncertainty on the expected correlations at low Q 4 .

5.
Mid-lived emitters. The extraction of the multipion QS correlations from the measured distributions in Eqs. 4-8 relies on the f 41 , f 42 , f 43 , f 44 coefficients in Ref. [19]. The default values were derived in the "core-halo" picture of particle production, for which there are only short and long-lived emitters. In general there are also mid-lived emitters (e.g. ω decays) which modify the f coefficients and can be estimated using the THERMINATOR model. The effect was found to be quite small [19] and leads to a 0.5% uncertainty at high Q 4 .

Renormalization.
To account for small normalization differences between two-, three-, and four-pion correlation functions, the expected correlations are re-normalized to the ones measured at high 7. Detector resolution. Numerous effects related to finite detector resolution were checked. The charge conjugated correlation functions were consistent within statistical uncertainties. Similarly, the polarity of the solenoidal magnetic field had a negligible effect on the correlation functions. We compared Pb-Pb data from two different data-taking periods which were known to have different tracking efficiencies. The measured and expected correlation functions differed by less than 0.5%. Finite momentum resolution is known to smear the correlation functions, decreasing the correlation strength at low relative momentum for all orders of correlation functions. We correct for finite momentum resolution using HIJING (Pb-Pb) and PYTHIA [46] (pp and p-Pb) data simulated with the ALICE detector response. The uncertainty on the momentum resolution at low p T is governed by the material budget uncertainty of the ALICE detector and is estimated to be less than 10%. The corresponding uncertainty on the measured and expected correlations is about 1%. Our pion purity is estimated to be about 96% for which the remaining 4% impurity is dominated by muon contamination. Simulations have shown that most of the muons in our sample originate from charged pion decays for which QS and FSI correlations are expected with primary pions. We apply muon corrections similar to Refs. [24,42]. We assign a 2% uncertainty to the muon correction procedure. The tracking efficiency of the ALICE detector decreases rapidly for p T < 0.2 GeV/c [28]. To estimate the potential bias caused by the tracking efficiency, we randomly discard pions in THERMINATOR according to the TPC reconstruction efficiency. We do not observe a bias on the measured nor expected correlation functions which could cause an artificial suppression.
In addition to the above mentioned sources of systematics, we also applied an additional uncertainty to cumulant correlation functions, c QS 4 . The cumulant correlations were found to be much more sensitive to effects induced by low statistics at low Q 4 . The additional uncertainty is several tens of percents for the lowest Q 4 bin.
Most of the systematic uncertainties were found to be similar in magnitude and highly correlated for both measured and expected correlations. As a consequence, the systematics largely cancel in the ratio of measured to expected. For the ratio, we apply the maximum difference of measured and expected systematics. The systematic uncertainties for the ratio are dominated by the interpolator and mid-lived emitter uncertainty at low Q 4 . At high Q 4 , the muon corrections and the extraction of T i j at high q dominate the uncertainties.

Possible origins of the suppression
A suppression of three-and four-pion Bose-Einstein correlations compared to the expectations from two-pion measurements has been observed in Pb-Pb collisions. Below we list our investigations into the origin of the suppression.
1. Quantum coherence. Incorporating the effects of quantum coherence can perhaps explain the four-pion suppression in Fig. 7 with a centrality averaged coherent fraction of 32% ± 3%(stat) ±9%(syst). However, the same coherent fraction fails to explain the suppression at the three-pion level in Fig. 8. In particular, the suppression at the lowest Q 3 and Q 4 intervals cannot be resolved with the same coherent fraction as needed at higher Q 3 and Q 4 intervals. The isospin effect for charged-pion coherent states [4,7,43,44] has not been calculated, since the expressions which incorporate isospin conservation do not exist at the four-pion level. For G = 32%, the isospin effect increases the intercept of two-and three-pion correlations by about 1% and 3%, respectively. The effect on the expected correlations at finite relative momentum has not been calculated.
2. Coulomb repulsion. Same-charge pions experience Coulomb and strong repulsion which is stronger for quadruplets than for pairs. The four-pion Coulomb corrections used in this analysis correspond to the asymptotic limit of the Coulomb wave function as mentioned before. Previous studies [47] have justified the use of such wave functions for the characteristic freeze-out volumes and relative momenta studied in this analysis. We have also shown that the cumulant (c QS 4 ) of mixed-charge correlations are near unity after FSI − − + + + + − − −+ − − ++ + + ++ Low K T3 ,K T4 1.07 ± 0.01 1.16 ± 0.02 1.6 ± 0.1 0.89 ± 0.02 1.17 ± 0.02 High K T3 ,K T4 1.06 ± 0.01 1.13 ± 0.02 1.2 ± 0.1 0.89 ± 0.02 1.09 ± 0.02 Table 1: The x factors used to modify the multipion FSI factor such that the suppression of same-charge correlations and the residues of mixed-charged cumulants are resolved. The multipion FSI factor is modified according to: K 3,4 → x|K 3,4 − 1| + 1. With − − −+ correlations, only K 4 was modified and not K 3 which is also used to isolate the cumulant. We further note that x is more Q 3 and Q 4 dependent for the case of + + + and + + ++. We find that for the lowest Q 3 and Q 4 bin, x is about 1.2. corrections. In the case that the genuine multipion Coulomb interactions are not negligible, we modify the three-and four-pion FSI correlations by an amount, x, needed to resolve the suppression (residue) of same-charge (mixed-charge) correlations. The FSI factors are modified as: K 3,4 → x|K 3,4 − 1| + 1. The x factors given in Tab. 1 demonstrate that if the suppression is solely caused by genuine multipion Coulomb effects, they should modify the two-body approximation by up to 20% at low relative momentum for the case of same-charge three-and four-pion correlations. Such large multibody Coulomb correlations are not expected from the arguments provided in Ref. [47]. 4. Background correlations. Event generators such as HIJING and AMPT [48] do not include the effects of QS nor FSI and may thus be used to estimate background correlations. We checked two-, three-, and four-pion correlation functions in the 5% most central events from HIJING and AMPT. All orders of correlation functions were consistent with unity.
5. Multipion phases. The expected correlations ignore the three-and four-pion Fourier transform phases [39]. The r 3 observable was extracted in ALICE [24] and THERMINATOR [19] and no significant Q 3 dependence was found. As the trend of G with Q 3 and Q 4 is opposite to that expected from the phases [41], we find them unlikely to explain the suppression.
6. Multipion distortions. At high freeze-out phase-space density, all higher order symmetrizations, which are usually neglected, can contribute significantly to all orders of correlation functions [49][50][51][52][53]. The distortions have been calculated for two-pion correlations and recently for three-and four-pion correlations [19]. The calculations suggest that the ratio of measured to expected correlations is robust with respect to this effect.

Summary
Three-and four-pion QS correlations have been measured in pp, p-Pb, and Pb-Pb collisions at the LHC. The measured same-charge multipion correlations are compared to the expectation from lower order experimental correlation functions. A significant suppression of multipion Bose-Einstein correlations has been observed in Pb-Pb collisions. The ratio of measured to expected same-charge four-pion correlations is about 6σ below unity in our lowest Q 4 interval.
In pp and p-Pb collisions, owing to background correlations at low multiplicity in two-pion correlation functions, we compare the measured four-pion correlations to the expectation from fits to three-pion correlations (E 4 (3) and e 4 (3)). Three-pion correlation functions contain substantially larger QS correlations and reduced background correlations, which makes them a preferred base for higher order expectations in pp and p-Pb collisions. We do not observe a significant suppression of four-pion correlations in pp nor p-Pb collisions. However, the more limited dimensionality and fitting procedure to three-pion correlations makes E 4 (3) and e 4 (3) expectations less accurate than E 4 (2). Nevertheless, despite the presence of the nonfemtoscopic background, we also performed the analysis in pp and p-Pb collisions with the first type of expected correlations (E 4 (2) and E 3 (2)). No significant suppression was observed in pp or p-Pb collisions, although the unknown strength of the non-femtoscopic background prevents an absolute statement.
Mixed-charge four-pion correlations have also been measured. They are used to demonstrate the effectiveness of the cumulant isolation via the event-mixing techniques as well as that of the FSI, muon, and momentum resolution corrections. The mixed-charge cumulant correlations are shown to be near unity although a finite residue exists with both types of mixed-charge correlations.
The suppression of same-charge three-and four-pion correlations in Fig. 7 and 8 cannot be unambiguously resolved with any of the possible origins discussed. For example, if genuine multipion Coulomb interactions are non negligible, a large increase of as much as 20% beyond the two-body approximation would be needed to account for the observed suppression. On the other hand, a coherent fraction of about 32% ± 3%(stat) ±9%(syst) could largely explain the four-pion suppression, but the same value cannot explain the three-pion suppression. There does not appear to be a significant centrality dependence to the extracted coherent fractions. The weak K T2 dependence of the coherent fractions does not favor the formation of Bose-Einstein condensates nor disoriented chiral condensates, which are expected to radiate mostly at low p T . The suppression observed in this analysis appears to extend at least up to p T ∼ 340 MeV/c.