Measurement of the Y(1S) pair production cross section and search for resonances decaying to Y(1S)μ + μ − in proton-proton collisions at s=13TeV

The ﬁducial cross section for Υ ( 1S ) pair production in proton-proton collisions at a center-of-mass energy of 13 TeV in the region where both Υ ( 1S ) mesons have an absolute rapidity below 2.0 is measured to be 79 ± 11 (stat) ± 6 (syst) ± 3 ( B ) pb assuming the mesons are produced unpolarized. The last uncertainty corresponds to the uncertainty in the Υ ( 1S ) meson dimuon branching fraction. The measurement is performed in the ﬁnal state with four muons using proton-proton collision data collected in 2016 by the CMS experiment at the LHC, corresponding to an integrated luminosity of 35.9 fb − 1 . This process serves as a standard model reference in a search for narrow resonances decaying to Υ ( 1S ) µ + µ − in the same ﬁnal state. Such a resonance could indicate the existence of a tetraquark that is a bound state of two b quarks and two b antiquarks. The tetraquark search is performed for masses in the vicinity of four times the bottom quark mass, between 17.5 and 19.5 GeV, while a generic search for other resonances is performed for masses between 16.5 and 27 GeV. No signiﬁcant excess of events compatible with a narrow resonance is observed in the data. Limits on the production cross section times branching fraction to four muons via an intermediate Υ ( 1S ) resonance are set as a function of the resonance mass. quark. A generic search for narrow resonances with mass between 16.5 and 27 GeV and decaying to a Υ ( 1S ) meson and a pair of muons is also presented. The ﬁnal state is the same as for the measurement of the Υ ( 1S ) pair production cross section, and a similar event selection is used. The Υ ( 1S ) pair production is a background to the resonance search.


Introduction
Quarkonium pair production is an important probe of both perturbative and nonperturbative processes in quantum chromodynamics. Experimental studies of this process can provide valuable information about the underlying mechanisms of particle production and improve our understanding of numerous physics processes that are treated as backgrounds in searches and measurements. Quarkonium pairs may originate from single-parton scattering (SPS) or double-parton scattering (DPS). These production mechanisms can be separated experimentally since the DPS production is characterized, among other features, by more forward and separated mesons. The analysis of nonperturbative effects is easier for quarkonium states composed of b quarks, as their large masses allow them to be approximated as nonrelativistic systems [1]. The CMS Collaboration observed for the first time the production of a pair of Υ(1S) mesons, using proton-proton data collected at a center-of-mass energy of 8 TeV [2]. This Letter presents a measurement of the Υ(1S) pair production cross section at a center-of-mass energy of 13 TeV. The cross section is measured in the fiducial region where both Υ(1S) mesons have an absolute rapidity below 2.0, using the final state with four muons. Additionally, the DPS contribution to the process is measured for the first time.
The Υ(1S) pair production can serve as a reference in searches for tetraquarks or generic resonances with masses close to twice the Υ(1S) meson mass. A light resonance decaying to a Υ(1S) meson and a pair of leptons might be the signature of a tetraquark characterized as a bound state of two b quarks and two b antiquarks, especially if its mass is below twice the η b mass [3][4][5][6][7][8][9][10][11][12][13]. In this Letter, in addition to the measurement of the Υ(1S) pair production cross section, we describe a search for tetraquarks with masses between 17.5 and 19.5 GeV, since bbbb tetraquarks would be expected to have a mass around four times that of the bottom quark. A generic search for narrow resonances with mass between 16.5 and 27 GeV and decaying to a Υ(1S) meson and a pair of muons is also presented. The final state is the same as for the measurement of the Υ(1S) pair production cross section, and a similar event selection is used. The Υ(1S) pair production is a background to the resonance search.
The LHCb Collaboration searched for bbbb tetraquarks using data collected at center-of-mass energies of 7, 8, and 13 TeV, without finding any hint of a signal [14]. This analysis probes a kinematic region that is not accessible with the LHCb detector and extends the covered mass range in the context of the generic search.
The Υ(1S) pair production fiducial cross section measurement and the resonance search are based on proton-proton collision data collected in 2016 at a center-of-mass energy of 13 TeV by the CMS experiment at the CERN LHC, corresponding to an integrated luminosity of 35.9 fb −1 .

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume, there 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. Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system [15]. 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. [16].
Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. Matching muons to tracks measured in the silicon tracker results in a relative p T resolution in the range 0.8-3.0% for muons with p T less than 10 GeV [17].

Simulated samples
The Υ(1S) pair production signal is simulated using the PYTHIA 8.226 generator [18], separately for the SPS and DPS mechanisms, under the assumption that the mesons are produced unpolarized. The DPS sample is produced by generating two hard interactions with color-singlet production of bottomonium states via gg → bb or color-octet production of bottomonium states via qq → bb. The invariant mass distribution of the meson pair and of the rapidity separation between the mesons are used to extract the fraction of DPS production, as detailed in Section 5. For this measurement, the distributions of these variables for the SPS process are taken from the next-to-leading-order (NLO*) calculation with a cutoff color-singlet mechanism (CSM) [19][20][21] using HELAC-ONIA 2.0.1 [22,23].
The signal of a narrow resonance decaying to a Υ(1S) meson and a pair of muons is modeled using different physics assumptions depending on the nature of the resonance: • a bottomonium state with the properties of the χ b1 (1P), assuming a phase-space decay to a Υ(1S) meson and a pair of muons, using the PYTHIA 8.226 generator; • a scalar particle produced in gluon fusion, using the JHUGEN generator [24][25][26][27]; • a pseudoscalar particle produced in gluon fusion, using the JHUGEN generator; • a spin-2 particle produced in gluon fusion, using the JHUGEN generator.
The signals are generated assuming the narrow-width approximation. The χ b1 (1P) sample is used to model the tetraquark signal, for which no dedicated generator exists. The other samples correspond to the signals in the generic search over an extended mass range. For each model, four resonance mass values are simulated: 14, 18, 22, and 26 GeV. Since the signal acceptance falls steeply around and below 14 GeV in the simulated samples, the probed mass range in this analysis is restricted to stay well above this mass threshold. The different mass points are used to interpolate and extrapolate the signal model over the whole mass range.
The PYTHIA generator with the tune CUETP8M1 [28] is used to model the parton shower and hadronization processes. Generated events are processed through a simulation of the CMS detector based on GEANT4 [29].
described in Ref. [17]. About 25% of simulated signal events and about 30% of data events have more than four such muons. Possible combinations of four muon tracks are refit with a constraint to come from a common vertex, and the χ 2 probability of the fit is determined. The combination of four muons with the largest χ 2 probability is chosen. For simulated signal events with more than four reconstructed muons, the correct muons are chosen in about 98% of cases. Among the four muons, at least three need to be associated with the trigger-level objects. At least two muons must be associated with the objects that passed the Υ mass compatibility and vertex criteria of the trigger, and they are paired together. If there are more than two such muons, which happens for 2 to 35% of simulated signal events depending on the resonance mass, those that have opposite-sign (OS) charges and an invariant mass closest to the worldaverage Υ(1S) mass [31] are paired together.
After selecting the best combination of four muons with p T > 2 GeV, the p T threshold is raised to 2.5 GeV for the selected muons. The final selection requiring p T > 2.5 GeV reduces the background from misidentified muons by about a factor of two. The muons are required to satisfy the medium muon identification criteria described in Ref. [17]. Both pairs of muons have to be composed of OS muons. The vertex fit χ 2 probability of the four muons is required to be greater than 5%, whereas that of the Υ(1S) candidate is required to be above 0.5%, similar to the requirement already imposed at trigger level. The muons are required to be separated from each other by at least ∆R = √ (∆η) 2 + (∆φ) 2 = 0.02, where ∆η and ∆φ are the differences in pseudorapidity and azimuthal angle between the muons. The positively (negatively) charged muon from one of the pairs can be paired with the negatively (positively) charged muon of the other pair to form so-called alternative pairs of OS muons. If one of these alternative pairs has an invariant mass compatible with a J/ψ particle within two standard deviations of the experimental resolution, which ranges between about 0.03 and 0.12 GeV depending on the muon pair kinematics, the event is discarded from the analysis. Events are also discarded if they contain two OS pairs of muons with invariant mass less than 4 GeV.
The selection criteria detailed above are common for the measurement of the Υ(1S) pair production cross section and the search for a resonant signal. The criteria that differ between the measurement and the search are described in the following. In the measurement of the Υ(1S) pair fiducial cross section, the reconstructed absolute rapidity of both muon pairs is required to be less than 2.0. In addition, for muons with |η| < 0.9, the p T threshold is raised to 3.5 GeV. Central muons with transverse momentum below 3.5 GeV have a high probability of being absorbed in the calorimeter or undergoing significant multiple scattering before reaching the muon detectors. This selection criterion reduces the systematic uncertainty in the muon reconstruction related to the detector simulation. It is, however, not used in the resonant search because it would strongly reduce the signal acceptance for the lower-mass signal range. In the resonance search, the invariant mass of the Υ(1S) candidate is required to be within two standard deviations of the experimental resolution from the Υ(1S) mass [31], where the resolution varies between about 0.06 and 0.15 GeV depending on the event.
The mass range of interest is known a priori for the search of a bbbb tetraquark signal. In this case, all the selection criteria described above have been determined and fixed in a blinded way, using simulation and without looking at data events with four muons having an invariant mass between 17.5 and 19.5 GeV.

Measurement of the Υ(1S) pair production cross section
The methodology used to measure the Υ(1S) pair production cross section is detailed in Section 5.1. After discussing the systematic uncertainties in Section 5.2, the results of the measurement of the inclusive Υ(1S) pair production fiducial cross section are presented in Section 5.3. Nonisotropic decays of the Υ(1S) mesons would change the measured cross section. Section 5.4 describes how the cross section would vary for nonzero values of the polarization parameters. Finally, the DPS and SPS mechanisms can be separated experimentally by measuring the Υ(1S) pair production cross section in bins of the rapidity difference between the mesons, |∆y(Υ(1S), Υ(1S))|, and of the invariant mass of the meson pairs, m Υ(1S)Υ(1S) . A measurement of the DPS-to-inclusive cross section ratio in the fiducial region is presented in Section 5.5.

Methodology
The Υ(1S) pair production cross section is measured in the fiducial region where both mesons have an absolute rapidity below 2.0. No other requirement is applied to define the fiducial region. The fiducial cross section, σ fid , can be expressed as: where N corr is the number of signal events corrected for the acceptance and efficiency of the selection, L is the integrated luminosity, and B stands for B(Υ(1S) → µ + µ − ) = (2.48 ± 0.05)% [31], which is the branching fraction of the Υ(1S) meson decay to a pair of muons.
To extract N corr from the data, we perform an extended unbinned two-dimensional (2D) maximum likelihood fit of the invariant mass distributions of two OS muon pairs, where all events are weighted for the acceptance and efficiency on an event-by-event basis by the weight ω, defined as: where the different terms are described below: • A, the probability for a Υ(1S) meson with an absolute rapidity below 2.0 and decaying to a pair of muons to have two muons in the geometrical acceptance of the detector (muon |η| < 2.4); No strong correlation between the acceptance values of the two mesons are found with a closure test described in Section 5.2, and the total acceptance is therefore computed as the product of the per-meson weights; • reco , the probability for a Υ(1S) meson with an absolute rapidity below 2.0 and decaying to a pair of muons each with |η| < 2.4 to have two reconstructed muons passing the identification and kinematic criteria listed in Section 4; • vtx , the probability for a Υ(1S) meson passing the acceptance reconstruction criteria outlined in items 2 and 3 to have a vertex fit χ 2 probability above 0.5%; • evt , the probability for an event where both Υ(1S) candidates pass all the criteria of items 2 and 3, and at least one of them passes the vertex fit χ 2 probability criterion of item 4, to pass the following event-level criteria: the trigger requirements, the fourmuon vertex fit χ 2 probability above 5%, and the absence of OS dimuon pairs with an invariant mass within two standard deviations of the world-average J/ψ meson mass [31].
The first three items in the above list are calculated as a function of the Υ(1S) rapidity and p T . The values of A, reco , and vtx , range between 0.47 and 1.00, 0.23 and 0.88, and 0.81 and 0.98, respectively, depending on the Υ(1S) rapidity and p T . The factor evt is calculated as a function of the p T of both Υ(1S) candidates, and ranges between 0.33 and 0.65. The subscript indices in Eq.
(2) indicate the Υ(1S) candidate to which the weight corresponds. The factor vtx enters the formula differently from the other acceptance and efficiency terms because the dimuon vertex fit χ 2 probability criterion needs to be satisfied by at least one of the two Υ(1S) candidates, but not necessarily by both. The weight ω is computed on an event-by-event basis, using the kinematic quantities of the reconstructed Υ(1S) candidates in data. They are estimated from simulation as efficiency maps and are similar for the SPS and DPS production modes, despite different correlations between the mesons. Data-to-simulation corrections for the trigger and muon identification efficiencies are taken into account in the computation of N corr .
In about 3% of cases, the four reconstructed muons are not correctly paired in the SPS and DPS Υ(1S) pair simulations. These events cannot be identified as part of the signal by the 2D fit since their distribution is similar to that of the floating combinatorial background. Therefore, the value N corr extracted from the fit is corrected by +3% to take into account these mispairings.
In the 2D fit, the muons are paired as described in Section 4, and the invariant masses of the two pairs are randomly denoted m 12 and m 34 . The signal model corresponds to Υ(1S) + Υ(1S) events, whereas the background model is the sum of the following physics processes: The shape of the invariant mass distribution for the Υ(1S) component is determined from a 2D fit of the two dimuon invariant masses in the Υ(1S) pair SPS simulation. The results are verified to be compatible with those of a fit performed using the simulated DPS events, even if the muon rapidity distributions differ between production modes. The m 12 and m 34 distributions are fitted with the sum of two same-mean Crystal Ball functions, which correspond to a power law tail added to a Gaussian core. This allows the radiative tails of the distributions to be well modeled. Figure 1 shows the projection of the 2D fit on the m 12 axis for Υ(1S)Υ(1S) simulated events. The projection on the m 34 axis is statistically identical and therefore not shown. The fitted mean of the Crystal Ball functions in simulation is compatible within one standard deviation with the world-average mass of the Υ(1S) meson, while the full width at half maximum is about 0.19 GeV, which is several orders of magnitude larger than the world-average width of the Υ(1S) meson [31] because of the limited detector resolution.
The contributions from Υ(2S) and Υ(3S) mesons are small, and the dimuon invariant mass distributions for these mesons are taken from a control region in data with events with two muons and two additional tracks that do not correspond to muon candidates. Both processes are modeled with a Gaussian function.
The combinatorial background components in the m 12 and m 34 distributions are modeled with second-order Chebychev polynomials with identical parameters. The number of degrees of freedom has been determined with a Fisher F-test [32], where the distribution of the combinatorial background is found by inverting the muon pair association in the signal region. The parameters of the polynomial are free to float in the 2D fit to data in the signal region, detailed in Section 5.3.
In the 2D fit to the data performed in the signal region, the free parameters are the normalizations of all the processes and the parameters of the combinatorial background mass distribution. The function parameters of the Υ(1S), Υ(2S), and Υ(3S) signal shapes are constrained within their uncertainties.

Systematic uncertainties
The normalization uncertainties that affect the measurement are the following: • 2.5% uncertainty in the integrated luminosity for the 2016 running period [33], which appears in Eq. (1).
• 0.5% uncertainty per muon in the efficiency of the muon identification and tracking, measured with a tag-and-probe method [17]. It sums up to 2% per event because the uncertainties are assumed to be correlated for the four muons since they mostly originate from the same source. This uncertainty is related to the term reco in the weight ω.
• 1% uncertainty in the vertex fit χ 2 probability criterion, determined by comparing background-subtracted observed and simulated distributions of the vertex fit χ 2 probability for events with a Υ(1S) meson and two nearby tracks. This uncertainty is related to the term vtx in the weight ω.
• 2% uncertainty per muon matched to trigger objects in the trigger efficiency, measured with a tag-and probe method, summing up to 6% per event because the uncertainties are assumed to be correlated for the three muons required at trigger level. This uncertainty is related to the term evt in the weight ω.
These normalization uncertainties propagate directly into identical uncertainties in the Υ(1S) pair production cross section. Additionally, the uncertainty of 2% in the B(Υ(1S) → µ + µ − ) branching fraction, which is used to compute N corr based on Eq. (1), results in a 4% uncertainty in the Υ(1S) pair production cross section measurement.
The parameters of the combinatorial background are freely floating, while the parameters of  The consistency of the method to obtain N corr is checked by applying the efficiency and acceptance weights to the events selected in simulation, and comparing the computed N corr to the number of events generated in the fiducial region before applying any selection criterion. This test is performed for both the SPS and DPS simulations using the correction maps derived from one sample, the other one, or their combination. Using the combined map, the weighted DPS yield has a deviation of (−1.3 ± 3.7)% with respect to the generated yield, and the corresponding deviation for the SPS sample is (−0.6 ± 1.5)%. The level of closure is similarly good for both production modes despite average event weights differing by more than a factor of 3 because of the kinematic differences. The weighted number of data events used to compute the Υ(1S) pair production cross section is increased by 1% to allow for a potential nonclosure, and an uncertainty of 1.5% is associated with this correction.
The systematic uncertainties are summarized in Table 1.

Measurement of the fiducial cross section
The 2D unbinned fit to the m 12 vs. m 34 distribution yields N corr = 1740 ± 240 for the Υ(1S)Υ(1S) process. The projections on both dimensions with all the fit components are shown in Fig. 2. This number of events can be translated into an inclusive cross section for the Υ(1S)Υ(1S) process in the fiducial region defined such that both Υ(1S) mesons have an absolute rapidity below 2.0. Taking into account the statistical and systematic uncertainties described in Section 5.2, and assuming unpolarized Υ(1S) mesons, the inclusive fiducial cross section is measured to be: where the last uncertainty comes from the uncertainty in the Υ(1S) dimuon branching fraction.
The CMS Collaboration previously measured, in the same fiducial region, the Υ(1S)Υ(1S) production cross section at a center-of-mass energy of 8 TeV to be 69 ± 13 (stat) ± 7 (syst) ± 3(B) pb [2]. Assuming all uncertainties are uncorrelated with those in the result presented in this Letter except that in the branching fraction of the Υ(1S) meson to muons, the measured ratio of the cross section at a center-of-mass of 13 TeV to that at 8 TeV is 1.14 ± 0.32, where the uncertainty includes both the statistical and systematic components. The PYTHIA generator   predicts a ratio of 2.1 for DPS production, and 1.6 for the SPS production. Taking the fraction of the DPS mechanism in the total cross section f DPS = (39 ± 14)% at a center-of-mass energy of 13 TeV, as measured in Section 5.5, the cross section ratio predicted by PYTHIA is 1.79 ± 0.27. Combining the uncertainties in quadrature, the prediction is within two standard deviations of the measurement.
Another unbinned extended maximum likelihood fit is performed to extract the number of Υ(1S)Υ(1S) events observed in data after the selection. The Υ(1S)Υ(1S) unweighted signal yield is obtained from a fit where all observed events have a weight of 1.0. For this fit, a separate signal shape is determined by fitting the m 12 and m 34 distributions in the unweighted simulation. The absence of weighting does not significantly modify the signal distribution. The unweighted event yields are given for all processes in Table 2. There is no evidence for the simultaneous production of two excited states of the Υ meson, but excesses with a significance lower than two standard deviations indicate the possible presence of Υ(1S)Υ(2S) and Υ(1S)Υ(3S) events. The number of events from data in the m 12 vs. m 34 distribution is shown in Fig. 3, along with the results of the fit to the signal+background model, using the color scale to the right of the plot.

Effect of the polarization
The acceptance and efficiency corrections have been computed assuming negligible polarization of the Υ(1S) mesons. A different assumption on the polarization can change the measured fiducial cross section. The polarization of the Υ(1S) states affects the angular distributions of the leptons produced in the Υ(1S) → µ + µ − decays through the following formula [34]: (1 + λ θ cos 2 θ + λ φ sin 2 θ cos 2φ + λ θφ sin 2θ cos φ), where θ and φ are the polar and azimuthal angles, respectively, of the positively charged muon with respect to the z axis of a polarization frame, and λ θ , λ φ , and λ θφ are the angular distribution parameters. To estimate the effect of the polarization on the measurement of the Υ(1S)Υ(1S) fiducial cross section, we choose to use the helicity frame, where the polar axis coincides with the direction of the Υ(1S) momentum. Measurements performed by the CMS and LHCb Collaborations on single Υ production indicate compatibility of all the angular distribution parameters with zero over a large phase space [35,36]. However, the same may not be true for Υ(1S) pair production. To estimate the effect of polarization on the Υ(1S) pair production cross section, simulated events are reweighted to have the angular distributions corresponding to various λ θ values, without changing the overall simulated yield. The same efficiency and acceptance corrections as in Eq. (2) are used to calculate N corr for these different polarization scenarios. The variations in the measured Υ(1S) pair production cross section are given for different λ θ coefficients in Table 3. The effect of different polarizations can be substantial, changing the measured cross section by −60 to +25%.

Measurement of the DPS-to-inclusive fraction
The DPS and SPS mechanisms lead to different kinematic distributions for the Υ(1S)Υ(1S) events. The DPS production is characterized by a larger separation in rapidity between the mesons, |∆y(Υ(1S), Υ(1S))|, as they are largely uncorrelated, and by a larger invariant mass of the meson pairs, m Υ(1S)Υ(1S) . The distributions of ∆φ(Υ(1S), Υ(1S)), ∆R(Υ(1S), Υ(1S)), and p T (Υ(1S)Υ(1S)) also differ for the SPS and DPS mechanisms, but they are very sensitive to the choice of model parameters in the simulation and are subject to large theoretical uncertainties [37]. Measuring the Υ(1S)Υ(1S) fiducial cross section in bins of |∆y(Υ(1S), Υ(1S))| or of m Υ(1S)Υ(1S) can give a measurement of the fraction of DPS events, f DPS , defined as: where σ DPS fid and σ SPS fid are, respectively, the DPS and SPS cross sections in the fiducial region. We measure the fiducial cross section in five bins of |∆y(Υ(1S), Υ(1S))| and five bins of m Υ(1S)Υ(1S) . The signal and background models are the same as for the inclusive measurement, except that the width of the function describing the Υ(1S) invariant mass shape is allowed to float between its best-fit values for the inclusive selection and for the selection in the relevant exclusive bin. This allows for a potential degradation (improvement) of the muon momentum resolution at high (low) pseudorapidity to be taken into account, since the muon pseudorapidity is correlated with both |∆y(Υ(1S), Υ(1S))| and m Υ(1S)Υ(1S) . The systematic uncertainties are identical to those presented in Section 5.2.
The extracted fiducial cross sections as a function of |∆y(Υ(1S), Υ(1S))| and m Υ(1S)Υ(1S) are compared to the expected distributions for SPS and DPS production, as obtained in the fiducial region using PYTHIA for the DPS process, and from HELAC-ONIA with the NLO* CSM predictions for the SPS process. The fraction f DPS is measured with a binned maximum-likelihood fit of these two simulated distributions with floating normalizations to the measured fiducial cross sections in bins of |∆y(Υ(1S), Υ(1S))| and m Υ(1S)Υ(1S) . As determined from pseudoexperiments, the best precision is expected to be achieved using |∆y(Υ(1S), Υ(1S))|. Theoretical uncertainties coming from the choice of parton distribution functions and the factorization and renormalization scales are taken into account for both the SPS and DPS predicted distributions. The fraction f DPS is measured to be (39 ± 14)% using |∆y(Υ(1S), Υ(1S))| as the discriminative distribution. This results includes both statistical and systematic uncertainties, where the former strongly dominates. The result using m Υ(1S)Υ(1S) is compatible with this measurement, but with much lower precision: (27 ± 22)%. The uncertainties are strongly dominated by the uncertainties in the measurements of the cross section in the |∆y(Υ(1S), Υ(1S))| and m Υ(1S)Υ(1S) bins, with theoretical uncertainties in the predicted SPS and DPS distributions playing a role at the percent level. The measured differential fiducial cross sections are shown in Fig. 4, together with the SPS and DPS predictions.
6 Search for resonances 6

.1 Methodology
We search for a narrow excess of events above an expected smooth four-muon invariant mass spectrum. Assuming that the resonant state decays into two muons and a Υ(1S) meson that further decays to a pair of muons, the signal mass resolution can be improved by using a massdifference observable [38]: where m 4µ is the invariant mass of the four leptons, m µµ the invariant mass associated with the Υ(1S) candidate, and m Υ(1S) the nominal mass of the Υ(1S) particle (9.46 GeV [31]). This estimated mass, denoted as m 4µ , has a resolution about 50% better than the four-muon invariant The results are extracted by performing an unbinned maximum-likelihood fit to the m 4µ spectrum. The signal and background components are modeled by several functional forms in the fit, as described in the next paragraphs.
The signal distributions are parameterized by the sum of two Gaussian functions with the same mean. The parameters are extracted for the four mass points available in simulation. The signal modeling needs to be interpolated for masses between 16.5 and 26 GeV and extrapolated to masses up to 27 GeV to search for narrow resonances with any mass between 16.5 and 27 GeV. This is done by fitting with polynomials the different parameters of the two Gaussian functions as a function of the generated resonance mass. The same procedure is repeated for every signal model. The full width at half maximum is about 0.2 GeV for a resonance mass of 18 GeV.
The background is separated into two components: the Υ(1S)Υ(1S) process, which was the signal in Section 5 and is characterized by a sharp rising edge in the m 4µ spectrum at twice the Υ(1S) meson mass, and the combinatorial background, which is described by a smooth function as explained below.
The m 4µ spectrum for the Υ(1S)Υ(1S) process is obtained from simulation, and is modeled as the product of a sigmoid function and an exponential function with a negative exponent. The nominal model for the Υ(1S)Υ(1S) background is taken as an average between the DPS and SPS templates, which is consistent with the measurement of the DPS fraction presented in Section 5.3. Figure 5 shows the m 4µ models obtained from simulated DPS and SPS events, together with the average fit model. The number of Υ(1S)Υ(1S) events in the signal region is extracted, as detailed in Section 5, using the selection designed for the resonance search and without applying the acceptance and efficiency corrections from Eq. (2). In this case, only events with 13 < m 4µ < 28 GeV are retained and no rapidity criteria are applied for the reconstructed Υ(1S) candidates. The yield is measured to be 78 ± 13 events. The requirement that the mass of a dimuon pair is compatible with the mass of a Υ(1S) meson within two standard deviations is enforced in the resonance search but is not applied to extract the yield because the 2D fit relies on the mass tails to estimate the combinatorial background. Since the efficiency of this criterion is 95% in both the SPS and DPS Υ(1S)Υ(1S) simulations, the Υ(1S)Υ(1S) yield in the signal region is corrected to 74 ± 13. The normalization of the Υ(1S) pair production process and its uncertainty are extracted from the same data as in the signal region of the resonance search, but this does not lead to a significant overconstraint of the uncertainty in the maximum-likelihood fit of the m 4µ distribution because the latter can determine the Υ(1S) pair normalization only with poor precision. The m 4µ spectrum for the combinatorial background is obtained in the fit to the data in the signal region. Several generic functions are used to parameterize this smooth background: • Chebychev polynomials of various orders; • the sum of a Gaussian function and a Chebychev polynomial; • the sum of a Breit-Wigner function and a Chebychev polynomial.
The widths of the Gaussian and Breit-Wigner functions are constrained to be above 2 GeV to avoid fitting narrow structures due to statistical fluctuations. We verify, using a control region where the vertex fit χ 2 probability of the four muons is in the range 10 −10 -10 −3 , that these three functional forms describe the smooth m 4µ spectrum of the combinatorial background with a good χ 2 probability. Muons with a vertex probability in this range are likely to be associated with processes from the same primary vertex, but can originate from decays in flight or displaced secondary vertices. This control region is shown in Fig. 6 for illustrative purposes. The parameters of the functions determined from the fit are not used in the signal region, where the parameters of the combinatorial background, as well as the choice of the functional form, are freely floating.

Systematic uncertainties
The systematic uncertainties are to a large extent similar to those used in the measurement of the Υ(1S)Υ(1S) cross section and introduced in Section 5.2. In this section, only the differences are highlighted. They arise from slightly different selection criteria, a different choice of ob- The uncertainty per muon in the muon identification and tracking is increased from 0.5% to 1% because poorly reconstructed muons with p T < 3.5 GeV in the barrel are included in the resonance search to increase the signal acceptance for light resonances. In addition, in the resonance search, the signal is affected by a 1% yield uncertainty related to the requirement that the Υ(1S) candidate has an invariant mass compatible with the nominal Υ(1S) meson mass within two standard deviations. This uncertainty is determined by comparing the dimuon invariant mass resolution distributions in Υ(1S)Υ(1S) simulated events and in Υ(1S) events in data. The modeling of the signal process with a resonance mass other than those for which simulated samples were generated leads to a 2% uncertainty in the signal normalization for the resonance search.
The discrete profiling method [39] is used to model the combinatorial background. This allows the choice of the fit functions among those provided to be considered as a discrete nuisance parameter. The parameters of these fit functions are freely floating.
The normalization of the Υ(1S)Υ(1S) background in the resonance search is extracted from the 2D unbinned fit to the invariant mass of the dimuon pairs in the Υ(1S) mass region. The uncertainty in the yield obtained from the fit is considered as a log-normal uncertainty in the fit to the m 4µ distribution. The m 4µ distribution of the Υ(1S)Υ(1S) background is allowed to float between the predictions for the SPS and DPS simulations.
Uncertainties in the m 4µ distribution of the resonant signal take into account the limited size of the simulated samples, and the limited precision of the description of the signal for masses not available in simulations. The uncertainty in the mean mass of the signal is 0.2%, corresponding to the uncertainty in the muon momentum scale. The other parameters describing the shape of the signal have an uncertainty between 5 and 15%, which leads to a combined impact on the final upper limits of less than 2%.
The uncertainty in the Υ(1S) dimuon branching fraction is not considered, since the limits are set on the product of the resonance production cross section and its branching fraction to four muons via an intermediate Υ(1S) resonance.

Results
The binned m 4µ distribution in the signal region of the resonance search is shown in Fig. 7. The background and example signal components are shown using their best-fit shapes and normalizations. Using the number of Υ(1S)Υ(1S) events observed in data as a reference, a resonance with a mass around 19 GeV and having a similar production cross section and branching fraction to four muons as the Υ(1S)Υ(1S) production, would produce about 100 events in our sample, given the similarity between the kinematic distributions of both processes. No significant narrow excess of events is observed above the background expectation. The largest excess is observed for a resonance mass of 25.1 GeV, and has a local significance of 2.4 standard deviations for the scalar signal hypothesis.  Upper limits on the product of the production cross section of a resonance and the branching fraction to a final state of four muons via an intermediate Υ(1S) resonance are set at 95% confidence level (CL) using the modified frequentist construction CL s in the asymptotic approximation [40][41][42][43][44], separately for each signal model. The upper limits are extracted using unbinned distributions. The cross section is defined in the entire phase space without fiducial requirements, and the branching fraction used is the product of the branching fraction of the resonant state to a Υ(1S) meson and two muons, and the branching fraction of the Υ(1S) meson to two muons. Masses between 17.5 and 19.5 GeV are probed in the context of the tetraquark search, using the bottomonium model, whereas the limits in the extended mass range 16.5-27 GeV are set for the generic search, using the JHUGEN models. The corresponding upper limits are given in Fig. 8. They range between 5 and 380 fb, depending on the mass and signal model. The patterns in the limits are broader for the spin-2 signal than for the scalar and pseudoscalar states because the signal is characterized by softer and more forward muons, leading to a worse m 4µ resolution.

Summary
The cross section for Υ(1S) pair production is measured in the fiducial region where both Υ(1S) mesons have an absolute rapidity below 2.0. The measurement is performed using protonproton collision data collected at a center-of-mass energy of 13 TeV by the CMS detector in 2016  and corresponding to an integrated luminosity of 35.9 fb −1 . Assuming that the Υ(1S) mesons are produced unpolarized, the fiducial Υ(1S) pair production cross section is determined to be 79 ± 11 (stat) ± 6 (syst) ± 3(B) pb, where the last uncertainty comes from the uncertainty in the Υ(1S) dimuon branching fraction. The result can change if the Υ(1S) mesons are produced with a nonzero polarization. Changing the polarization coefficient λ θ from −1 to +1, the resulting Υ(1S) pair production cross section measurement varies by −60 to +25%.
The contribution of double-parton scattering to the total inclusive Υ(1S) pair production cross section is determined for the first time. It is measured to be (39 ± 14)% in the same fiducial region as described above, where the uncertainty includes both statistical and systematic components, with the statistical uncertainty dominating.
The results of a search are also presented for a light narrow resonance, such as a tetraquark or a bound state beyond-the-standard model, decaying to a Υ(1S) and a pair of opposite-sign muons. No excess of events compatible with a signal is observed in the four-muon invariant mass spectrum. Upper limits at 95% confidence level on the product of the signal cross section and branching fraction to four muons via an intermediate Υ(1S) resonance are set for different signal models, expanding the kinematic and mass coverage of previous searches.

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: