Search for Higgs boson decays into Z and J / ψ and for Higgs and Z boson decays into J / ψ or Y pairs in pp collisions at √

Decays of the Higgs boson into a Z boson and a J / ψ or ψ ( 2S ) meson are searched for in four-lepton ﬁnal states with the CMS detector at the LHC. A data set of proton-proton collisions corresponding to an integrated luminosity of 138fb − 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% conﬁdence level are placed on the corresponding branching fractions ( B ). Assuming longitudinal polarization of the Higgs boson decay products, 95% conﬁdence level observed upper limits for B ( H → ZJ / ψ ) and B ( H → Z ψ ( 2S )) are 1 . 9 × 10 − 3 and 6 . 6 × 10 − 3 , respectively. © 2022 The Author(s). Published by Elsevier B.V. This is an open access article


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.E-mail address: cms -publication -committee -chair @cern .ch. quark coupling contributions to the H → ZQ decay, where Q represents a quarkonium state and q is a quark, which in the leftmost diagram is either of charm or bottom flavor.The diagrams represent Higgs boson decays into quarkonium pairs when replacing the bottom section with the upper half in each.
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.The first diagram in Fig. 1 represents 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 one-loop 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][18][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 am-

The CMS Collaboration
Physics Letters B 842 (2023) 137534 plitudes.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 decays of the Higgs boson into ZJ/ψ and Zη c have been searched for by the ATLAS Collaboration in hadronic final states, reaching 95% confidence level (CL) upper limits on the branching fraction of the Higgs boson that exceed 100% [21].Recently, the CMS Collaboration published upper limits on the branching fraction for H → Zρ and H → Zφ at 95% CL that are larger than the expected SM branching fraction by more than a factor of 700 [22].
A second related class of processes is the Higgs boson decay into pairs of quarkonia.The Feynman diagrams are variants of the ones depicted in Fig. 1: in each diagram, the on-shell Z boson in the lower part is replaced by a quarkonium decay, similar to the process depicted in the upper part.In the rightmost diagram, both vector bosons could be also gluons, in which case additional soft-gluon exchange occurs.The importance of the measurement of such decays has been pointed out in Refs.[23][24][25][26].Using a phenomenological approach for the direct H-qq coupling, Ref. [23] finds that the dominant quarkonium pair decay mode is H → ΥΥ with an estimated branching fraction of O (10  −5 ).More recently, Ref. [27] predicts values of B(H → J/ψJ/ψ) = 1.5 × 10 tudes.Inclusion of the mechanism where the Higgs boson couples directly to charm or bottom quarks, which then hadronize to heavy quarkonia, in the calculation in Ref. [28] leads to an increase by an order of magnitude in the related B(H → J/ψγ ).The Higgs boson is expected to couple to quarkonium pairs that include radially excited states with comparable strength [29].Recently, a first search for the decays H → J/ψJ/ψ and H → Υ(nS)Υ(mS)(n, m = 1, 2, 3) was performed by the CMS Collaboration [30].Related to these two classes of processes is the decay of the Higgs boson into a photon and a vector meson [17,28,31].The 95% CL upper limits on the branching fractions of the Higgs boson into γ J/ψ, γ ρ, and γ φ are found to be two orders of magnitude larger than their expected values in the SM [32][33][34].For the γ ψ(2S) and γ Υ(nS) decays, the corresponding upper limits are, respectively, three and five orders of magnitude larger than the expected SM branching fractions.
This Letter presents the first search for decays of the 125 GeV Higgs boson into a Z boson and a J/ψ or ψ(2S) meson in fourlepton final states.The Z boson is reconstructed from its decays into μ + μ − or e + e − , the J/ψ meson from its decay into μ + μ − , and the ψ(2S) from its inclusive decay into J/ψ X (feed-down), where X , which is mostly ππ, is not reconstructed.Furthermore, an update of Higgs boson searches in J/ψJ/ψ and Υ(nS)Υ(mS) decay channels with the full available data sample is presented.New channels are accessed via the inclusive decay of ψ(2S) into a J/ψ meson.For the Υ(nS) states the possibilities that they are the result of feed-down transitions from higher Υ states before decaying into muon pairs are included.Finally, this Letter also presents the search for decays of the Z boson into quarkonium pairs.The SM prediction for B(Z → J/ψJ/ψ) calculated in the framework of nonrelativistic quantum chromodynamics (QCD) and leading twist light cone models is of the order of 10 −12 [35].The Z decay into a J/ψ meson and a lepton pair, which is dominated by the electromagnetic fragmentation process, was observed at a rate consistent with SM predictions [36].
The results presented in this Letter are based on proton-proton (pp) collision data recorded in 2016-2018 with the CMS detector at a center-of-mass energy of in the quarkonium pair channels, where the second number is slightly smaller due to a delayed trigger deployment.

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 allsilicon 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 physicsobject 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 datataking 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].

The CMS Collaboration
Physics Letters B 842 (2023) 137534 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.4generator [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 fewz3.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 insideout, 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 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 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.2GeV 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 are 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.Non-resonant Fig. 3.The four-muon invariant mass distribution of J/ψJ/ψ candidates (error bars for empty bins are not shown).The result of the maximum likelihood fit to the background (Bkg) is superimposed (solid blue line).For illustrative purposes, the plots show the distributions for simulated signals of the Higgs and Z boson decaying into J/ψJ/ψ (dashed and dashed-dotted red lines).The signals for the Higgs boson decays H → ψ(2S)J/ψ (dotted magenta line) and H → ψ(2S)ψ(2S) (dashed-dotted black line) are also shown, where ψ(2S) decays into J/ψ .Each signal is normalized to their observed 95% CL upper limit branching fraction from this analysis.background in the dilepton invariant mass distributions is found to be negligible.Fig. 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.Fig. 4 upper and lower 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 Fig. 4. The four-muon invariant mass distributions for Υ(nS)Υ(mS) (upper) and Υ(1S)Υ(1S) (lower) candidates (error bars for empty bins are not shown).The result of the maximum likelihood fit to background (Bkg) is superimposed (solid blue line).For illustrative purposes, the plots show the distributions for simulated Higgs and Z boson signals (dashed and dashed-dotted red lines) normalized to their observed 95% CL upper limit branching fractions from this analysis.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 is 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.Fig. 4 shows the observed m 4μ distribution with the fit superimposed.
Separate fits are performed to the m 4 distributions for the different signal hypotheses.Signal shapes for the Higgs boson in decays involving the inclusive transition from ψ(2S) to J/ψ meson are modeled with a combination of the same functions as used for the Higgs boson directly decaying into ground state mesons (direct signal).For the fits to the feed-down channels, the background functions are identical and parameters are fixed to the ones from the previous direct signal fits.
Systematic uncertainties originate from imperfect knowledge of the detector and imperfect signal modeling.Systematic uncertainties considered in this analysis are listed below.
i) The integrated luminosities for the 2016, 2017, and 2018 datataking years have 1.2-2.5% individual uncertainties [67][68][69], while the overall uncertainty for the 2016-2018 period is 1.6%.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.v) The theoretical uncertainties in the production cross section for the H (Z boson) are ±3.2%(±1.7%) due to the choice of the PDF and the value of the strong coupling constant [7,48,70], and +4.6 −6.7 % (±3.5%) due to the renormalization and factorization scale choice [70][71][72][73].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.viii) The uncertainties in the Z, J/ψ, and Υ branching fractions to lepton pairs, and B(ψ(2S) → J/ψ + X) are taken from Ref. [66].

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.Figs. 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/ψ → μ 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].

Summary
This Letter presents the first search for decays of the Higgs boson (H) into a Z boson and a J/ψ meson in four-lepton final states.Data from proton-proton collisions at The observed upper limits for the Higgs and Z boson decays for longitudinally polarized mesons are B(H → ZJ/ψ) < 1.9 × 10 and B(Z → Υ(1S)Υ(1S)) < 1.8 × 10 −6 .The observed upper limit branching fraction for H → ZJ/ψ is about 800 times the value predicted by the standard model [17][18][19].For H → Υ(nS)Υ(mS) the upper limit is about one order of magnitude higher than predicted by earlier standard model calculations [23].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

√ s = 13
TeV, corresponding to an integrated luminosity of about 138 fb −1 , are used.Using the same data, decays of the Higgs and Z boson into quarkonium pairs are also searched for.No excess of a Higgs or Z boson signal above background is found in any of the searched channels and upper limits on branching fractions (B) at the 95% confidence level for various polarization scenarios are set.The Higgs boson decay is also searched for in channels where, before decaying into muon pairs, one or both J/ψ mesons could be the result of an inclusive ψ(2S) to J/ψ transition, and the Υ(nS)(n = 1, 2) mesons could be the result of inclusive transitions from Υ(nS)(n = 2, 3) mesons.

Table 1
Exclusion limits at 95% CL on the branching fractions of the H and Z boson decays.The second column lists the observed limits for the case that both intermediate particles are longitudinally polarized (λ θ = −1) as described in the text.The third column shows the median expected limits with the upper and lower bounds in the expected 68% CL intervals.The last two columns list observed upper limits for unpolarized (λ θ = 0) and transversely polarized (λ θ = +1) intermediate particles.
the Hungarian Academy of Sciences, the New National Excellence Program -ÚNKP, the NKFIH research grants K 124845, K 124850, K 128713, K 128786, K 129058, K 131991, K 133046, K 138136, K 143460, K 143477, 2020-2.2.1-ED-2021-00181, and TKP2021-NKTA-64 (Hungary); the Council of Science and Industrial Research, India; the Latvian Council of Science; the Ministry of Education and Science, project no.2022/WK/14, and the National Science Center, contracts Opus 2021/41/B/ST2/01369 and 2021/43/B/ST2/01552 (Poland); the Fundac ¸ão para a Ciência e a Tecnologia, grant CEECIND/01334/2018 (Portugal); the National Priorities Research Program by Qatar National Research Fund; MCIN/AEI/10.13039/501100011033,ERDF "a way of making Europe", and the Programa Estatal de Fomento de la Investigación Científica y Técnica de Excelencia María de Maeztu, grant MDM-2017-0765 and Programa Severo Ochoa del Principado de Asturias (Spain); the Chulalongkorn Academic into Its 2nd Century Project Advancement Project, and the National Science, Research and Innovation Fund via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation, grant B05F650021 (Thailand); the Kavli Foundation; the Nvidia Corporation; the SuperMicro Corporation; the Welch Foundation, contract C-1845; and the Weston Havens Foundation (USA).