Search for Higgs boson decays into Z and J/ ψ and for Higgs and Z boson decays into J/ ψ or Υ pairs in pp collisions at √ s = 13 TeV

Decays of the Higgs boson into a Z boson and a J/ ψ or ψ ( 2S ) meson are searched for in four-lepton final states with the CMS detector at the LHC. A data set of proton-proton collisions corresponding to an integrated luminosity of 138 fb − 1 is used. Using the same data set, decays of the Higgs and Z boson into quarkonium pairs are also searched for. An observation of such decays with this sample would indicate the presence of physics beyond the standard model. No evidence for these decays has been observed and upper limits at the 95% confidence level are placed on the corresponding branching fractions ( B ). Assuming longitudinal polarization of the Higgs boson decay products, 95% confidence level observed upper limits for B ( H → ZJ/ ψ ) and B ( H → Z ψ ( 2S )) are 1.9 × 10 − 3 and 6.6 × 10 − 3 , respectively. Published


Introduction
A boson with a mass of about 125 GeV was discovered by the ATLAS and CMS Collaborations at the CERN LHC in 2012 [1][2][3]. Comprehensive studies in various decay channels and production modes followed, and combined measurements show that the properties of the new boson are, so far, consistent with the standard model (SM) predictions for the Higgs boson (H) [4][5][6]. Rare exclusive decays of the H to mesons provide experimentally clean final states to study deviations of Yukawa couplings from SM predictions that cannot be obtained with inclusive measurements. Several models beyond the SM predict enhanced Yukawa couplings to fermions [7], leading to branching fractions (B) that are enhanced by up to three orders of magnitude. Examples include the Giudice-Lebedev model of quark masses [8], the two Higgs doublet model [9], the single Higgs doublet model with the Froggatt-Nielsen mechanism [10], and Randall-Sundrum models [11,12]. The required sensitivity for observing Yukawa couplings to the second-and first-generation fermions has not yet been reached, although recently the CMS Collaboration published evidence for the Higgs boson decay to a pair of muons [13,14]. The observed upper limit for B(H → cc ) times the cross section for H production in association with vector bosons as measured by the ATLAS and CMS experiments is found to be 26 and 14 times the SM expectation [15,16], respectively.
The first class of processes that is considered here is the decay of the Higgs boson into a Z boson and a vector meson quarkonium state (Q) [7,17,18]. The relevant SM Feynman diagrams for the decays H → ZQ are shown in Fig. 1 contributing amplitudes at leading order (LO), where the Higgs boson directly couples to a quark-antiquark pair that radiates a Z boson and forms the meson. The last two diagrams depict indirect contributions to the decay amplitude. The last diagram corresponds to oneloop amplitudes as indicated by the blob, as well as tree-level effective vertices [17]. The branching fractions of these processes in the SM are B(H → ZJ/ψ) = 2.3 × 10 −6 [17-19] and B(H → Zψ(2S)) = 1.7 × 10 −6 [19]. The main source of background events in these final states is from production of a Z boson in association with a genuine [20] or misidentified meson candidate. New physics could affect the direct boson couplings or could enter through loops and alter the interference pattern between the amplitudes. Any of those possibilities can enhance branching fractions with respect to the SM predictions. For the rare decays H → ZJ/ψ and H → Zψ(2S), the H → Zγ * amplitude contributes significantly. Hence, these rare decays provide complementary information to the decay H → Zγ, both in and beyond the SM [17, 19].

The CMS detector
The CMS apparatus [37] is a multipurpose, nearly hermetic detector, designed to trigger on and identify electrons, muons, photons, and hadrons. The CMS superconducting solenoid provides 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 (ECAL), 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. 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. [37].
The silicon tracker measures charged particles within the range |η| < 2.5. It consists of pixel and strip detector modules. An entirely new pixel detector has been installed during a technical stop between the 2016 and 2017 data-taking periods, featuring an all-silicon device with four layers in the barrel and three disks in the endcaps [38], providing four pixel detector measurements within a range |η| < 3.
Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The single-muon trigger efficiency exceeds 90% over the full η range, and the efficiency to reconstruct and identify muons is greater than 96%. Matching muons to tracks measured in the silicon tracker results in relative transverse momentum (p T ) resolutions of 1% in the barrel and 3% in the endcaps for muons with p T up to 100 GeV, and better than 7% in the barrel for muons with p T up to 1 TeV [39].
An electron is reconstructed by combining an energy measurement in the ECAL with a momentum measurement in the tracker. The ECAL consists of 75,848 lead tungstate crystals, which provide coverage of |η| < 1.48 in the barrel region and 1.48 < |η| < 3.00 in the two endcap regions (EE). Preshower detectors consisting of two planes of silicon sensors interleaved with a total of three radiation lengths of lead are located in front of each EE detector. The momentum resolution for electrons with p T ≈ 45 GeV from Z → e + e − decays ranges from 2 to 5%. It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [40].
The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [41,42] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector p T sum of those jets. Other collision vertices in the event are considered to have originated from additional inelastic pp collisions in each bunch crossing, referred to as pileup (PU). The average number of PU interactions during the 2016 data-taking period was 23, and increased to 32 during the 2017 and 2018 data-taking periods.
Events of interest are selected using a two-tiered trigger system. 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 fixed latency of about 4 µs [43]. 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 [44].
Dedicated triggers were deployed to enhance the selection of events of interest for the present study. Single-lepton (muon or electron) triggers are used to study the ZJ/ψ channel. The single muon trigger requires an isolated muon with p T > 27 GeV. The single electron trigger requires an isolated electron having p T > 27 (35, 32) GeV during the year 2016 (2017, 2018). The mesonspecific triggers are used to study quarkonium channels. They require the presence of at least three muons with p T > 2 GeV. Two of those must be oppositely charged and originate from a common vertex with a probability greater than 0.5%, as determined by a Kalman vertex filter [45], thus suppressing random combinations of two muons. The J/ψ-(Υ-) specific trigger requires a dimuon system invariant mass to be between 2.95 and 3.25 (8.5 and 11.4) GeV. The dimuon system p T is required to exceed 3.5 GeV for the years 2017 and 2018.

Simulated samples
Simulated samples of the H and Z boson signals are used to estimate the expected signal yields and model the distribution of signal events in the four-lepton invariant mass. The SM Higgs boson signals are simulated at next-to-leading order (NLO) in perturbative QCD with the POWHEG v2.0 Monte Carlo event generator [46,47], which includes the gluon-gluon fusion and vector boson fusion production processes. The parton distribution function (PDF) set used is NNPDF3.1 [48]. The JHUGen 7.1.4 generator [49,50] is used to decay the Higgs boson into Z bosons and Q mesons. The generator is interfaced with PYTHIA 8.226 [51] for parton showering, hadronization, and underlying event simulation using the CUETP8M1 [52] 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 [7].
The Z bosons are produced and decayed with the PYTHIA 8.226 generator [51] which implements the LO matrix element calculation interfaced with parton showering, hadronization, and underlying event simulation for which the tune CUETP8M1 [52] is used. For the calculation of the branching fraction, the inclusive SM Z boson production cross section is obtained according to Ref. [53], where the prediction includes the next-to-next-to-leading order (NNLO) QCD contributions, and the NLO electroweak corrections from FEWZ 3.1 [54] calculated using the NLO PDF set NNPDF3.0. The generated events are then reweighted to match the p T -spectrum of the Z boson predicted at NLO [46, 47,55].
The generated events are processed through a detailed simulation of the CMS detector based on GEANT4 [56]. Simulated events include additional pp interactions. Events are then reweighted to match the PU profile observed in data.
The acceptance of the final states changes with the angular distribution of leptons in the decay. The distribution of the decay angle θ, defined as the angle between the positive lepton momentum in the rest frame of the intermediate particle (J/ψ meson or Υ meson or Z boson) with respect to the direction of this intermediate particle in the rest frame of the parent particle (H or Z boson), is proportional to (1 + λ θ cos 2 θ), where λ θ is the average polar anisotropy parameter [57]. According to Refs. [27,35,58], longitudinal polarization is expected for the Z boson and assumed for mesons. In this Letter, the nominal results are obtained using a signal acceptance calculated for the longitudinally polarized case (λ θ = −1). Two other extreme scenarios where the intermediate particles are both either fully transversely polarized (λ θ = +1) or unpolarized (λ θ = 0) are also considered. No azimuthal anisotropies are considered.

Event reconstruction
In a first step, events with at least 2µ and two additional ℓ (ℓ = e or µ) are selected. The J/ψ and Υ candidates are built from µ + µ − pairs, and Z candidates from µ + µ − or e + e − pairs.
Muons are reconstructed by combining information from the silicon tracker and the muon system [39]. 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 detector layers 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 detector layer. 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. To suppress muons originating from non-prompt 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 (along the longitudinal axis). Events with at least four such muons with p T > 3 GeV and |η| < 2.4 and zero sum of charges are accepted. To measure the isolation of the leading (highest p T ) 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 p T sum of the reconstructed innerdetector tracks originating from the primary pp interaction vertex within the cone has to be less than 50% of the muon's p T . The leading muon p T is excluded from the sum and the subleading muon p T is also excluded if this muon falls within the isolation cone of the leading muon.
The energy of electron candidates is determined from a combination of the track momentum at the main interaction vertex, the corresponding ECAL cluster energy, and the energy sum of all bremsstrahlung photons attached to the track [40]. The dielectron mass resolution for Z → e + e − decays when both electrons are in the ECAL barrel is 1.9%, and is 2.9% when they are in the endcaps. To reduce contamination from particles incorrectly identified as electrons, reconstructed electrons are required to pass a multivariate electron identification discriminant. This discriminant, which is described in Ref. [40], combines information about the quality of tracks, the shower shape, kinematic quantities, and hadronic activity in the vicinity of the reconstructed electron. Isolation variables are also included among the discriminant's inputs. Therefore, no additional isolation requirements are applied. The selection based on the multivariate identification discriminant has an electron identification efficiency of 90% while the rate of misidentifying other particles as electrons is 2-5%. Events with at least two oppositely charged electrons with p T > 3 GeV and |η| < 2.5 and two oppositely charged muons are accepted.
Further selection criteria are applied to the different four-lepton final states to achieve the lowest expected upper limit at 95% CL [59][60][61]. This optimization is performed with data with the signal regions removed [62] and replaced by simulated events, and with the simulated signal shape.
The selection in the ZJ/ψ channel requires each dilepton resonance to have p T > 5 GeV, and the candidate invariant mass to lie in the region 80-100 GeV for the Z boson (3.0-3.2 GeV for the J/ψ meson). The dimuon mass resolution is about 1%. Each dilepton must fit to a common vertex with probability greater than 1%, determined by a vertex fit probability. The four-lepton candidate must have p T > 5 GeV, and the fit to a common four-lepton vertex must have a probability of greater than 1%. A total of 230 (177) candidate events are found in the 4µ (2e2µ) invariant mass window of 112-142 GeV. The lower limit of the four-lepton invariant mass (m 4ℓ ) range is chosen to exclude the region close to the ZJ/ψ threshold where m 4ℓ changes rapidly, which is difficult to model. The selection criteria for the decay H → Zψ(2S), where ψ(2S) decays inclusively into J/ψ, are identical. The respective m 4ℓ distributions are shown in Fig. 2. Non-resonant background in the dilepton invariant mass distributions is found to contribute about 20%. These events are an additional source of background which does not concentrate in a specific region of m 4ℓ , and thus is accounted for by an empirical parameterization of the background in the fit as described in the next section.   In the case of the J/ψ pair channel, each dimuon has to be fit to a common vertex with a probability greater than 0.5%. In addition, the J/ψ candidate must have p T > 3.5 GeV, matching the trigger requirement, and the invariant masses of the higher and lower-p T J/ψ candidates have to be within 0.10 and 0.15 GeV, respectively, of the nominal mass of the J/ψ meson. The mass window of the subleading J/ψ is wider to allow monitoring the reduction of the sideband population as the selection progresses. In order to suppress contributions from non-prompt hadrons, separately produced J/ψ mesons and muons from other sources, the four-muon vertex fit probability of J/ψ pairs must be greater than 5%. Finally, the absolute value of the difference in rapidity (|∆Y|) between the two J/ψ candidates has to be less than 3. This criterion marginally affects the signal, while removing about 20% of the background events. After the selection, 720 events are found in data in the 40-140 GeV four-muon invariant mass (m 4µ ) range. Nonresonant background in the dilepton invariant mass distributions is found to be negligible. Figure 3 shows the m 4µ distribution of the J/ψJ/ψ candidates.
An Υ pair candidate event must have at least four muons each with p T > 4 GeV. The Υ(nS)(n = 1, 2, 3) and Υ(1S) candidates are formed with oppositely charged muon pairs with p T > 5 GeV, and the dimuon invariant mass within the range 9.0-10.7 GeV and 9.0-9.7 GeV, respectively. In order to suppress random combinations, dimuon and four-muon objects are required to have a vertex fit probability greater than 1%. Between two candidate dimuons, the |∆Y| value has to be less than 2.3 and the azimuthal angle difference has to be greater than 1 radian. The four-muon combination must have p T > 5 GeV and an absolute rapidity of less than 1.7. After applying the selection criteria in data, 59 Υ(nS)Υ(mS) (18 Υ(1S)Υ(1S)) candidate events are found in the 40-140 GeV m 4µ range. Figure 4 left and right show the m 4µ distributions for Υ(nS)Υ(mS) and Υ(1S)Υ(1S) candidates, respectively.
The differences in the efficiencies between data and simulation for the trigger, offline lepton reconstruction, identification, and isolation requirements are corrected by reweighting the simulated events with data-to-simulation correction factors, which are obtained with the "tag- and-probe" method [63] using J/ψ → µ + µ − and Z → ℓ + ℓ − events. The total signal efficiency, including kinematic acceptance, trigger, reconstruction, identification, and isolation efficiencies, for the ZJ/ψ → 4µ (ZJ/ψ → 2e2µ) decays with longitudinally polarized J/ψ and Z, is found from simulation to be around 30 (24)%. For the Higgs boson decays into Υ(nS)Υ(mS) and J/ψJ/ψ, the corresponding total efficiencies are about 31 and 30%, respectively. For the Z boson, the corresponding values are about 28 and 32%.

Signal extraction and systematic uncertainties
Unbinned extended maximum likelihood fits [64] to the m 4ℓ distributions are performed. Yields for signals and backgrounds are free parameters in the fit. The background shapes in the m 4ℓ distributions are obtained from data and are described by an exponential plus constant function. For the ZJ/ψ → 4µ channel, the H signal is parameterized with a sum of a Gaussian and a Crystal Ball function [65], and for the ZJ/ψ → 2e2µ channel with two Crystal Ball functions. Similarly, the H signal in the J/ψ pair channel is described with a double-Gaussian function, and in the Υ pair channel with a combination of a Gaussian and Crystal Ball functions. In all combined functions the mean is a common parameter. The simulated Z boson signal is described with a Voigtian function with the resonance width fixed to the world-average value [66]. The mass resolution and the mean are taken from the fit to the simulation, and they are fixed in the fit to data. The background function from the fit to data in the ZJ/ψ and J/ψJ/ψ channels are superimposed as solid blue lines in Figs. 2 and 3, respectively. In the Υ pair sample, no events are observed above the m 4µ of 80 GeV. The m 4µ distribution below 80 GeV is well described solely by an exponential function. Figure 4 shows the observed m 4µ distribution with the fit superimposed. ii) The differences between data and simulation for the trigger, offline muon reconstruction, identification, and isolation efficiencies are corrected by reweighting the simulated events with data-to-simulation correction factors, which are obtained with the tag-and-probe method using J/ψ (or Z) → µ + µ − events. The resulting scale factors in muon identification, isolation, and trigger efficiencies are observed to deviate from unity by less than 2(2), 0.5(0.5) and 1(3)% for the ZJ/ψ (QQ) channel. Analogously, the uncertainty in the electron reconstruction, identification and trigger efficiency is found to be about 2% [40].
iii) The relative difference in the four-lepton vertex criterion between data and simulation is evaluated with ZJ/ψ (J/ψ pair) event samples. It is found to be less than 2(3)% for the ZJ/ψ (QQ) channel. iv) Differences in the lepton momentum scale and resolution in data and simulation are estimated from J/ψ and Z dilepton signals and extrapolated to the four-lepton signals. The systematic uncertainty is estimated as the relative change in the upper limit when varying the signal mass mean and width by these differences. They are found to be less than 1(3)% in the 4µ (2e2µ) channel.
vi) A common parameterization for each signal model is used for the entire run period. The relative uncertainty in the signal model due to the change in detector conditions in each year is determined to be 1(2)% for the ZJ/ψ (QQ) channels.
vii) The background is alternatively parameterized with a second order Chebyshev polynomial or a power law function. The relative uncertainty due to the choice of the background function is found to be negligible.

Results
No evidence for the Higgs or Z boson signal is found in any of these channels. The results of this analysis are presented as upper limits on the branching fractions and are set at the 95% CL. Limits are determined with the modified frequentist CL s criterion, in which the profile likelihood ratio modified for upper limits is used as the test statistic [59][60][61]. Systematic uncertainties are incorporated in the likelihood as nuisance parameters. The observed and median expected exclusion limits for the branching fractions at 95% CL for the H and Z boson decays are listed in Table 1. Figures 2-4 show the distributions for simulated boson signals in different search channels as dashed lines. The signals are normalized to their observed 95% upper limit branching fractions.
The results for B(H → ZJ/ψ) and B(H → Zψ(2S)) each are obtained by combining the channels with Z → e + e − and Z → µ + µ − . The values for B(J/ψ → µ + µ − ), B(Z → ℓ + ℓ − ) and B(ψ(2S) → J/ψ) are taken from Ref. [66]. This analysis does not distinguish between the three Υ(nS) states. To calculate their contribution to the corresponding H and Z boson branching fractions, the coupling strength of the bosons to any Υ(nS) pairing is assumed to be the same. All Υ states can directly decay into muon pairs with the different branching fractions taken from Ref. [66]. In addition, it is assumed that the Υ states could be the result of a one step transition Υ(3S) → Υ(2S), Υ(3S) → Υ(1S), or Υ(2S) → Υ(1S) before decaying into muons [66]. Consequently, in the Υ(1S)Υ(1S) channel, the feed-down transitions from Υ(3S) and Υ(2S) to Υ(1S) are included. The observed upper limit branching fractions at 95% CL agree with the expected limits. In the case of H → ZJ/ψ, the upper limit branching fraction is about 800 times higher than the SM prediction [17]. In the H → Υ(nS)Υ(mS) channel, the observed upper limit branching fraction is found to be about one order of magnitude higher than SM predictions that assume dominance of direct quark couplings [23]. The factors are larger for all other channels. Tabulated results are available in the HEPData record for this analysis [74].