Search for Higgs and Z boson decays to J/$\psi$ or $\Upsilon$ pairs in proton-proton collisions at $\sqrt{s}=$ 13 TeV

A search for decays of the Higgs and Z boson to pairs of J/$\psi$ or $\Upsilon$(nS) (n=1, 2, 3) mesons, with their subsequent decay to $\mu^+\mu^-$ pairs, is presented. The analysis uses data from proton-proton collisions at $\sqrt{s}=$ 13 TeV, collected with the CMS detector at the LHC in 2017 and corresponding to an integrated luminosity of 37.5 fb$^{-1}$. While an observation of such a decay with this sample would indicate the presence of physics beyond the standard model, no significant excess is observed. Upper limits at 95% confidence level are placed on the branching fractions of these decays. In the J/$\psi$ pair channel, the limits are 1.8$\times$10$^{-3}$ and 2.2$\times$10$^{-6}$ for the Higgs and Z boson, respectively, while in the combined $\Upsilon$(nS) pair channel, the limits are 1.4$\times$ 10$^{-3}$ and 1.5$\times$10$^{-6}$, respectively, when the mesons from the Higgs and Z boson decay are assumed to be unpolarized. When fully longitudinal and transverse polarizations are considered the limits reduce by about 22-29% and increase by about 10-13%, respectively.


Introduction
A new boson with a mass of 125 GeV was discovered by the ATLAS and CMS Collaborations at the CERN LHC in 2012 [1][2][3][4][5][6][7]. Comprehensive studies in various decay channels and production modes have shown that the properties of the new boson are consistent, so far, with expectations for the standard model (SM) Higgs boson (H) [7][8][9]. Recently, the Higgs boson couplings to top and bottom quarks have been directly measured [10][11][12][13]. Couplings to lighter quarks are still not observed directly. Rare exclusive decays of the Higgs boson to mesons provide experimentally clean final states to study Yukawa couplings to quarks and physics beyond the SM (BSM). Examples of diagrams for decays of the Higgs boson into quarkonium pairs according to Ref. [14] are displayed in Fig. 1 (three leftmost plots). The symbol Q refers to charmonium and bottomonium states.  Figure 1: The three leftmost plots show Feynman diagrams for H → QQ with Q = charmonium or bottomonium states according to Ref. [14]. In the two leftmost diagrams, the virtual particles are Z bosons. In center-right diagram, quarks are the main contribution to the loops and the virtual particles are either photons or gluons. In the latter case additional soft-gluon exchange occurs. The rightmost diagram shows the leading order Z → QQ decay diagram according to Ref. [15].
The importance of the measurement of such decays has been pointed out by Ref. [16][17][18][19]. Using a phenomenological approach for the H-qq coupling, Ref. [16] finds that the dominant quarkonium pair decay mode is H → ΥΥ and estimates its branching fraction (B) to be at the level of 10 −5 . The early calculations of Higgs boson decays into a pair of heavy quarkonia states did not include relativistic corrections caused by the internal motion of quarks [14]. The importance of the latter corrections is underlined by the fact that the predicted e + e − → J/ψη c cross section increases by an order of magnitude [20][21][22] when these effects are included, in agreement with measurements by the Belle and BaBar experiments [23,24].
With emphasis on amplitudes where the Higgs boson couples indirectly to the final state mesons, such as represented by the two leftmost diagrams in Fig. 1, Ref. [14] arrives at values of about B(H → J/ψJ/ψ) = 1.5 × 10 −10 and B(H → ΥΥ) = 2 × 10 −9 for the Higgs boson. The mechanism where the Higgs boson couples directly to charm or bottom quarks, which then hadronize to heavy quarkonia, was considered in a recent calculation [25] leading to an increase of an order of magnitude in B(H → J/ψγ). Recently, this decay has been searched for by the ATLAS and CMS collaborations [26,27]. The Higgs boson decay to the J/ψ pair could also occur when the photon in the J/ψγ decay is virtual and transforms into a J/ψ meson. This Letter also presents the first search for decays of the Z boson into quarkonium pairs. A leading order Feynman diagram is shown in Fig. 1 (rightmost plot). The SM prediction for B(Z → J/ψJ/ψ) is of the order of 10 −12 in nonrelativistic QCD and leading twist light cone models [15].
New physics could affect the direct H-qq couplings or could enter through loops, and alter the interference pattern between the amplitudes. Any of those possibilities enhance branching fractions with respect to the SM predictions. Many BSM theories predict substantial modifications of the Yukawa couplings of the Higgs boson to quarks, such as models with Higgsdependent Yukawa couplings [28], the minimal flavor violation framework [29], the Froggatt-Nielsen mechanism [30], and the Randall-Sundrum family of models [31]. An overview of models can be found in Ref. [32]. In the related quarkonium-γ channels, deviations of the Hqq couplings from the SM predictions can change the interference between direct and indirect amplitudes, resulting in substantial modifications of the branching fractions, particularly in the Υ channel, where the increase is up to several orders of magnitude [25]. The observation of a Higgs or Z boson signal in the quarkonium pair decay modes with the available LHC data sets would indicate the presence of BSM physics.
This Letter presents the first search for the Higgs and Z boson decays into J/ψ or Υ meson pairs, where Υ stands for the combined contribution of the Υ(nS) states with n = 1,2,3. The subsequent decay of these meson pairs to the 4µ final state offers a very clean experimental signature that is used in this analysis. For the J/ψ meson pairs, feed-down from higher charmonium states are not taken into account. For the Υ(nS) meson pairs, decays from higher to lower mass Υ(nS) states are included. The results presented in this Letter are based on proton-proton (pp) collision data recorded in 2017 with the CMS detector at a center-of-mass energy of √ s = 13 TeV, amounting to an integrated luminosity of 37.5 fb −1 .

The CMS detector
A 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. [33]. The central feature of the CMS apparatus is a superconducting solenoid, 13 m in length and 6 m in internal diameter, providing an axial 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. 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. They are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers.
An entirely new pixel detector has been installed after 2016, featuring an all-silicon device with four layers in the barrel and 3 disks in the endcaps [34], providing four pixel detector measurements and reduced material budget in front of the calorimeters.
Events of interest are selected using a two-tiered trigger system [35]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
Dedicated triggers were deployed to enhance the events of interest for the present study. They require the presence of at least three muons with p T greater than 2 GeV. Two of these must be oppositely charged and have to originate from a common vertex with a probability greater than 0.5%, as determined by a Kalman vertex fit [36]. The J/ψ-specific trigger requires a dimuon system's invariant mass to be between 2.95 and 3.25 GeV and its p T to be greater than 3.5 GeV. The trigger used to select the Υ sample requires two of the three muons to have p T greater than 3.5 GeV, and one muon p T greater than 5 GeV. The invariant mass for one oppositely charged muon pair must lie in the interval 8.5-11.4 GeV. Both triggers gave an efficiency exceeding 85% to select events satisfying the selection criteria used in the analysis.

Signal and background modelling
Simulated samples of the Higgs and Z boson signals are used to estimate the expected signal yields and model the distribution of signal events in the four-muon invariant mass. For the H → J/ψJ/ψ and H → ΥΥ samples the Higgs boson is produced with the POWHEG v2.0 Monte Carlo (MC) event generator [37,38], which includes the gluon-gluon fusion (ggF) and vector-boson fusion production processes. The parton distribution function (PDF) set used is NNPDF3.1 [39]. The JHUGen 7.1.4 generator [40,41] is used to decay the Higgs boson into two vector mesons taking into account their helicity. To produce the decay for unpolarized quarkonia, the JHUGen generator is configured to model a uniform muon helicity angle distribution. The generator is interfaced with PYTHIA 8.226 [42] for parton-showering and hadronization according to the CUETP8M1 [43] tune. The total SM Higgs boson production cross section for the calculation of branching fractions is taken from the LHC Higgs cross section working group [32].
The Z → J/ψJ/ψ and Z → ΥΥ samples are produced with the PYTHIA 8.226 generator [42], tune CUETP8M1 [43]. The SM Z boson production cross section includes the next-to-next-toleading order (NNLO) QCD contributions, and the next-to-leading order (NLO) electroweak corrections from FEWZ 3.1 [44] calculated using the NLO PDF set NNPDF3.0. The Z boson p T is reweighted to match the NLO calculation [37,38,45]. The total cross section is obtained with In the J/ψ and Υ pair channels backgrounds are assumed to originate from prompt nonresonant pair production, which in pp collisions dominantly occurs via ggF [47-50]. Initially, the two mesons are color-octet bound states that then radiate soft gluons to become real mesons. Event samples are generated according to this model [49].
The generated events are processed through a detailed simulation of the CMS detector based on GEANT4 [51]. The high instantaneous luminosity of the LHC results in multiple pp interactions per bunch crossing. Simultaneous pp interactions that overlap with the event of interest, i.e. pileup, are included in simulated samples. The distribution of the number of additional interactions per event in simulation corresponds to that observed in the data.
The acceptance of the final states changes with the angular distribution of the muons in the quarkonium decay. The distribution of the decay angle θ, defined as the angle between the positive muon direction of flight in the rest frame of the quarkonium with respect to the quarkonium direction in the boson rest frame, is proportional to (1 + λ θ cos 2 θ). In this Letter, the nominal results are obtained using a signal acceptance calculated for the unpolarized case (λ θ = 0). Two extreme scenarios have also been considered, where the J/ψ and Υ mesons are either fully transversely polarized, λ θ = +1, or fully longitudinally polarized, λ θ = −1. No azimuthal anisotropies have been considered. According to Refs. [14, 15] the J/ψ and Υ mesons produced in the decays of both bosons are expected to be dominantly longitudinally polarized.

Data reconstruction and selection
Muons are reconstructed by combining information from the silicon tracker and the muon system [52]. The matching between tracks reconstructed in each of the subsystems proceeds either outside-in, starting from a track in the muon system, or inside-out, starting from a track provided by the silicon tracker. In the latter case, tracks that match track segments in only one or two stations of the muon system are also considered in the analysis to collect very low-p T muons that may not have sufficient energy to penetrate the entire muon system. The muons are selected from the reconstructed muon track candidates that match with at least one segment in any muon station in both x and y. The number of silicon tracker layers with hits used in the muon track candidate has to be greater than 5 and include at least one pixel detector layer. Matching muons to tracks measured in the silicon tracker results in a relative transverse momentum (p T ) resolution of 1% in the barrel and 3% in the endcaps for muons with p T up to 100 GeV. The p T resolution in the barrel is better than 7% for muons with p T up to 1 TeV [52].
The reconstructed vertex with the largest value of summed charged particle p 2 T is taken to be the primary pp interaction vertex. To suppress muons originating from nonprompt hadron decays, the impact parameter of each muon track, computed with respect to the position of the primary pp interaction vertex, is required to be less than 0.3 (20.0) cm in the transverse plane (longitudinal axis). Events with at least four such muons with p T > 3 GeV and |η| < 2.4 are accepted. To isolate the leading muon candidate from other hadronic activity in the event, a cone of size ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 is constructed around its momentum direction, where φ is the azimuthal angle in radians. The sum of the p T of the reconstructed inner-detector tracks originating from the primary pp interaction vertex within the cone has to be less than 50% of the muon's p T . The transverse momentum of the leading muon is subtracted from the sum and the subleading muon p T is also subtracted, if this muon falls within the isolation cone of the leading muon.
The J/ψ and Υ candidates are built from pairs of oppositely charged muons. Each muon pair must fit to their common vertex with a probability greater than 0.5%. The J/ψ candidate p T has to be greater than 3.5 GeV, and the invariant mass of the higher and lower-p T J/ψ candidates have to be within 0.1 and 0.15 GeV, respectively, of the nominal mass of J/ψ. The dimuon mass resolution is about 30 MeV. To suppress contributions from nonprompt hadrons, separately produced J/ψs and muons from other sources, the four-muon Kalman vertex fit probability of J/ψ pairs has to be greater than 5%. Finally, the absolute value of the difference in rapidity between the two J/ψ candidates has to be less than 3. This criterion marginally affects the signal while removing about 20% of the selected events. After the selection, 189 events are found in data set in the 40-140 GeV four-muon invariant mass range. Figure 2 (left) shows the four-muon invariant mass distribution.
For the selection of Υ pair candidates, the same event selection criteria are applied, except that the Υ candidate p T has to be greater than 5 GeV, and the invariant mass has to fall within the range 8.5-11 GeV. Furthermore, the four-muon Kalman vertex fit probability has to be greater than 1% to suppress random combinations. The nonprompt background is negligible in this channel. After applying the selections, 106 events are found in data set in the 20-140 GeV fourmuons invariant mass range. Figure 2 (right) shows the four-muon invariant mass distribution.
The differences in efficiencies between data and simulation for the trigger, offline muon reconstruction, identification, and isolation are corrected by reweighting the simulated events with data-to-simulation correction factors, which are obtained with the "tag-and-probe" method [53] using J/ψ → µ + µ − events. The scale correction factors are observed to deviate from unity by less than 3%. The difference in the four-muon Kalman vertex fit efficiency between data and simulation is evaluated with J/ψ pair event samples and found to be less than 3%. The total signal efficiency, including kinematic acceptance, trigger, reconstruction, identification, and isolation efficiencies, for the J/ψJ/ψ decays with unpolarized J/ψ is approximately 23% for both bosons. For the ΥΥ decays the corresponding efficiency is about 27%.

Results
Unbinned extended maximum-likelihood fits [54] to the four-muon invariant mass distributions M 4µ are performed. Yields for signals and backgrounds are free parameters in the fit. For the Higgs boson the invariant mass distribution obtained from simulation is described with two Gaussian functions with a common mean. The simulated Z signal is described with a Voigtian function with the world-average value for the resonance width [46]. The mass resolution and mean are taken from the fit to the simulation, and they are fixed in the fit to data.
The four-muon invariant mass distribution up to 140 GeV is described by an exponential plus constant function. The relative contribution and decay constant of the exponential function are varied in the fit to data. The values of both parameters are found to be in close agreement between observation and simulation [49]. The result of the fit is shown as a solid blue line in Fig. 2 (left).
In the Υ pair sample, no events are observed above the four-muon invariant mass of 40 GeV. The four-muon invariant mass distribution is modeled analogously to the J/ψ pair channel. The M 4µ distribution below 40 GeV is well described by an exponential function. The decay constant of the exponential function is also varied in the fit. The same function describes an event sample generated with the pair production model [49]. The observed and median expected exclusion limits for the branching fractions at 95% confidence level (CL) for the H and Z boson decays listed in Table 1.
The relative changes in the upper limits on the Higgs boson decay branching fractions with respect to the case of unpolarized decay mesons are about −22% for fully longitudinally polarized J/ψ and Υ mesons, and +10% for fully transversely polarized mesons. For the Z boson the relative changes are about −29 (−26)% for fully longitudinally polarized J/ψ (Υ) mesons and +13 (+12)% for fully transversely polarized mesons.

Summary
In summary, this Letter presents the first search for decays of the Higgs and Z boson to pairs of J/ψ or Υ(nS) (n=1,2,3) mesons, with their subsequent decay to µ + µ − pairs. Data from pp collisions at √ s = 13 TeV, corresponding to an integrated luminosity of 37.5 fb −1 are used. No excess has been observed above a small background in the J/ψ pair and with vanishingly small background in the Υ pair channels. The observed upper limits at 95% confidence level on the branching fractions for the Higgs boson decays for unpolarized mesons are B(H → J/ψJ/ψ) < 1.8 × 10 −3 and B(H → ΥΥ) < 1.4 × 10 −3 . The observed upper limits on the branching fractions for the Z boson decay in the unpolarized scenario are B(Z → J/ψJ/ψ) < 2.2 × 10 −6 and B(Z → ΥΥ) < 1.5 × 10 −6 , where all three Υ(nS) states are considered. Extreme polarization scenarios give rise to variations in the observed boson decay branching fractions between −(22-29)% for fully longitudinally polarized J/ψ and Υ mesons and +(10-13)% for fully transversely polarized mesons.
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: [6] ATLAS Collaboration, "Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment", Eur. Phys. J. C 76 (2016) 6, doi:10.1140/epjc/s10052-015-3769-y, arXiv:1507.04548.