Two-particle Bose–Einstein correlations in pp collisions at √ s = 13 TeV measured with the ATLAS detector at the LHC

ThispaperpresentsstudiesofBose–Einsteincor-relations (BEC) in proton–proton collisions at a centre-of-mass energy of 13 TeV, using data from the ATLAS detector at the CERN Large Hadron Collider. Data were collected in a special low-luminosity conﬁguration with a minimum-bias trigger and a high-multiplicity track trigger, accumu-lating integrated luminosities of 151 µ b − 1 and 8 . 4 nb − 1 , respectively. The BEC are measured for pairs of like-sign charged particles, each with | η | < 2 . 5, for two kinematic ranges: the ﬁrst with particle p T > 100 MeV and the second with particle p T > 500 MeV. The BEC parameters, charac-terizing the source radius and particle correlation strength, are investigated as functions of charged-particle multiplicity (up to 300) and average transverse momentum of the pair (up to 1 . 5 GeV). The double-differential dependence on charged-particle multiplicity and average transverse momentum of the pair is also studied. The BEC radius is found to be independent of the charged-particle multiplicity for high charged-particle multiplicity (above 100), conﬁrming a previous observation at lower energy. This saturation occurs independent of the transverse momentum of the pair. quadrature, systematic uncertainties include track efﬁciency, Coulomb correction, non-closure and multiplicity-unfolding uncertainties


Introduction
Particle correlations play an important role in the understanding of multiparticle production in hadronhadron, hadron-nucleus and nucleus-nucleus collisions. Correlations between two identical bosons, called Bose-Einstein correlations, (BEC), are a well known phenomenon in high-energy and nuclear physics (for reviews see Refs. [1][2][3][4][5][6][7][8]). The BEC are often considered to be the analogue of the Hanbury Brown and Twiss effect [9][10][11] in astronomy, describing the interference of incoherently emitted identical bosons [12][13][14][15]. The BEC constitute a sensitive probe of the space-time geometry of the hadronization region, and allow the determination of the size and shape of the source from which particles are emitted. High-multiplicity data in proton-proton interactions can serve as a reference for studies of nucleus-nucleus collisions [16]. The effect of BEC is reproduced in both the hydrodynamical/hydrokinetic [17][18][19] and the Pomeron-based [20,21] approaches for modelling of high-multiplicity hadronic interactions. The BEC between a pair of particles depend on the particle multiplicity in the event and on the average transverse momentum of the pair. The dependence on these parameters is particularly sensitive to the space-time features of the hadronization region, as discussed by the authors of various models [22][23][24][25][26][27].
The BEC can be studied for two, three, or more identical bosons, assuming a hadronization region parametrized by one, two, or three size parameters. At the Large Hadron Collider (LHC), the BEC with one size parameter, the source radius, have been studied by the ATLAS Collaboration in collisions at the centre-of-mass energies √ = 0.9 and 7 TeV [28], and in +Pb collisions at √ NN = 5.02 TeV per nucleon-nucleon pair [29]. The BEC source radius parameter, , was observed by ATLAS to saturate, i.e. to reach a plateau level, at high charged-particle multiplicity, in √ = 7 TeV collisions [28]. The BEC in collisions have also been studied in the ALICE [30][31][32][33][34], CMS [35][36][37][38] and LHCb [39] experiments. The CMS results at 13 TeV also show a plateau in at high multiplicity [38].
Pomeron-based models predict a plateauing of at intermediate multiplicities as a consequence of the overlap of the colliding protons, followed by the return to a lower value of at high multiplicities [20,21]. The IP-Glasma model [40] predicts some tendency to saturation at high multiplicities for interactions.
In a different approach, two-particle correlations have been analysed at √ = 7 TeV by ATLAS in the framework of a phenomenological model based on ordered hadron chains. This approach provides an alternative interpretation of the correlations commonly attributed to Bose-Einstein interference [41].
In this paper, studies of the BEC in collisions at √ = 13 TeV using the ATLAS detector [42] are presented. The same methodology is used in this analysis as was used in the previous ATLAS BEC studies [28], and the same minimum-bias (MB) trigger is used as was used in an earlier study of charged-particle distributions at 13 TeV [43,44]. In addition, the study is extended to the region of high-multiplicity collisions by using the high-multiplicity track (HMT) trigger previously used for studies of long-range elliptic azimuthal anisotropies [45]. This paper is organized as follows. The analysis method is described in Section 2. The ATLAS detector is described in Section 3. The Monte Carlo (MC) and data samples are presented in Section 4. The track selection criteria and data correction procedure are presented in Section 5 and the reference sample is discussed in Section 6. The systematic uncertainties are presented in Section 7. The BEC source radius and correlation strength dependence on, separately, charged-particle multiplicity and average transverse momentum of the pair is shown in Section 8. This section also presents the double-differential dependence of the BEC parameters on multiplicity and average transverse momentum. A comparison with previous, lower-energy ATLAS results is presented, as well as a comparison with CMS results at 13 TeV. The summary and conclusions are given in Section 9.

Analysis method
The BEC are observed at small values of the square of the four-momentum difference, 2 , between two identical bosons, where 2 = −( 1 − 2 ) 2 = 2 12 − 4 2 , 1 and 2 are the four-momenta of the particles, 12 is the invariant mass of the particle pair and is the particle mass. The BEC parameters are measured in terms of a two-particle correlation function, where ( ) is the distribution of formed from all like-sign charged-particle pairs, assuming in this analysis that all charged particles are pions, and 0 ( ) is a reference distribution specially constructed to exclude the effect of BEC. In this analysis, 0 ( ) is constructed from unlike-sign charged-particle pairs, as discussed in Section 6. The distributions ( ) and 0 ( ) are normalized to unity, i.e. they are probability density functions.
To account for the effects of resonances (for example , , , , * , ), the 2 ( ) correlation function is corrected by dividing by the corresponding quantity constructed from MC simulations that contain resonances but not the BEC: where data 2 ( ) and MC 2 ( ) are the 2 ( ) functions reconstructed from data and MC samples, respectively.
The BEC effect is usually described by a function, Ω, with two parameters: the effective parameter, and the strength parameter, , where the latter is also called the incoherence or chaoticity parameter [46]. If measurement can be made down to a very small value of ( 10 MeV), the parameter may also be sensitive to a "halo" contribution due to long-lived resonances [5], but this is not possible in this analysis. A typical functional form is: In a simplified scheme for fully coherent emission of identical bosons, = 0, while for incoherent (chaotic) emission, = 1. Various parameterizations of the Ω( , ) function in Eq. (3) can be found in the literature (see for example Refs. [3,47]), each assuming a different shape for the particle-emitting source. The parameter 0 is a constant depending on the normalizations of ( ) distributions and is chosen such that 2 ( ) is unity at large .
A Lévy-type parameterization for Ω introduces a further free parameter, : with 0 < ≤ 2. The Lévy distribution is a simple exponential when = 1 and a Gaussian when = 2. Fits to 2 ( ) with free usually return a value between 1 and 2, and the fit quality, in terms of 2 / , is generally better than the fit using an exponential or Gaussian. However, the interpretation of the R parameter in terms of the distribution of space-time points is then less clear. Also, in studies for the 7 TeV analysis it was found that varied when 2 ( ) was studied more differentially in terms of charged multiplicity or pair transverse momentum. So, in common with the other LHC experiments, this analysis constrains to be 1 or 2. Investigation of the Gaussian parameterization is found to give a significantly poorer description of the data in the low-region. (See Section 8.1). So, in the study presented here, the BEC parameters are extracted using an exponential parameterization, of a static source assuming a radial Lorentzian distribution of the source. The fitted parameter in Eq. (3) takes into account possible long-distance correlations that have not been removed fully by ( ).

The ATLAS detector
ATLAS is a multipurpose particle physics experiment operating at one of the beam interaction points at the LHC. The ATLAS detector 1 covers almost the whole solid angle around the collision point with layers of tracking detectors, calorimeters and muon chambers. It is designed to study a wide range of physics topics at LHC energies. For the measurements presented in this paper, the tracking devices and the trigger system are of particular importance.
The innermost part of the ATLAS detector is the inner tracking detector (ID), which has full coverage in and covers the pseudorapidity range | | < 2.5. The ID is immersed in the 2 T axial magnetic field of a superconducting solenoid and measures the trajectories of charged particles. It consists of a silicon pixel detector (Pixel), a silicon microstrip detector (SCT) and a straw-tube transition radiation tracker (TRT), each split into a barrel and two endcap components. The Pixel, SCT and TRT are located around the interaction point spanning radial distances of 33-150 mm, 299-560 mm and 563-1066 mm, respectively. The barrel (each endcap) consists of four (three) pixel layers, four (nine) double layers of silicon microstrips and 73 (160)  For Run 2 a new innermost pixel layer, the insertable B-layer (IBL) [48], was installed around a new, narrower and thinner beam pipe. The IBL is composed of 14 lightweight staves arranged in a cylindrical geometry, each made of 12 silicon planar sensors in its central region and 2 × 4 three-dimensional sensors at the ends. The IBL pixel dimensions are 50 m in the -direction and 250 m in the -direction (compared with 50 m by 400 m for the other pixel layers). The intrinsic spatial resolution of the IBL readout is 10 m in ( , )-position and 75 m in -position [49]. The smaller radius and the reduced pixel size result in improvements in both the transverse and longitudinal impact parameter resolutions [43,44]. The services for the existing pixel detector were upgraded, significantly reducing the amount of material in the region | | > 1.5, in particular at the boundaries of the active tracking volume.
A track from a charged particle traversing the barrel detector typically has 12 silicon measurement points (hits), of which 4 are Pixel and 8 are SCT, and more than 30 TRT straw hits. Requirements on an IBL hit and on impact parameters strongly suppress the number of tracks from secondary particles.
The ATLAS detector has a two-level trigger system: the first-level (L1) trigger and the high-level trigger (HLT) [50]. Events used in this analysis were required to satisfy L1 triggers using the minimum-bias trigger scintillators (MBTS). These are mounted at each end of the detector in front of the liquid-argon endcapcalorimeter cryostats at = ±3.56 m, and are segmented into two rings in pseudorapidity (2.07 < | | < 2.76 and 2.76 < | | < 3.86) with 8 azimuthal sectors in the inner ring and 4 in the outer. The MB events were selected on the basis of the MBTS alone, while the HMT events satisfied additional HLT criteria. Further details are given in Section 4.1.
An extensive software suite [51] is used in the reconstruction and analysis of real and simulated data, in detector operations and in the trigger and data acquisition systems of the experiment.

Data samples and event selection
This study uses the -collision datasets at √ = 13 TeV that were used in previous ATLAS studies of MB events [43,44,52] and HMT events [53]. The MB and HMT datasets were taken in a special configuration of the LHC in June 2015, with low beam currents and reduced beam focusing, producing a low mean number of interactions per bunch-crossing, : in the range 0.003-0.007 for MB events and 0.01-0.04 for HMT events. This configuration and the event selection used are documented in detail in the ATLAS √ = 13 TeV MB analysis [43,44]. For the MB data sample, events were selected from colliding proton bunches using the MBTS trigger. This trigger selection required one counter above threshold from either side of the detector, and is referred to as a single-arm trigger. The HMT trigger is a combination of an L1 requirement on the number of hits in the MBTS and an HLT requirement on the multiplicity of HLT-reconstructed charged-particle tracks. This trigger required more than 900 SCT space-points and more than 60 reconstructed good-quality charged tracks with T > 400 MeV associated with the reconstructed vertex with the highest multiplicity in the event. These selections correspond to integrated luminosities of 151 b −1 and 8.4 nb −1 for the MB and HMT triggers, respectively.
In the offline selection, track candidates are reconstructed in the Pixel and SCT detectors and extended to include measurements in the TRT [54]. The selected tracks satisfy the following criteria: T > 100 MeV and | | < 2.5; at least one Pixel hit and an IBL hit if expected; 2 at least two, four or six SCT hits for T < 300 MeV, T < 400 MeV or T > 400 MeV, respectively, in order to account for the dependence of track length on T . All events are required to have at least one vertex [55], formed from a minimum of two tracks. The vertex position is required to be consistent with the beam-spot position within the ATLAS detector [54]. If more than one vertex is reconstructed, the primary vertex (PV) is defined as the vertex with the largest number of associated reconstructed tracks, weighted by the square of their transverse momenta [54]. Selected tracks are then also required to satisfy | BL 0 | < 1.5 mm, where the transverse impact parameter BL 0 is calculated with respect to the measured beam line (BL) [56]; and | BL 0 × sin | < 1.5 mm, where BL 0 is the difference between the longitudinal coordinate of the track along the beam line at the point where BL 0 is measured and the longitudinal coordinate of the PV, and is the polar angle of the track. High-momentum tracks with mismeasured T are removed by requiring the track-fit 2 probability to be larger than 0.01 for tracks with T > 10 GeV. A special configuration of the track reconstruction algorithms is used in this analysis to reconstruct low-momentum tracks [54]. The 2 A hit is expected if the extrapolated track crosses a known active region of a pixel module. If an innermost pixel layer hit is not expected, a next-to-innermost pixel layer hit is required if expected.
track reconstruction efficiency at T ≈ 100 MeV is about 28% and increases rapidly to almost 75% at T ≈ 300 MeV [44]. For the MB trigger, the low-T selected dataset consists of 9.3 × 10 6 events with a total of 2.4 × 10 8 tracks with T > 100 MeV. The high-T selection requires in addition at least two tracks with T > 500 MeV. This yields 8.9 × 10 6 events with 1.1 × 10 8 tracks with T > 500 MeV. For the HMT trigger, which requires at least 60 tracks with T > 400 MeV at the trigger level, the low-T selection requires at least 60 tracks with T > 100 MeV and results in 9.1 × 10 6 events with 9.8 × 10 8 tracks with T > 100 MeV. For T > 500 MeV the requirement is at least 40 tracks with T > 500 MeV; the sample consists of 9.1 × 10 6 events with 5.3 × 10 8 tracks with T > 500 MeV.
To reduce contamination from multiple proton-proton interactions (pile-up) in the colliding bunches, events with a second vertex containing four or more tracks are removed. Such events represent less than 0.3% of the sample. Secondary vertices can be reconstructed with high efficiency provided the distance, Δ , between the -coordinates of the primary and pile-up vertices satisfies |Δ | ≥ 4 mm. To estimate any possible residual influence of multiple interactions in MB-and HMT-triggered data, a dedicated study of the Δ distribution was performed, using a sample with no pile-up rejection. This distribution was extrapolated to estimate the level of pile-up that will survive the maximum track-multiplicity limit on a second vertex. For the MB events, the fraction of events with pile-up vertices in the region |Δ | ≤ 4 mm around the PV is about 0.006%. For the HMT events, it is about 0.05%. 3 The average track multiplicities per pile-up vertex for the MB and HMT samples are 9.4 and 23, respectively. For the primary vertices these multiplicities are 26 and 108. The fractions of surviving pile-up tracks in the MB and HMT samples are approximately 0.002% and 0.01% per event, respectively, and these levels have a negligible influence on the BEC studies.
The contributions from beam-gas collisions and from non-collision background (cosmic rays and detector noise) were investigated in Ref. [57] and found to be negligible.

Monte Carlo samples
The P 8 MC generator [58] was used to calculate the MC 2 ( ) correlation functions (Eq. (2)). The generator was used with parameter values set to the A2 tune [59] and with the MSTW2008 PDF set [60]. P 8 models the effects of colour coherence, which is important in dense parton environments and effectively reduces the number of particles produced in multiple parton-parton interactions. In P 8, the simulation is split into non-diffractive and diffractive processes. The non-diffractive process is dominated by -channel gluon exchange, comprising approximately 80% of the selected events. The diffractive processes, which are further split into single-and double-diffractive processes, are described by a Pomeron-based approach [61]. The contributions from non-diffractive, single-diffractive and double-diffractive processes were included in proportion to the cross sections predicted by P 8 with the A2 tune. Large MC samples of MB and HMT events were generated and passed through the ATLAS simulation program [62], which is based on G 4 [63], and the reconstruction chain, which is exactly that used for collision dataset. In addition, all selection criteria applied to the data were also applied as appropriate to the MC samples.
In order to estimate the systematic uncertainty due to the choice of MC model, large MC samples, comparable to the size of the collision dataset, were also generated with three other programs. These were the P 8 Monash tune [64], the EPOS MC generator [EPOS] and H ++ [65] using the UE-EE-5 3 Since is much smaller than 1 for both the MB and HMT samples, the ratio of the number of events with more than one vertex to the number of events with at least one vertex is ≈ /2. The requirement |Δ | ≤ 4 mm reduces this fraction by a further factor of ≈ 100.
tune [66]. The Monash underlying-event tune is based on the NNPDF2.3 PDF [67] and incorporates updated fragmentation parameters. EPOS implements a parton-based Gribov-Regge [68] theory, which is an effective field theory describing both hard and soft scattering. The EPOS generator was used with the LHC tune [69]. The H ++ UE-EE-5 tune is based on the CTEQ6L1 PDF. All sets of generated events were passed through the ATLAS simulation program and reconstruction chain.
The purity of the analysis sample of particle pairs is defined as the fraction of identical particle pairs ( ± ± , ± ± , and¯¯) in the sample relative to all like-charge particle (LCP) pairs. The purity of the analysis sample in terms of identical boson pairs is estimated from P 8 with the A2 tune to be about 74% (where about 73% are ± ± and about 1% are ± ± ), using the selection criteria described above for MB T > 100 MeV. In studies for the 7 TeV analysis [28] it was found that the impurity reduces the value of by a factor that is approximately equal to the impurity of the sample, where a 100% purity has no effect, but the parameter is unaffected.

Data correction procedure
To correct for the inefficiencies in the trigger and in the reconstruction of tracks and vertices, event-by-event and track-by-track weighting procedures are applied as appropriate and as used in the previous ATLAS measurements [43,44].
The efficiency of the MB trigger was studied with an independent control trigger [43,44]. The control trigger selects events randomly at L1, which are then filtered in the HLT by requiring at least one reconstructed track with T > 200 MeV. The MB trigger efficiency increases from 97% for events with low track multiplicity, sel = 2, to 100% for events with higher multiplicity, sel ≥ 7, where sel is the number of tracks reconstructed offline and satisfying the selection criteria given in Section 4.1. The efficiency of the HMT trigger is studied using MB-triggered events as the control sample. The efficiency is found to be 15% for events with sel 85 and reaches a plateau at 100% for sel 110. The measured MB and HMT trigger inefficiencies as a function of sel are used to correct each event. These corrections are found to have a negligible impact on the extraction of the BEC parameters discussed in Section 8. The effect of the vertex-reconstruction efficiency on the BEC parameters is also negligible.
In addition, each track is assigned a weight that corrects for the track reconstruction efficiency, for the fraction of secondary particles, for the effect of detector resolution smearing primary particles 4 into and out of the fiducial region ( T > 100 MeV, | | < 2.5) and for the fraction of fake tracks formed by a random combination of hits [43,44]. The effect of the loss of reconstructed pairs at small opening angles is negligible in comparison with the track reconstruction corrections.
The resolution is found to be similar for all MC generators. For the BEC region < 0.2 GeV, where the distribution rises steeply, the resolution is better than 5 MeV and much smaller than the bin size. For > 0.2 GeV, where the distribution is essentially flat, the resolution increases linearly with , rising to 40 MeV for 2 GeV. Overall, no unfolding of is required.
The track reconstruction inefficiency is such that sel differs significantly from the true charged-particle multiplicity, ch . In order to present and study 2 ( ) (Eq. (2)) as a function of the true charged-particle multiplicity, ch , the correlation functions 2 ( , sel ) must be unfolded. In the earlier ATLAS MB study at 13 TeV [43, 44] an unfolding matrix that reflects the probability of reconstructing sel charged tracks in an event with a generated charged-particle multiplicity ch was constructed, using MC simulation with the P 8 ATLAS A2 tune. For the present analysis this unfolding matrix was applied to both the MB and HMT data. It is found that, on average, an event with sel = 100 reconstructed tracks corresponds to ch = 127.6 ± 0.1 charged particles at particle level.
As a check of the entire analysis chain, the distribution of 2 ( ) (Eq. (1)) using the generated particles, MC,gen 2 ( ), was compared with the MC,rec 2 ( ) distribution using the simulated detector response to reconstruct the tracks and vertices, which were then corrected for inefficiencies as described above, and unfolded to the particle level. This study was carried out for several ch and T intervals, where T is the average transverse momentum of the pair: In the small-region the difference between the unfolded and particle-level distributions varies, depending on ch and T , but is a few percent on average. The effect is significant for the first two -bins used, 0.02 ≤ < 0.04 GeV and 0.04 ≤ < 0.06 GeV, and for these bins it is accounted for by adding it in quadrature to the statistical uncertainty in the 2 ( ) distribution. It is therefore not included in the systematic uncertainties discussed in Section 7 and listed in Table 1. A large part of the difference in the first -bin, 0.02 ≤ < 0.04 GeV, is due to photon conversions in the detector material, which are included in the simulated detector response but are not part of the set of generated particles. Since 2 is a ratio of like-sign pairs to unlike-sign pairs, the denominator in MC,rec 2 ( ) contains the contribution of pairs from photon conversions, which show up at small , whereas MC,gen 2 ( ) does not. For MB events the difference for the first bin is 4%, of which 3.5% is the contribution of photon conversions. For HMT events, the corresponding numbers are 3.25% and 2%. Conversion photons are less significant for the HMT events, because the number of track combinations increases quadratically with ch while the number of conversion photons increases only linearly. The effect of photon conversions in the region > 0.04 GeV is negligible.
The long-range Coulomb force causes a momentum shift between particles in the LCP and unlike-charge particle (UCP) pairs. The density distributions are corrected for this effect by applying the Gamow penetration factor per track pair with a weight 1/ ( ) such that corr ( ) = ( )/ ( ) [70][71][72] (for a review see Ref. [73]). The Gamow factor is given by with the dimensionless parameter defined as = ± / , where is the electromagnetic fine-structure constant and is the pion mass. The corrected two-particle density function, corr ( ), is normalized to unity. The sign of is positive for LCP pairs and negative for UCP pairs. The Gamow factor particularly affects pairs with < 0.1 GeV, and the correction is 24% for 0.02 ≤ < 0.04 GeV, 14% for 0.04 ≤ < 0.06 GeV and 10% for 0.04 ≤ < 0.06 GeV. The Coulomb interaction is not included in the generated MC event samples, and so the correction is applied only to data events.

Reference sample
A good choice of reference sample is important in obtaining an unbiased estimate of the BEC signal. Ideally, 0 ( ) should include all momentum correlations present in the signal dataset except for those  Figure 1: Comparison of single-ratio two-particle correlation functions, data 2 ( ) and MC 2 ( ), with the two-particle double-ratio correlation function, 2 ( ), for the high-multiplicity track (HMT) events using (a) the opposite hemisphere (OHP) like-charge particles pairs reference sample for T -interval 1000 < T ≤ 1500 MeV and (b) the unlike-charge particle (UCP) pairs reference sample for T -interval 1000 < T ≤ 1500 MeV. arising from BEC. Several different methods of constructing an appropriate reference sample have been investigated by various collaborations and authors.
The main methods are: • take UCP pairs from the same event; • take LCP pairs from the same event, but invert the momentum of one of the selected particles: . This is the opposite hemisphere (OHP) technique [7]; • take LCP pairs from the same event, but invert the transverse momentum of one of the selected particles: ( x , y , z → − x , − y , z ). This is the -rotation method [7]; • take LCP pairs from different events, randomly mixed; this is the 'mixed event' (MIX) technique [74].
In addition, most of the recently used methods compare the data with a reference sample having global event parameter values that are the same as, or similar to, those in data -such as the event multiplicity, the invariant mass of the pair, or the rapidity of the pair. For the previous ATLAS analysis [28] the UCP, OHP and MIX procedures were investigated in detail using dedicated MC studies, and it was observed that MIX and OHP behaved similarly. For this analysis, the UCP and OHP were studied.
As the MC models do not implement BEC, the single-ratio correlation functions MC 2 ( ) are expected to be close to 1, and to be a rather flat function of , particularly in the low-, BEC-sensitive region. This is indeed the behaviour observed in the earlier analysis [28] over much of the parameter space in ch and T for the MIX or OHP technique, and for the OHP technique in this analysis. However, in high-multiplicity and high-T intervals, which are relevant for this analysis, MC 2 ( ) constructed with the MIX or OHP technique exhibits an increase in the low-, BEC-sensitive region.
As an example, Figure 1(a) shows the correlation functions using the OHP technique for the reference sample. The correlation functions shown are the single ratios data 2 ( ) and MC 2 ( ), and the double ratio 2 ( ), all for HMT events, and for the average transverse momentum of the pair 1000 < T ≤ 1500 MeV.
The MC 2 ( ) exhibits a rise similar to that seen in data 2 ( ) in the BEC-sensitive region, although no such effect should be present in MC 2 ( ). Such a rise would have a marked effect on the BEC parameters extracted from 2 ( ). Figure 1(a) also shows that 2 ( ) is less than 1 in the region 0.7 < < 1.2 GeV.
The use of UCP pairs for the reference sample was also investigated in detail. As in the other techniques, MC 2 ( ) behaves satisfactorily over much of the parameter space in ch and T , although resonances distort the region ≈ 0.5 GeV, particularly at low multiplicity or at high T . An example is shown in Figure 1(b), using the UCP reference sample, and for the same ch and T intervals as in Figure 1(a). Although there is some distortion in the resonance region, where MC 2 ( ) and data 2 ( ) do not match perfectly, the resultant 2 ( ) matches well the BEC region of data 2 ( ). Accordingly, the UCP method is used in this analysis, as it was for the previous ATLAS analysis [28], and the other methods are not suitable.
The UCP reference sample used in this analysis is taken from the same events used to form pairs of like-charge particles, i.e. 0 ( ) ≡ ( ± 1 , ∓ 2 ). By construction, this UCP sample has the same topology and global parameters as the LCP sample ( ± 1 , ± 2 ), but is naturally free of any BEC effect. However, the UCP sample contains oppositely charged hadron pairs from the decay of resonances such as , , , , * , and , which are not present in the LCP combinations. Details concerning resonances are given in Appendix A. These resonances contribute to the low-region, and can give a spurious BEC signature with a large effective [75][76][77][78][79][80][81]. The effect of the resonant states in the variable was checked in MC studies, and the most affected regions are excluded from the analysis of the BEC effect, as described in Section 8.1.

Systematic uncertainties
In this section the systematic uncertainties in values of the BEC parameters and are summarized. These parameter values are obtained from a fit to the exponential model (4). This fit is applied to the 'inclusive' distributions, corresponding to 2 ≤ ch ≤ 250 for MB and 100 < ch ≤ 300 for HMT, and all T . The fit is also applied to selected ranges of ch or T , and to the double-differential ( ch , T )-distributions. The discussion in this section is mainly in terms of the 'inclusive' fits, but the systematic uncertainties have been evaluated for all selections applied to the data.
The systematic uncertainties resulting from the track reconstruction efficiency, which are parameterized as a function of T and , were determined in earlier analyses [43,44] and these are implemented as uncertainties in the track weights for particle pairs in the distributions that enter the correlation functions. The effects of track splitting and merging have been studied in detail. They are largest for values below 0.04 GeV, but even for the first used bin (0.02 ≤ < 0.04 GeV) the effect is only 0.4% and is negligible in comparison with other uncertainties.
The effects of uncertainties in the unfolding matrix on the BEC parameters and are found to be negligible for ch < 200. For the highest-multiplicity MB (201 ≤ ch < 250) and HMT (231 ≤ ch < 300) events this systematic uncertainty is very small and not more than 0.01% for and 0.03% for .
The leading source of uncertainty in and is the choice of MC generator used in the calculation of 2 ( ).    As can be seen in this figure, there are some regions of ch and T where there are pronounced and significant differences between the results obtained using the different generators. In order to assign the systematic uncertainty due to the MC modelling, a conservative approach is adopted: for each generator the weighted deviation is evaluated with respect to P 8 A2, where the weight is inversely proportional  to the statistical uncertainty for that generator. The systematic uncertainty is taken as the largest (in magnitude) of these weighted deviations, treating positive and negative deviations separately, which leads to asymmetric uncertainties. This procedure is followed for MB and HMT events separately.
An important source of systematic uncertainty for the determination of the parameter arises from the Coulomb correction. This uncertainty is estimated by varying the Gamow factor ( ) by ±0.15, which is the magnitude of the effect on UCP or LCP pairs at = 0.02 GeV. This is the standard variation for a conservative estimate of the systematic uncertainty in the Coulomb correction [28,35], and accounts also for effects that include the extended size of the emission source, as discussed in Ref. [82]. The resulting systematic uncertainties for the inclusive fit of the parameter are the same for both triggers and are equal to 1.8%. For the parameter the resulting systematic uncertainties are small compared to the cumulative systematic uncertainties, and equal to 0.1% and 0.3% for MB and HMT events, respectively.
The effect of the upper limit of the fit range in is estimated by changing it from the nominal 2 GeV by ±0.1 GeV, which corresponds to a total variation of about five times the nominal -resolution for ≈ 2 GeV. This gives an estimate of the residual uncertainty due to long-range correlations accounted for by , the parameter in the linear term of Eq.  The same sources of uncertainty are considered for the single-and double-differential measurements in ch and T . The results for the ch distribution are also shown in Table 1. The pattern of uncertainties is similar to the inclusive case, but shows a wider spread. Away from the BEC region, the MC description of 2 improves as the multiplicity increases, and this is reflected in Table 1 by smaller systematic uncertainties for the HMT results. This is also true for the T distribution and the double-differential distribution.

Results
The BEC parameters, and , are extracted from a fit to the 2 ( ) distribution in intervals of ch and T . In addition, there are two possible selections on the T threshold of the final-state particles. In order to give a reasonably complete picture of all of this information, the parameters describing the dependencies of and on variables of interest are themselves fitted to simple functions (usually exponentials), and the parameters of these fits are presented.
In Section 8.1 the procedure used to fit the 2 ( ) distribution is described. The dependence on multiplicity is presented in Section 8.2. The comparison with previous ATLAS results at 7 TeV is discussed in Section 8.3. The behaviour with T is presented in Section 8.4. The 'double-differential' dependence of the BEC parameters on multiplicity and T at 13 TeV is explored in Section 8.5, by looking at the dependence on one variable for different intervals of the other. Finally, a comparison with CMS results at 13 TeV is shown in Section 8.6. Figure 3 shows examples of the 2 ( ) correlation functions measured at √ = 13 TeV for two different multiplicity intervals, and the fits to the Gaussian and exponential forms. In the lower panels of Figure 3 the differences between the measured 2 ( ) and the fit to the exponential form, normalized to the experimental uncertainty ( ), Δ 2 ( )/ ( ) = [ data 2 ( ) − fit 2 ( )]/ ( ), are shown for both multiplicity intervals. It is evident from Figure 3 that the Gaussian function does not give a good description of the low-, BEC-sensitive region, ≤ 0.1 GeV, while the exponential function matches the data well in this region. Therefore, the quantity Δ 2 ( )/ ( ) is not shown for the Gaussian case, and the Gaussian fit is not considered further. Figure 3 shows structure in 2 ( ), compared with the smooth fitting functions, in the regions ≈ 0.25 GeV and 0.4 ≤ ≤ 0.9 GeV. If the MC modelling of resonance production, particularly the , and , does not agree very well with data, then this region will be affected, as appears to be the case. Accordingly, the regions 0.2 ≤ ≤ 0.3 GeV and 0.4 ≤ ≤ 0.9 GeV are excluded from the fits to the 2 ( ) correlation functions for the MB and HMT events. For the MB events, similar arguments lead to the additional exclusion of the region 1.0 ≤ ≤ 1.16 GeV, as it appears to be influenced by the reflection of the 0 (980) and 2 (1270) mesons in the case of 100 ≤ T ≤ 200 MeV and in the multiplicity region ch ≤ 40. These excluded regions lie outside the BEC-sensitive region and the main effect of the exclusions is to markedly reduce the 2 without significantly changing the fitted BEC parameters. The results for the BEC parameters obtained using the exponential fit with the UCP reference sample for the MB and HMT samples, summed over all multiplicities for T > 100 MeV, are = In both cases the uncertainties are the statistical and systematic uncertainties summed in quadrature. The overall quality of the fits is not good, as the high number of events in the inclusive distributions exposes the imperfections in the models used. As noted in Section 2, the fits would be better using a Lévy distribution with an additional free parameter. For all fits throughout Section 8, the statistical uncertainty is multiplied by the factor √︁ 2 / if this factor is significantly greater than 1 [83]. This applies in relatively few cases, see e.g. 4 of the 8 entries in Table 3.

Multiplicity dependence
The dependence of and on ch is studied in multiplicity intervals that are chosen to be well populated, and comparable to those used by other LHC experiments [30,31,35,36,38].
To compare the BEC parameters from different T -thresholds, a scaled charged-particle multiplicity is introduced as follows. At 13 TeV, the average charged-particle multiplicity per unit of pseudorapidity in the region | | < 0.2, d ch /d | | |<0.2 , is 6.5 for particles with T > 100 MeV [44] and 3.24 for particles with T > 500 MeV [43]. A scaled charged-particle multiplicity for | | < 2.5, i. e. 5 units of pseudorapidity is then defined as: Using the exponential fit, the extracted BEC parameters ( ch ) and ( ch ) at 13 TeV are shown as a function of scaled multiplicity in Figure 4 for T > 100 and 500 MeV. The MB and HMT data agree well where they overlap.
For T > 100 MeV the parameter decreases with increasing ch , and the dependence is well described by the exponential function: The parameter increases with multiplicity up to about ch 3.08. For ch ≤ 1.9, the ch dependence of is well described by: which is also the function used in heavy-ion studies [8,31]. A linear increase of with 3 √ ch follows naturally if the hadronization volume is proportional to the number of produced particles.
For higher multiplicities, the measured parameter is observed to be independent of multiplicity. This saturation of the BEC parameter at high multiplicity is similar to that observed by ATLAS at √ = 7 TeV [28] and by CMS at √ = 13 TeV [38]. The dependence of for the region ch > 3.08 is then fitted simply by a constant: ( ch ) = .
For the combined MB and HMT data at √ = 13 TeV and T > 100 MeV the parameters of the ch dependence of and are given in the upper part of Table 2, with statistical and asymmetric systematic uncertainties summed in quadrature.
For T > 500 MeV the variation of the BEC parameters is fitted to the same functional forms as used for T > 100 MeV (Eqs. (5), (6), (7)), and the results are shown in the lower part of Table 2.
From Figure 4(a) it can be seen that ( ch ) shows only a weak dependence on multiplicity for T > 500 MeV. It is constant within measurement uncertainties, indicating a relatively high level of incoherence. Figures 4(b) and 4(c) show similar behaviour of ( ch ) for both T thresholds. There is a linear rise with 3 √ ch up to 1.2, after which there is some tapering off towards a plateau for 3 for T > 100 MeV.

Comparison of 2 ( ) functions at 13 and 7 TeV
The BEC parameters and were measured previously by ATLAS at 0.9 and 7 TeV [28]. The data sample at 7 TeV was large enough to explore high-multiplicity events, and a clear plateau in was observed for ch ≥ 50. This plateau at high multiplicity is confirmed by the present analysis at 13 TeV, as shown in Figures 4(b) and 4(c). However, a direct comparison of the BEC parameters obtained at the two centre-of-mass energies is complicated by small but significant differences between the Monte Carlo generators used to construct MC 2 ( ), and hence 2 ( ), in the two analyses (see Eqs. (1) and (2)). The BEC parameters obtained at 7 TeV are systematically lower than those at 13 TeV by about 1.5 This finding is supported by a comparison between the quantities data 2 ( ) at the two centre-of-mass energies. Figures 5(a) and 5(b) show data 2 ( ) for MB events at 7 and 13 TeV. Figure 5(a) shows a representative plot for a scaled-multiplicity interval (3.09 < ch ≤ 3.86) and Figure 5(b) shows a representative plot for a T interval (400 < T ≤ 500 MeV). In Figure 5 the uncertainties shown are the statistical and systematic uncertainties combined in quadrature. The systematic uncertainties include the effects of the track efficiency, Coulomb correction, non-closure and multiplicity unfolding. The plots shown in Figures 5(a) and 5(b) are representative of all multiplicity and T intervals. The data from the two centre-of-mass energies are in agreement everywhere within the uncertainties, and the ratio of 7 TeV to 13 TeV data is 1 to better than 1% except in the very first one or two bins.

Dependence on the pair average transverse momentum
The BEC are studied in T intervals from 0.1 to 1.5 GeV that are chosen so that they are populated roughly equally. The fitted dependences of ( T ) and ( T ) for the two T thresholds are shown in Figure 6 for T > 100 and 500 MeV. In the small-, BEC-sensitive region, the minimum T is highly correlated with the T threshold, and so T thresholds of 100 and 500 MeV are applied for the T thresholds of 100 and 500 MeV, respectively. The values of ( T ) and ( T ) decrease with increasing T . The decrease of ( T ) with T is described by the function: The ( T ) parameter shows similar behaviour: The fitted parameters in Eqs. (8) and (9) are given in Table 3.

Double-differential dependence on ch and T
There is sufficient data at 13 TeV to study the double-differential dependence: ( ch , T ) and ( ch , T ). Figure 7 shows the two-dimensional (2D) dependence of and on ch and T , for T > 100 MeV. In order to convey a clear picture of the 2D dependence, the values are presented without uncertainties, but they are generally less than 10%.
The following features can be seen in Figure 7: • The parameter decreases monotonically with T in all ch intervals, similar to the behaviour seen in Figure 6(a).
• The variation of with ch is rather flat for T 500 MeV, which may be compared with the clear decrease of as a function of ch seen in Figure 4(a) for T > 100 MeV. This is discussed below.   • The variation of with ch for T > 500 MeV indicates a slight rise with ch similar to what is seen in Figure 4(a) for T > 500 MeV.
• The parameter also decreases monotonically with T in all ch intervals, with the decrease being more pronounced in lower ch intervals. A similar decrease is seen in Figure 6(b).
• increases with ch in all T intervals, starting from a lower value at larger T .
• This increase of flattens off to a plateau level in all T intervals. This is similar to the plateau seen in Figures 4(b) and 4(c). Within the precision of the present data, the onset of the plateau occurs at  Figure 7: The two-dimensional dependence on ch and T for T > 100 MeV for (a) the correlation strength, , and (b) the source radius, , obtained from the exponential fit to the 2 ( ) correlation functions using the MB sample for ch ≤ 3.08 and the HMT sample for ch > 3.08. ch ∼ 3 in all T intervals.
These features can be seen more clearly in selected 1D projections in Figures 8 and 9 for and , respectively. MB and HMT data agree where they overlap.

Figures 8(c) and 8(d)
show the essentially flat behaviour of as a function of ch for low T , and the slight rise at higher T , which are described above. As also noted above, Figure 4(a) shows a slight but clear decrease of with ch for T > 100 MeV when integrated over all T . This arises because particle T and hence T of the pair increase with ch , but decreases as T increases. This leads to the overall decrease with ch seen in Figure 4(a) for T > 100 MeV.
The 2D behaviour can be parameterized according to Eqs. (8) and (9): and The parameters , , and are shown in Figures 10 and 11 for both the T > 100 MeV and T > 500 MeV selections.
The parameters , , and can themselves be fitted as functions of ch . For the parameters and T > 100 MeV:  where the fit for applies for ch > 3.0.
The 2D dependence has also been studied for the T > 500 MeV selection and is shown in Figure 12. As explained above, the T selection constrains the BEC parameters to the region T > 500 MeV.
The main features are similar to those observed in Figure 7. In particular, the value of tends to a plateau at large ch in all T intervals, but the plateau level is systematically lower than for T > 100 MeV, in line with what is observed in Figures 4(b) and 4(c). Selected 1D projections are shown in Figures 13 and 14 for and , respectively. As with the lower T selection, MB and HMT data agree where they overlap.
The 2D behaviour is parameterized according to Eqs. (10) and (11). The fit parameters , , and for T > 500 MeV are shown in Figures 10 and 11. As before, the parameters , , and are themselves  where the fit for applies for ch > 2.8.   Figure 12: The two-dimensional dependence on ch and T for T > 500 MeV for (a) the correlation strength, , and (b) the source radius, , obtained from the exponential fit to the 2 ( ) correlation functions using the MB sample for ch ≤ 3.08 and the HMT sample for ch > 3.08.

Comparison with CMS
As noted in the Introduction, BEC have been studied extensively at the LHC. In this section we compare our results on the BEC radius parameter with the only other published study at 13 TeV. The CMS Collaboration has studied BEC at 13 TeV in the kinematic region T > 200 MeV, | | < 2.4 and T < 1 GeV [38]. In the CMS publication comparison is made to the ATLAS 7 TeV results. To do so, CMS has adjusted its results for to match the ATLAS kinematic region, T > 100 MeV, | | < 2.5, by extrapolating the ch values to the lower T threshold and excluding the same regions as ATLAS in the fit to 2 ( ). These adjusted CMS values are shown together with the ATLAS results in Fig. 15. The behaviour of in the two experiments is qualitatively similar. However, there is a clear difference in the results, particularly at low ch . Results from this analysis ( Fig. 4(b)) and earlier studies for the lower-energy ATLAS analysis [28] indicate a direct dependence of on the T threshold. So the effect of the difference in kinematic coverage as well as the differing approaches to constructing the reference sample requires further investigation.

Summary and conclusions
Results are presented of measurements of two-particle BEC of like-charge hadrons with track T -thresholds of 100 and 500 MeV and | | < 2.5 produced in collisions at √ = 13 TeV recorded by the ATLAS detector at the CERN LHC. The study is carried out using data collected with the MB and HMT triggers. The integrated luminosities are 151 b −1 and 8.4 nb −1 for the MB and HMT data samples, respectively.
The studies are performed using the 2 ( ) correlation function, which is the ratio of the data 2 ( ) correlation function obtained from the data to the MC 2 ( ) calculated using MC events with no BEC effect. The 2 function is itself a ratio of like-charge particle pairs to a reference sample, and the reference sample is constructed using unlike-charge particle pairs. A clear indication of BEC is observed in the region of small four-momentum difference, 0.2 GeV. The BEC parameters, and , are studied as functions of a scaled charged-particle multiplicity ch , the pair average transverse momentum ( T ) and in ( ch , T )-intervals.
The parameter ( ch ) is found to increase as 3 √ ch for low multiplicity for both T thresholds. For T > 100 MeV, = 2.54 +0.12 −0.22 (total) fm, for scaled multiplicity up to ch ≈ 2. For ch ≥ 3, the source radius saturates at a value = 3.35 +0.20 −0.09 (total) fm, confirming the previous observation of a high-multiplicity plateau by ATLAS at √ = 7 TeV. For T > 500 MeV, the behaviour of is similar, but systematically lower. The parameter ( ch ) decreases with multiplicity for the lower T threshold and is lower for the higher T threshold but increases slightly with multiplicity. The parameters ( T ) and ( T ) both decrease with increasing pair average transverse momentum T . Considering the combined ( ch , T ) dependence, decreases with T in all ch intervals, for both T thresholds. The decrease is more pronounced at lower multiplicity. The radius parameter also decreases with T for T > 100 MeV, decreasing more strongly in the lower ch intervals. For T > 500 MeV the behaviour of is rather flat with T . As a function of ch , increases and reaches a plateau at large multiplicity in all ch intervals and for both T thresholds. Within the uncertainties of this analysis, the turn-on of the plateau occurs at the same value of ch in all T intervals and for both T thresholds. A comparison of the parameter as a function of ch with CMS results at 13 TeV shows qualitatively similar behaviour and the differences at low ch require further investigation.
The measurements presented here complement the ATLAS measurements at √ = 0.9 and 7 TeV by extending the studies to higher multiplicities and pair average transverse momenta, and provide a more detailed study of the BEC parameters by exploring the double-differential ( ch , T )-dependence. The present results confirm the saturation of the source radius at higher multiplicities with ch 100.

Appendix A Effect of resonances on 2 ( ) correlation function
The potential influence of resonances on the 2 ( ) is evident by studying the distributions of the UCP pairs (ℎ ± ℎ ∓ ) where both particles come from the same species of resonance; this is shown in Figure 16. The presented distributions are for track pairs with T > 100 MeV. The reflections of the resonance peaks can be seen clearly, and arise mainly from contributions from the , , and * 0 resonances, with much smaller contributions from the and mesons. The contributions of 0 (980) → + − and 2 (1270) → + − are insignificant.  The ATLAS Collaboration