Observation of Upsilon(1S) pair production in proton-proton collisions at sqrt(s) = 8 TeV

Pair production of Upsilon(1S) mesons is observed at the LHC in proton-proton collisions at sqrt(s) = 8 TeV by the CMS experiment in a data sample corresponding to an integrated luminosity of 20.7 inverse femtobarns. Both Upsilon(1S) candidates are fully reconstructed via their decays to mu+ mu-. The fiducial acceptance region is defined by an absolute Upsilon(1S) rapidity smaller than 2.0. The fiducial cross section for the production of Upsilon(1S) pairs, assuming that both mesons decay isotropically, is measured to be 68.8 +/- 12.7 (stat) +/- 7.4 (syst) +/- 2.8 (B) pb, where the third uncertainty comes from the uncertainty in the branching fraction of Upsilon(1S) decays to mu+ mu-. Assuming instead that the Upsilon(1S) mesons are produced with different polarizations leads to variations in the measured cross section in the range from -38% to +36%.


Introduction
The measurement of quarkonium pair production in proton-proton (pp) collisions provides insight into the underlying mechanism of particle production. The production mechanisms of heavy-quarkonium pairs such as J/ψ and Υ are especially valuable as they probe these interactions in both perturbative and nonperturbative regimes [1]. In this paper, the first cross section measurement of Υ(1S) pair production, reconstructed in the four-muon final state, in pp collisions from the CERN LHC at √ s = 8 TeV is reported.
Proton collisions can be described by parton models [2] in quantum chromodynamics (QCD). In this framework, each colliding hadron is characterized as a collection of free elementary constituents. Because of the composite nature of hadrons, in a single hadron-hadron collision two partons often undergo a single interaction (single-parton scattering, SPS). However, it is also possible that multiple distinct interactions (multiple-parton interactions, MPIs) occur, the simplest case being double-parton scattering (DPS). The SPS mechanism for heavy-quarkonium pair production can be described by a nonrelativistic effective theory [3] of QCD. However, contributions from the DPS mechanism are not addressed in a simple way within the framework of perturbative QCD [4], and DPS or MPIs are widely invoked to account for the observations that cannot be explained otherwise, such as the rates for multiple heavy-flavor production [5].
There are still large uncertainties owing to possible higher-order SPS contributions and the limited knowledge of the proton transverse profile [6]. Heavy-quarkonium final states are expected to probe the distribution of gluons in a proton since their production is dominated by gluon-gluon (gg) interactions [7,8]. Cross section measurements of quarkonium pair production are essential in understanding SPS and DPS contributions and the parton structure of the proton.
The first quarkonium pair production measurement dates back to 1982 when the NA3 experiment at CERN observed direct production of two J/ψ mesons [9] in π − -platinum interactions with beam energy of 150 GeV and 280 GeV, where the main contribution is quark-antiquark annihilation. Early theoretical calculations of the pair production cross section assumed only color-singlet states [10][11][12]. However, because of the high parton flux and high center-of-mass energy at the LHC, DPS is expected to play a significant role in quarkonium pair production [13].
Quarkonium pair production in pp collisions via DPS is assumed to result from two independent SPS occurrences [7,13]. Several DPS production processes, including final states with associated jets, are commonly described by an effective cross section (σ eff ) that characterizes the transverse area of the hard partonic interactions [14,15]. It is estimated as the product of single quarkonium production cross sections divided by the corresponding DPS quarkonium pair cross section [16]. Assuming the parton distribution functions are not correlated, σ eff is expected to be independent of the final state and the center-of-mass energy.

The CMS detector and muon reconstruction
The central feature of the CMS apparatus is a 13 m long superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Extensive forward calorimetry complements the coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The main subdetectors used for the present analysis are the muon detection system and the silicon tracker.
Charged particle trajectories are reconstructed by the silicon tracker in the pseudorapidity range |η| < 2.5. The innermost component of the tracker consists of three cylindrical layers of pixel detectors in the barrel region and two endcap disks at each end of the barrel. The strip tracker has 10 barrel layers and 12 endcap disks at each end of the barrel. There are 66 million 100 × 150 µm 2 silicon pixels and more than 9 million silicon strips. Strip pitches vary between 80 µm and 120 µm in the inner barrel, while the inner disk sensors have a pitch between 100 µm and 141 µm. For charged particles of transverse momentum 1 < p T < 10 GeV and |η| < 1.4, the relative track resolution is typically 1.5% in p T and 25−90 (45−50) µm in the transverse (longitudinal) impact parameter [25].
The outermost part of the CMS detector is the muon system, which covers |η| < 2.4. It consists of four layers of detection planes made of drift tubes, cathode strip chambers, and resistiveplate chambers. Muon reconstruction and identification begins within the muon system using hits in the muon chambers. A matching algorithm makes the best association between a track segment in the muon system and a track segment in the inner tracker. A track fit is performed combining all hits from both subsystems to select high-purity muon candidates.
Data are collected with a two-level trigger system. The first level of the trigger system, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events that pass the trigger requirements in a fixed time interval of about 4 µs, of which 1 µs is available for data processing. The second stage, the high-level trigger (HLT), is a processor farm that further reduces the event rate from around 100 kHz to less than 1 kHz before data storage. The HLT has full access to all the event information, including tracking, and therefore selections based on algorithms similar to those applied offline are used.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in ref. [26].

Data, simulation, and event selection
The data sample used in this analysis was collected with the CMS detector at the LHC in proton-proton collisions at a center-of-mass energy of 8 TeV and corresponds to an integrated luminosity of 20. . The CASCADE event simulation uses an off-shell matrix-element description of the g*g* → Υ(1S)Υ(1S) production and includes the k T -dependent parton distribution function JH-2013-set2 [29]. The parton evolution follows the CCFM prescription [30][31][32], which forms a bridge between the DGLAP and BFKL parton resummation models. The CMS detector response is simulated with the GEANT4 toolkit [33]. These samples are simulated with the appropriate conditions for the analysed data, including the effects of alignment, efficiency, and pileup.
For the rest of the paper, the notation ΥΥ will be used to denote any pair-wise combination of Υ(1S), Υ(2S), and Υ(3S) mesons. Candidate ΥΥ events are required to have four or more muons that pass the first-level muon trigger, of which at least three must be identified as muons by the HLT. There must be at least one pair of oppositely charged muons that have an invariant mass M µµ between 8.5 and 11 GeV, and a vertex fit chi-squared (χ 2 ) probability greater than 0.5%, as determined from a Kalman-filter algorithm [34].
The offline selection of ΥΥ candidates begins with a search for four muons with the sum of their charges equal to 0. Selected muons are further required to pass the following quality criteria: tracks of muon candidates must have at least six hits in the silicon tracker, at least one hit in the pixel detector, and they must match at least one muon segment in the muon detector. Loose selection cuts are applied on the longitudinal and transverse impact parameters of the muons to reject muons from cosmic rays and weakly decaying hadrons. To ensure a nearly uniform single-muon acceptance and a well-defined kinematic region, the muon candidates must have p µ T > 3.5 GeV and |η µ | < 2.4. To reconstruct the Υ(1S)Υ(1S) → 4µ decay, a kinematic fit [35] is performed on oppositely charged muon pairs and then on the four muons. The fit incorporates the decay and production kinematic properties to reconstruct the full decay chain of the user-defined process. To be compatible with the trigger requirements, each reconstructed Υ candidate is required to have a vertex fit χ 2 probability larger than 0.5%, as determined by kinematic fitting, and an invariant mass between 8.5 and 11 GeV. To suppress contributions from pileup events, all four muons are also required to be associated with a common vertex having a fit χ 2 probability larger than 5%.
In order to measure the cross section in a well-define kinematic region, both Υ candidates must have |y Υ | < 2.0, where y is the rapidity. Some of the events containing more than four muons result in multiple ΥΥ candidates per event since multiple dimuon pairs can pass the trigger and selection requirements. These events constitute about 4% of the selected events and are discarded from further analysis. After all selection criteria are applied, a total of 313 ΥΥ candidates are found. The major contribution of background comes from miscombined muons from Drell-Yan production and semileptonic b → µ + X or c → µ + X decays that pass the trigger and selection requirements.
To extract the signal yield of Υ(1S)Υ(1S) events, a two-dimensional (2D) unbinned extended maximum-likelihood fit to the invariant masses of the two reconstructed µ + µ − combinations is used. Two kinematic variables are employed to discriminate the Υ(1S)Υ(1S) signal from the background: the invariant mass of the higher-mass Υ candidate (M (1) µµ ) and the invariant mass of the lower-mass Υ candidate (M (2) µµ ). This choice of variables helps to resolve the ambiguity of having a dimuon pair in which an Υ(2S) candidate appears. There is no visible signal for Υ(3S), which is not considered further in this study. Figure 1 shows the distribution of M  The contributions from other categories (combinatorial-Υ(1S) and combinatorial-Υ(2S)), where two unassociated muons form the higher-mass Υ candidate, are found to be consistent with zero. There is no indication of Υ(2S)Υ(2S) events.
Each opposite-sign dimuon invariant mass distribution contains the contributions from Υ candidates and a nonresonant background. The Υ signal is modeled as a sum of two Crystal Ball functions [36] with a common mean, labeled DCB, to capture the variation in dimuon resolution as a function of y Υ . The DCB parameters are extracted from the signal SPS simulated sample and fixed in the fit to the data. The combinatorial background components are modeled with first-order Chebyshev polynomials (Poly). The distributions of the complete combinatorial background are found to be compatible with those of like-sign dimuon combinations. The yields of the five event categories and the two slope parameters for the combinatorial background are the free parameters in the fit to data. The likelihood function k for an event k (ignoring the normalisation term for simplicity) is given as: The yields N i for the five categories are determined using the RooFit package [37] by mini-  [38], and is found to be well above five standard deviations. The statistical significance of the Υ(2S)-Υ(1S) peak, evaluated using the same procedure, is 2.6 standard deviations. The invariant mass distributions of the higher-and lower-mass muon pairs in the data, with the result of the 2D fit superimposed, are shown in figure 2. The fit is validated by applying it to simulated samples generated using the probability distribution functions (pdf) that describe the data. No bias in the yields is found, and the likelihood value from the fit to the data closely agrees with the mean likelihood value obtained from the simulated samples.

Efficiency and acceptance
The muon acceptance is defined by the geometry of the CMS detector and the kinematic requirements of the analysis, which are established based on studies of simulated Υ(1S)Υ(1S) events.
The kinematic variable distributions of the Υ(1S)Υ(1S) system change with the underlying production mechanism. Therefore, to minimize the model dependence, the acceptance and efficiency corrections are calculated on an event-by-event basis using the measured Υ meson and muon momenta. This method was utilized in a similar CMS analysis [23], and is intended to minimize the model dependence of the correction factors.
The probability of an event passing the muon acceptance requirements is obtained by repeatedly generating the four-momenta of the four decay muons of the measured Υ pairs in the event. The acceptance correction for each selected event is calculated separately by simulating a large sample of such decays. The acceptance correction factor, a i for an event i, is defined as the number of times all four muons pass the acceptance criteria divided by the total number of trials. The angular distribution of the decay muons with respect to the direction of flight of the parent Υ(1S) meson, in the Υ(1S) meson rest frame, is assumed to be isotropic. Different polarization scenarios of Υ(1S) mesons are considered and discussed later.
The efficiency correction for the four-muon system is determined with a data-embedding method that repeatedly substitutes the measured muon four-momenta into different simulated events, which are then subjected to the complete CMS detector simulation and reconstruction chain. Each event is simulated with the same pileup scenario as in the corresponding data event. The efficiency correction factor i for the i th event in the data is the number of simulated events that pass the offline reconstruction and trigger requirements, divided by the total number of simulated events. For each data event, 1000 different simulated events are generated. The number of efficiency-corrected signal events, N true , is then given as ∑ i 1/ i , where the summation is over the data events that satisfy all the selection requirements.
Because of the finite µ + µ − mass resolution, typically near the kinematic acceptance of the detector, the measured momenta of reconstructed Υ candidates will be distorted from the original sample. To account for this, scaling factors are determined for both the SPS and DPS simulated samples. First, an average efficiency using a data-embedding method, 1 , is calculated as the number of events surviving the trigger and reconstruction requirements, divided by the number of efficiency-corrected events N true . Alternatively, an average efficiency using the generator information, 2 , is determined as the number of events satisfying the trigger and reconstruction criteria, divided by the number of events generated within the muon and Υ meson acceptance of the detector. Compared to the first method, the latter efficiency is based on the originally generated muon momenta. Then, the scaling factor is defined as the ratio 2 / 1 . Scaling factors of 98% and 96% are determined for the SPS and DPS models, respectively. The average of these scaling factors is applied to the data-embedding efficiency.

Systematic uncertainties
Several sources of systematic uncertainty in the cross section measurement are considered and discussed below.
To estimate the uncertainty caused by the imperfect modeling of the data, several variations are checked: (i) the parameters for the Υ(1S) resonance shape that are fixed for the maximumlikelihood fit are varied by their uncertainty, as determined from the simulated samples; (ii) the Υ(1S) resonance shape is alternatively parameterized with a relativistic Breit-Wigner convolved with a Gaussian function; (iii) Υ(1S) mesons are alternatively parameterized using signal parameters obtained from the DPS simulation samples. The largest difference in the signal yield among these fits of 7.9% is taken as the corresponding systematic uncertainty.
Evaluation of the trigger and reconstruction efficiencies for Υ pair production rely on detector simulation. To estimate the systematic uncertainty owing to the detector simulation, the fourmuon efficiency is also calculated from the single-muon efficiencies measured in data with a "tag-and-probe" method [39] that utilizes data samples containing a single J/ψ meson decaying into a muon pair. To calculate the four-muon efficiency using the tag-and-probe method, the correlations among the muons is separately derived from the simulation. A total correlation factor, determined from simulation as a function of the p T and rapidity of the four-muon system, is multiplied by the product of the four single-muon efficiencies, also measured versus η and p T . The difference of the corrected signal yield with the tag-and-probe method and the data-embedding method of 4.9% is taken as the systematic uncertainty from the simulation.
The efficiency correction method is evaluated using the full detector simulation and reconstruction event samples generated according to the two production mechanisms: SPS and DPS. The systematic uncertainty from using the efficiency-correction method, a data-embedding efficiency correction method is applied to obtain the per-event efficiency, i . The efficiencycorrected yield is compared to the number of generated events, and the fractional difference between the two is used to estimate the uncertainty in the correction method. This procedure is repeated for both simulated samples separately. The systematic uncertainty is then taken to be the larger of the differences in the two samples and is found to be 3.7%.
Signal SPS and DPS simulated samples are used to estimate the uncertainty in the event-byevent acceptance correction method. A sample of N gen simulated events are subjected to the acceptance criteria. The event-based acceptance correction is then applied to the events passing the full selection to arrive at a corrected yield, N gen . The uncertainty is taken as the relative deviation of N gen from N gen . The larger deviation of 2.8% between the SPS and DPS samples is taken as a systematic uncertainty.
The luminosity delivered to CMS is measured based on cluster counting from the pixel detector. The overall systematic uncertainty in the integrated luminosity is estimated to be 2.6% [40]. Table 1 provides a summary of the individual systematic uncertainties. The total systematic uncertainty in the Υ(1S)Υ(1S) cross section measurement of 10.7% is calculated as the sum in quadrature of the individual components. In this analysis, Υ(1S) mesons are assumed to decay isotropically. To check the sensitivity of the fiducial cross section measurement on the Υ(1S) meson decay angular distribution, studies are performed for several extreme polarization scenarios. With θ * defined as the angle between the µ + direction in the Υ(1S) meson rest frame and the direction of the Υ(1S) in the lab frame, the angular distribution is parameterized as 1 + λ θ cos 2 (θ * ), where λ θ is the polar anisotropic parameter [41]. The extreme cases λ θ = 0 or λ θ = +1 (−1) corresponds to an isotropic Υ(1S) meson decay or 100% longitudinal (transverse) polarization. The polar anisotropic parameter for each Υ(1S) candidate, labeled as λ θ1 and λ θ2 , are varied simultaneously, and Table 2 shows a summary of this study. Compared to an isotropic Υ(1S) meson decay, the corrected yield is 38% lower for λ θ1 = λ θ2 = −1 and 36% higher for λ θ1 = λ θ2 = +1. This variation is not included in the total systematic uncertainty, but is quoted separately.

Results
The fiducial cross section of Υ(1S) pair production, σ fid , is measured for events in which both Υ(1S) mesons have |y Υ | < 2.0. It is derived from the measured Υ(1S) pair yield, N signal , using the expression: where ω is the average correction factor that accounts for the inefficiencies stemming from reconstruction, identification, and trigger requirements, as well as the fraction of events lost because of the muon acceptance criteria, L is the integrated luminosity, and B(Υ(1S) → µ + µ − ) = (2.48 ± 0.05)% is the world-average branching fraction of the Υ(1S) to µ + µ − [42].
The factor ω is the inverse product of per-event efficiencies, i , and acceptances, a i , averaged over the events that pass selection cuts. For the calculation of the total cross section, values of L = 20.7 fb −1 , N signal = 38 ± 7, and ω = 23.06 are used. Under the assumption that both Υ(1S) mesons decay isotropically, the total cross section for Υ(1S) pair production is measured to be σ fid = 68.8 ± 12.7 (stat) ± 7.4 (syst) ± 2.8 (B) pb, where the uncertainties are statistical, systematic, and related to the Υ(1S) → µ + µ − branching fraction, respectively.
In quarkonium pair production, the measurement of the effective cross section σ eff depends on the fraction of DPS, which is usually estimated either as a residual to the SPS prediction or as the result of a fit to the rapidity or azimuthal angle difference between quarkonia pairs. To estimate σ eff using the fiducial Υ(1S) pair production cross section σ fid , the following formula is used [7]: where σ(Υ) is the cross section for single-Υ(1S) production [43] extrapolated to the fiducial region (|y Υ | < 2.0) and f DPS is the fraction of the DPS contribution. To estimate the effective cross section, we use σ(Υ) = 7.5 ± 0.6 nb and a value of f DPS ≈ 10% [8], which gives σ eff ≈ 6.6 mb. The cross section for Υ(1S) pair production, assuming SPS with feed-down from higher Υ states, is 48 pb [8]. This result, combined with our measurement, gives f DPS 30%, and with this estimation, σ eff is calculated to be ≈2.2 mb. These two estimates of σ eff are consistent with the range of values from the heavy-quarkonium measurements (2−8 mb) [21,22], but are smaller than that from multijet studies (12−20 mb) [17-20]. The quarkonium final states are dominantly produced from gluon-gluon interactions, while the jet-related channels that correspond to higher values of σ eff are produced by quark-antiquark and quark-gluon parton interactions. The trend in the measured effective cross section might indicate that the average transverse distance between gluons in the proton is smaller than between quarks, or between quarks and gluons.

Summary
The first observation of Υ(1S) pair production has been performed with the CMS detector in proton-proton collisions at √ s = 8 TeV from a sample corresponding to an integrated luminosity of 20.7 fb −1 . A signal yield of 38 ± 7 Υ(1S)Υ(1S) events is measured with a significance exceeding five standard deviations. Using the measured muon kinematic properties in the data, an event-based acceptance and efficiency-correction method is applied, minimizing the model dependence of the cross section measurement. The fiducial cross section of Υ(1S) pair production measured within the single-Υ acceptance |y Υ | < 2.0 is found to be σ fid = 68.8 ± 12.7 (stat) ± 7.4 (syst) ± 2.8 (B) pb, where the Υ(1S) decay into muons is assumed to be isotropic. Assuming instead that the Υ(1S) mesons are produced with different polarizations leads to variations in the measured cross section in the range from −38% to +36%. An effective cross section is also estimated using our result that is in agreement with the values from heavy-quarkonium measurements, but is smaller than that from multijet studies.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: