Search for low-mass dilepton resonances in Higgs boson decays to four-lepton final states in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for low-mass dilepton resonances in Higgs boson decays is conducted in the four-lepton final state. The decay is assumed to proceed via a pair of beyond the standard model particles, or one such particle and a Z boson. The search uses proton-proton collision data collected with the CMS detector at the CERN LHC, corresponding to an integrated luminosity of 137 fb$^{-1}$, at a center-of-mass energy $\sqrt{s} =$ 13 TeV. No significant deviation from the standard model expectation is observed. Upper limits at 95% confidence level are set on model-independent Higgs boson decay branching fractions. Additionally, limits on dark photon and axion-like particle production, based on two specific models, are reported.


Introduction
Following the discovery of the Higgs boson (H) by the ATLAS and CMS Collaborations [1][2][3] at the CERN LHC, a thorough program of precise measurements [4][5][6] has been carried out to uncover possible deviations from the standard model (SM) or to decipher the nature of the Higgs sector. In particular, various exotic decays of the Higgs boson have been considered, in which small deviations in the Higgs boson decay width or discovery of exotic decay modes could constitute evidence of beyond the SM (BSM) physics. This paper describes a search for exotic decays of the Higgs boson H → ZX or H → XX in the four-lepton (electrons or muons) final state, using a sample of proton-proton collision data at a center-of-mass energy of 13 TeV recorded by the CMS experiment in [2016][2017][2018]. The analyzed data sample corresponds to an integrated luminosity of 137 fb −1 . Here X represents a possible BSM particle that could decay into a pair of opposite-sign, same-flavor (OSSF) leptons. In this paper, we consider two specific BSM models. In both models, leptonic decays of X and Z to either two muons or electrons give rise to the 4 (where 4 may denote 4µ, 2e2µ, or 4e) final states. Assuming narrow-width approximation decays of X, only the mass range m X < m H − m Z ≈ 35 GeV (m X < m H /2 ≈ 62.5 GeV) is kinematically possible for H → ZX (H → XX), where m H and m Z are the Higgs boson mass and Z boson mass, respectively. The decay channel pp → H → 4 has a large signal-to-background ratio. This channel allows a complete reconstruction of the kinematics of the Higgs boson based on final-state decay particles. In this analysis, a mass range of 4.0 < m X < 35.0 GeV (62.5 GeV) is considered.
The first model considered, hereby referred to as the "hidden Abelian Higgs model" (HAHM), concerns theories with a hidden "dark" sector [7][8][9][10][11], with the X particle identified as the dark photon (Z D ), which mediates a dark U(1) D gauge symmetry, which is spontaneously broken by a dark Higgs mechanism. Interactions of the dark sector with SM particles can occur through a hypercharge portal via the kinetic-mixing parameter ε, or through a Higgs portal via the Higgs-mixing parameter κ, as shown in Fig. 1. Details of this theory and subsequent phenomenological implications can be found in Ref. [7]. Several searches for Z D were previously performed by collider experiments, for example ATLAS [12,13] and LHCb [14]. Other experiments, such as beam dump experiments, fixed target experiments, helioscopes, and cold dark matter direct detection experiments, provide complementary sensitivities to Z D . A summary of the experimental coverage of the HAHM model can be found in Refs. [15,16]. o Y F P Q 5 9 G 5 V S W z s m I 0 B e R L z q + 2 L t B D M f A T C t t l Y a G S 3 I P G L N E r C b R S t J Z S W z 5 C u 2 g B G R m h 6 G m n m a i u J r + Z e Q 6 L m + D H Q s j f y Z T n m K u R a a Q F i g v G E M x J h f X + Z H b o 5 k i E + f S j c o W K Z e I r 4 G 7 j T s I 2 v u l c X 3 Y c s t u N k E 0 T x I J 8 A m u m z r q d m 9 i h l i l o K a Z o 3 b g R u 5 f z p a P h R L c r 8 Z 1 0 P r Y W m j D d b G 3 5 L 3 Q 6 t 7 C F / 0 / Z W e l 9 8 C 8 K x c j o s 7 q L 3 b v Q 6 H N X M 7 p 5 u e K M s F b 5 R V 2 E T 3 g w 5 k T 3 r h + r 3 6 H R m P Q 3 A 7 a Q b X Q s h F O j G 1 v s o 4 H z b + 9 o S B F B l w T h p U 6 D 4 N c 9 w 2 W m h K X s V c o y O 2 5 w C M 4 t y b H G a i + q X 7 l E u 1 Y Z Y g S I e 3 D N a r U 2 Q i D M 6 X G W W w 9 M 6 x T t c i c e B s 7 L 3 T y u m 8 o z w s 7 S X J V K C m Y O 6 n u X k B D K o F o N r Y G J p L a X h F J s c R E 2 9 t j r o q 7 I L K 8 d I M J F 8 e w b J x G 7 T B o h + + D 7 c M 3 k x F t e M + 8 5 9 6 u F 3 q v v E P v y D v 2 T j x S + 1 z 7 W v t W + 1 7 / U v 9 R / 1 n / d e W 6 v j a J e e r N r f r v f + T B m U A = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 Z o K 1 q h n 3 k 4 l u 0 v h d v y O K p 3 Z B u w = " > A A A F H n i c h Z T d b t M w F M e z t c A o H + v g k h u L d d J A U Z W k m + B m 0 g R c 7 H J I b J 1 Y q 8 p x T x p r j h 3 Z z q R i 5 Z r X 4 A W 4 h T f g D n E L L 8 B z Y G d d 1 4 9 V W I p 0 8 v + d r x w 7 j n N G l Q 6 C P 2 v r t f q d u / c 2 7 j c e P H z 0 e L O 5 9 e R U i U I S O C G C C X k W Y w W M c j j R V D M 4 y y X g L G b Q j S / e O t 6 9 B K m o 4 B / 0 O I d + h k e c J p R g b a X B V g 3 1 Y h h R b p I s S S i D 0 i Q w 5 i O J 8 7 R s o B l Y S S / L 3 X A / 8 O 3 z o m G p 5 Z Y o L T E d p X o q M E i 0 o Y F P Q 5 9 G 5 V S W z s m I 0 B e R L z q + 2 L t B D M f A T C t t l Y a G S 3 I P G L N E r C b R S t J Z S W z 5 C u 2 g B G R m h 6 G m n m a i u J r + Z e Q 6 L m + D H Q s j f y Z T n m K u R a a Q F i g v G E M x J h f X + Z H b o 5 k i E + f S j c o W K Z e I r 4 G 7 j T s I 2 v u l c X 3 Y c s t u N k E 0 T x I J 8 A m u m z r q d m 9 i h l i l o K a Z o 3 b g R u 5 f z p a P h R L c r 8 Z 1 0 P r Y W m j D d b G 3 5 L 3 Q 6 t 7 C F / 0 / Z W e l 9 8 C 8 K x c j o s 7 q L 3 b v Q 6 H N X M 7 p 5 u e K M s F b 5 R V 2 E T 3 g w 5 k T 3 r h + r 3 6 H R m P Q 3 A 7 a Q b X Q s h F O j G 1 v s o 4 H z b + 9 o S B F B l w T h p U 6 D 4 N c 9 w 2 W m h K X s V c o y O 2 5 w C M 4 t y b H G a i + q X 7 l E u 1 Y Z Y g S I e 3 D N a r U 2 Q i D M 6 X G W W w 9 M 6 x T t c i c e B s 7 L 3 T y u m 8 o z w s 7 S X J V K C m Y O 6 n u X k B D K o F o N r Y G J p L a X h F J s c R E 2 9 t j r o q 7 I L K 8 d I M J F 8 e w b J x G 7 T B o h + + D 7 c M 3 k x F t e M + 8 5 9 6 u F 3 q v v E P v y D v 2 T j x S + 1 z 7 W v t W + 1 7 / U v 9 R / 1 n / d e W 6 v j a J e e r N r f r v f + T B m U A = < / l a t e x i t > s Z D Z D h` < l a t e x i t s h a 1 _ b a s e 6 4 = " V S / X 0 N Z e + K 8 9 w 4 v K W W k e j + + e w H k = " > A A A E 8 X i c d V T N b h M x E N 4 2 A U r 4 a Q t H L h Z J p Y J W 0 e 6 W C j h U q o B D j 0 U i T U U S R d 7 N b N a q 1 1 7 Z 3 k j B 2 g N X e A J u i C t P x A P w B j w A 9 m a T J k 1 j y d L 4 + z 7 P z M 6 M N 8 w o k c r z / m x t 1 + p 3 7 t 7 b u d 9 4 8 P D R 4 9 2 9 / S c X k u c i g k 7 E K R e X I Z Z q v d r X 2 r f a 9 7 q s / 6 j / r P + a S b e 3 q j t P n Z V V / / 0 f j i i L G g = = < / l a t e x i t > Figure 1: Feynman diagrams for Higgs boson decay via the kinetic-mixing (left) or Higgsmixing mechanism (right) [7]. The symbol h represents the Higgs boson, and s represents the dark Higgs boson. The symbol ε represents the kinetic-mixing parameter while κ represents the Higgs-mixing parameter.
The second model involves axion-like particles (ALPs), with X being a pseudoscalar gauge singlet a. Axions were originally proposed to address the strong CP problem [17]. Recently, ALPs were proposed to explain the observed anomaly in the magnetic moment of the muon [18]. Theoretical overviews of the ALP models can be found in Refs. [19,20]. The models are formulated as an effective field theory of ALPs coupled to various SM particles. In particular, the theory allows the coupling between the Higgs boson, Z boson, and the ALP field, or the Higgs boson and the ALP field. These couplings are represented by the Wilson coefficients C ZH /Λ and C aH /Λ 2 , respectively, where Λ is the decoupling energy scale in the effective field theory, or the mass scale of new physics. The former (latter) coefficient gives rise to the exotic decay of H → Za (aa). Various experimental searches for H → aa have been performed [21][22][23][24][25][26]. Recently a direct search for H → Za has been performed targeting a signature with a light and hadronically decaying resonance a with m a < 4 GeV [27]. The present search provides complementary coverage of the phase space of the ALP model with mass greater than 4 GeV. This paper is organized as follows. Section 2 describes the CMS detector and event reconstruction algorithms. Section 3 outlines the collision data used and various software packages used to generate the samples of simulated events. Section 4 summarizes the selection criteria and the categorization of signal events, and Section 5 describes the reducible background estimation method. Section 6 describes the various sources of systematic uncertainties in the search. Finally, results and interpretations are detailed in Section 7, and a summary is given in Section 8. Tabulated results are provided in HEPDATA [28].

The CMS detector and event reconstruction
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 are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), 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. [29].
Events of interest are selected using a two-tiered trigger system [30]. 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 time interval of about 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.
The candidate vertex with the largest value of summed physics-object p 2 T (where p T is the transverse momentum) is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [31, 32] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
The particle-flow (PF) algorithm [33] aims to reconstruct and identify each individual particle in an event (PF candidate), with an optimized combination of information from the various elements of the CMS detector. The energy of photons is obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.
The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [34]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.
Muons in the four lepton final state 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 a relative p T resolution, for muons with p T up to 100 GeV, of 1% in the barrel and 3% in the endcaps [35].
Electrons in the four lepton final state with p T > 7 GeV and |η| < 2.5 are identified by a multivariate discriminant, which is constructed by observables related to the bremsstrahlung along the electron trajectory, ECAL energy measurements, electromagnetic showers, missing pixel detector hits, and the photon conversion vertex fit probability [36]. The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7 to 4.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. The dielectron mass resolution for Z → ee decays when both electrons are in the ECAL barrel (endcap) is 1.9% (2.9%).
This analysis focuses on promptly produced signal processes. To reduce the contributions from leptons arising from hadron decays within jets, a requirement is imposed on each lepton candidate using a variable defined as: where the sums are over the PF candidates within a cone of radius R = √ ∆η 2 + ∆φ 2 < 0.3 (where φ is the azimuthal angle in radians), p i T represents transverse momenta from each particle i, where i represents either charged hadrons, neutral hadrons, photons, or particles originating from overlapping proton-proton interactions (pileup) [37]. For muons, the isolation is required to be I µ < 0.35. For electrons, this variable is included in the multivariate discriminant for datasets in 2017 and 2018, while for the dataset in 2016, an isolation requirement I e < 0.35 is imposed on each electron candidate. In addition, the three-dimensional impact parameter of electrons and muons is required to be consistent with the primary collision vertex. The requirement implies a negligible acceptance to signal models with long-lived X.
An algorithm is utilized to correct for effects arising from final-state radiation (FSR) from leptons. PF-reconstructed photons are considered as FSR candidates if they satisfy the requirement p γ T > 2 GeV and I γ < 1.8, where I γ is calculated similarly to the lepton isolation variable. Then each FSR candidate is assigned to the closest lepton in the event. The candidates are further required to have ∆R(γ, )/(p γ T ) 2 < 0.012 GeV −2 and ∆R(γ, ) < 0.5. These candidates are excluded from the calculation of the lepton isolation variables.
Lepton reconstruction and selection efficiencies are measured in data by a "tag-and-probe" technique with an inclusive sample of Z boson events [38]. The difference between the efficiencies in data and simulation are observed to be around 1-4%, depending on p T and η of the lepton considered. The differences are used to correct lepton efficiencies in simulation.  [7] at leading order. Cross sections for each Z D signal are calculated by multiplying the next-to-next-to-next-to-leading order (NNNLO) Higgs production cross section [42] by the branching fraction of H → ZZ D and H → Z D Z D , respectively [7]. Final states with τ leptons are neglected as their contribution to the signal region yield is below 1%. Signal contributions from vector-boson fusion and associated production with a top quark pair or a vector boson are also omitted.

Data and simulated samples
The SM Higgs boson simulation samples, which include gluon fusion, vector boson fusion, and associated production with a top quark pair or a vector boson, and the simulated ZZ background from quark-antiquark annihilation are generated at next-to-leading order (NLO) in perturbative quantum chromodynamics with POWHEG v2 [43][44][45][46]. The cross section for the dominant production mode, gluon fusion, is taken at NNNLO [42].
Decays of the Higgs boson to four leptons are simulated with JHUGEN 7.0.2 [47,48]. The nonresonant process of gg → ZZ process is simulated at LO with MCFM 7.0.1 [49]. NLO correction factors [50] are applied to the gg → ZZ process.
Minor backgrounds from ttZ and triboson production processes are also simulated at LO and NLO, respectively, with the MADGRAPH5 aMC@NLO 2.  [54,55]. The response of the CMS detector is modeled using the GEANT4 program [56,57]. Simulated events are reweighted according to a specified instantaneous luminosity and an average number of pileup events. par

Event selection
In the trigger system, events are required to have more than two leptons. The overall trigger efficiency is measured in data using a sample of 4 events from single-lepton triggers and agreements are observed with simulation within 5%, and is found to be larger than 99%.
A set of requirements is applied to maximize the sensitivity of the search for a potential signal in the ZX and XX event topologies. In both cases, at least four well-identified and isolated leptons from the primary vertex are required, possibly accompanied by an FSR photon. Each muon (electron) is required to have p T > 5 GeV (7 GeV). All four leptons must be separated from each other by ∆R( i , j ) > 0.02. The leading (subleading) lepton p T is required to satisfy p T > 20 GeV (10 GeV). The four-lepton invariant mass m 4 is required to be within 118 < m 4 < 130 GeV.
To further suppress background contributions from hadron decays in jet fragmentation or from the decay of low-mass resonances, all opposite-charge leptons pairs, regardless of lepton flavor, are required to satisfy m + − > 4 GeV.
For each event in the ZX and XX searches, dilepton pair candidates are formed by considering all OSSF leptons. The dilepton invariant mass m + − for each candidate is required to be within 4 < m + − < 120 GeV, however the mass window around the Υbb bound states (8.0 < m Υ < 11.5 GeV) is also excluded.
Two dilepton candidates are then paired to form a ZX or XX event candidate. For the ZX search, Z 1 is the OSSF dilepton pair with an invariant mass closest to the Z boson mass [58] (representing Z in ZX), and Z 2 is the other pair (X). For the XX search, Z 1 is the OSSF dilepton pair with the larger invariant mass, and Z 2 is the lower-mass pair. For the ZX search, m Z 1 is required to be larger than 40 GeV. For the XX search, m Z 1 and m Z 2 must lie between 4 and 62.5 GeV. For events with more than four selected leptons, the combination of four leptons with m Z 1 closest to the Z boson is used for the ZX candidate, while the combination with the least value of (m Z 1 − m Z 2 )/(m Z 1 + m Z 2 ) is used to select XX candidates with similar invariant masses.
Four final-state lepton categories can be defined as 4µ, 2µ2e, 4e, 2e2µ, where the order of lepton flavors corresponds to Z1 and Z2 flavors. For the 4µ and 4e final states, one alternative pairing of the four leptons is possible, labelled by Z a and Z b . For the ZX search, events with m Z b < 12 GeV and m Z a closer to the Z boson mass than Z 1 are discarded to suppress background contributions from on-shell Z and low-mass dilepton resonances. For the XX search, the XX candidate with the smallest value of (m Z 1 − m Z 2 )/(m Z 1 + m Z 2 ) is chosen.

Irreducible background estimation
Irreducible backgrounds for this search come from processes including a SM Higgs boson, as well as nonresonant production of ZZ via quark-antiquark annihilation or gluon fusion, and rare backgrounds such as tt + Z and triboson production. These backgrounds are estimated from simulation. Details of the simulation used for each of the backgrounds are described in Section 3.

Reducible background estimation
The reducible backgrounds in the 4 final state can arise from the leptonic decays of heavyflavor hadrons, in-flight decays of light mesons within jets, charged hadrons misidentified as electrons when in proximity of a π 0 , and photon conversions. These backgrounds primarily arise from the Z + jets process. Additional physics processes with kinematics similar to the signal include tt, Zγ, and WZ.
Two dedicated control regions are used to estimate the contribution from these backgrounds. The first (second) control region consists of events with two (three) leptons passing the lepton identification and isolation requirements and two (one) leptons failing the requirements, and is denoted as the 2P2F (3P1F) region. Backgrounds with only two prompt leptons, such as Z + jets and tt, are estimated by the 2P2F region, while backgrounds with three prompt leptons, such as WZ and Zγ with the photon converting to an electron pair, are estimated by the 3P1F region. Other than the lepton requirements, the 3P1F and 2P2F regions follow the same event selection and alternative pairing algorithms as in the signal region to closely mimic its kinematics.
The lepton misidentification rates f µ and f e are measured as a function of lepton p T and η with a sample which includes a Z candidate, formed by a pair of leptons passing the selection requirement of the analysis, and an additional lepton passing a relaxed requirement. These rates are measured separately in the data samples from 2016, 2017, and 2018. In addition, the mass of the Z candidate is required to satisfy the condition |m Z 1 − m Z | < 7 GeV to reduce contributions from WZ and tt processes, and p miss T is required to be less than 25 GeV.
To estimate the background contribution in the signal region, events in the 3P1F and 2P2F control regions are reweighted by lepton misidentification probabilities. Each event i in the 3P1F region is weighted by a factor f i 4 /(1 − f i 4 ), where f i 4 corresponds to the lepton misidentification rate of the failed lepton in the event. Physics processes in the 2P2F control region can contribute to the 3P1F region and are estimated by reweighting 2P2F events with f i , where f i 3 and f i 4 correspond to the lepton misidentification rates of the two failed leptons in the event. A minor contribution from ZZ events to the 3P1F control region is estimated from simulation and subtracted. The expected yield for the signal region can then be estimated as: where each sum is over all 3P1F and 2P2F events, respectively.
Furthermore, dedicated validation regions, which include adjacent m 4 regions to the signal region (70 < m 4 < 118 GeV, 130 < m 4 < 200 GeV), are defined to inspect the level of agreement between data and predictions.

Systematic uncertainties
Experimental sources of the systematic uncertainties applicable to all final states include the integrated luminosity uncertainty and the lepton identification and reconstruction efficiency uncertainty. The integrated luminosities of the 2016, 2017, and 2018 data-taking periods are individually known with uncertainties in the 1.2-2.5% range [59][60][61], while the total Run 2 (2016-2018) integrated luminosity has an uncertainty of 1.6% [62], the improvement in precision reflecting the (uncorrelated) time evolution of some systematic effects. Lepton efficiency uncertainties are estimated in bins of lepton p T and η using the tag-and-probe method, as described in Section 2. These uncertainties on each lepton candidate lead to variations from 2.5 to 16.1% on event yields, dependent on final-state lepton categories. In addition, the systematic uncertainties in the lepton energy scale are determined by fitting the Z → mass distribution in bins of lepton p T and η with a Breit-Wigner parameterization convolved with a doublesided Crystal Ball function [63]. Systematic uncertainties in the estimation of the reducible background are derived from the level of agreement between data and predictions in the validation regions in each lepton category (23-48% depending on data taking period), arising from different background compositions between signal and control regions (30-38% depending on lepton category), and from misidentification rate uncertainties (35-100% depending on lepton category).
Theoretical uncertainties that affect both the signal and background estimation include uncertainties in the renormalization and factorization scales and the choice of the PDF set. The uncertainty from the renormalization and factorization scales is determined by varying these scales between 0.5 and 2 times their nominal value while keeping their ratio between 0.5 and 2. The uncertainty from the PDF set is determined by taking the root-mean-square of the variation when using different replicas of the default NNPDF set [64]. An additional uncertainty of 10% in the K factor used for the gg → 4 prediction is included [37]. To estimate the effect of the interference between the signal and background processes, three types of samples are generated using the MADGRAPH5 aMC@NLO 2.4.2 [39][40][41] generator: inclusive sample (H → ZZ * → 4 , H → ZX/XX → 4 ), signal-only sample H → ZX/XX → 4 and background-only sample H → ZZ * → 4 . The inclusive sample contains background, signal, and interference contributions. The effect of the interference on the normalization of the signal is estimated by taking the difference of the inclusive sample cross section and the sum of the cross sections of the signal and background samples. This difference is at 1-2% after the final event selection. Theoretical values of branching fractions B(Z D → ee or µµ) are calculated in Ref. [7]. The calculations are based on experimental measurements of the ratio of the hadronic cross section to the muon cross section in electron-positron collisions R µµ /R had up to m Z D = 12 GeV and a next-to-leading order theoretical calculation for m Z D > 12 GeV. To account for uncertainties in these theoretical estimates, a conservative 20% (10%) uncertainty is assigned to them for m Z D < 12 GeV (m Z D > 12 GeV) [7]. Differences in the kinematic properties between the HAHM and ALP model have been inspected. For the determination of model-independent exclusion limits, differences in acceptances are included as systematic uncertainties, ranging from 10% (m X ∼ 4 GeV) to 30% (m X ∼ 35 GeV for ZX, m X ∼ 60 GeV for XX), while they are used to correct signal yields for the determination of ALP exclusion limits.
In the combination of the three data taking periods, the theoretical uncertainties and experimental ones related to leptons are correlated across all data taking periods, while all others from experimental sources are taken as uncorrelated. The sensitivity of this analysis is dominated by data statistical uncertainty rather than systematic uncertainties.

Results and interpretation
Dilepton mass distributions for the ZX and XX selections are shown in Figs. 2 and 3, respectively. The dilepton mass variable for the XX selection shown in Fig. 3 is m Z12 = (m Z 1 + m Z 2 )/2, which should peak at m X in case of a signal H → XX. In all cases, the observed distributions agree well with standard model expectations within the assigned uncertainties.
These results are further interpreted as upper limits on model-independent branching fractions and model parameters for the dark photon and ALP models. For interpretations of the results of the ZX selection, 351 mass hypotheses are considered. Each mass hypothesis m i is defined with an incremental step of 0.5%, as m i = 4.20 × 1.005 i , where i = 0, 1, 2, . . . , 424, excluding the mass points around the Υ bb bound states between 8.0 < m Υ < 11.5 GeV (i = 130, 131, . . . , 201). The incremental step is chosen so as not to miss any potential signal contribution due to detector resolution in m Z 2 . For each mass hypothesis, the counting experiments are performed on the m Z 2 distribution, with the bin centered at each mass hypothesis. Because of the finite mass resolution of m Z 2 , the choice of the bin width needs to be defined such that most of the signal contribution is included in the bin. The bin width is defined as 0.04 (0.10) × m i for the 4µ and 2e2µ (4e and 2µ2e) categories. This width is chosen as two times the m Z 2 resolution and includes ≈95% of signal events. The normalization of the Higgs background is allowed to float freely in the likelihood fit. For each mass hypothesis, events outside the mass window are included as a sideband to constrain the normalization parameter. No significant deviation with respect to the SM prediction is observed.
For interpretations of the results of the XX selection, 462 mass hypotheses are considered instead. In contrast to the ZX interpretations, the counting experiments are performed by con- L m,SB = ∏ Pois(n |µ Higgs n Higgs, where the function Pois(n|x) is the Poisson probability to observe n events, when the expectation is x. The symbol m represents a particular mass hypothesis. The likelihood term L m,SR (L m,SB ) corresponds to the event yields within (outside) the mass window. The symbol µ is the signal strength parameter, µ Higgs represents the free floating normalizing parameter on the SM Higgs boson process, represents each lepton category, b represents each background process, s represents a particular signal process and n i,m, represents the yield in a mass window associated with the mass hypothesis m, from a source i and the lepton category . In Equation 5, the symbols n Higgs, and n b, represent the yields of the SM Higgs boson and other backgrounds b outside the mass window for the lepton category . Systematic uncertainties are included and profiled as nuisance parameters ρ [65].
For each interpretation, 95% exclusion limits are obtained with an asymptotic formulation of the modified frequentist CL s criterion as described in Refs. [65][66][67][68] with the ZX selection and full CL s approach for the XX selection.  Figure 3: Event yields against m Z12 = (m Z 1 + m Z 2 )/2 with the XX selection for the 4µ, 2e2µ and 4e final states. Numbers in the legend show the total event yields with the XX selection corresponding to data, and the expected yields for each background and signal processes, along with the corresponding statistical uncertainty coming from the amount of simulated data.

Model-independent limits
Upper limits at 95% confidence level (CL) are derived on model-independent branching fractions with the ZX and XX selections assuming three decay channels: a flavor symmetric decay of X to a muon or an electron pair, exclusive X decays to a muon pair, and exclusive X decays to an electron pair. Acceptance effects arising from different signal models are included as systematic uncertainties in the signal yields after event selection. Little model dependence is expected as the event selection is defined without using angular information between the lep-tons. Figures 4 and 5 show the exclusion limits on the model-independent branching fractions with the ZX and XX selections, respectively. The weaker observed limit in the XX selection at m X ≈ 18 GeV is due to one observed data event and does not represent a significant statistical deviation from the background hypothesis. Kinematic differences between the dark photon and ALP models are included as systematic uncertainties, as detailed in Section 6. → ee) assuming X decays to dielectrons only, and B(H → ZX)B(X → ee or µµ) assuming a flavor symmetric decay of X to dimuons and dielectrons. The dashed black curve is the expected upper limit, with one and two standard-deviation bands shown in green and yellow, respectively. The solid black curve is the observed upper limit. The red curve represents the theoretical cross section for the signal process H → ZX → 4 . The discontinuity at 12 GeV in the uncertainty is due to the switch from experimental to theoretical uncertainty estimates of B(Z D → ee or µµ), as described in Ref. [7]. The symbol ε is the kinetic-mixing parameter. The grey band corresponds to the excluded region around the bb bound states of Υ.

Limits on dark photon model parameters
Upper limits at 95% CL are obtained on the Higgs-mixing parameter κ and B(H → Z D Z D ) with the XX selection, as shown in Fig. 6, assuming κ ε. The LHC provides unique sensitivity to the parameter κ due to the presence of the Higgs boson. In addition, this analysis provides some sensitivity to ε, but the upper limits are almost an order of magnitude weaker than those from the Drell-Yan search and from the LHCb Collaboration [14], and hence are not reported in this paper.  2 assuming X decays to dimuons only, B(H → XX)B(X → ee) 2 assuming X decays to dielectrons only, and B(H → XX)B(X → ee or µµ) 2 assuming a flavor symmetric decay of X to dimuons and dielectrons. The dashed black curve is the expected upper limit, with one and two standarddeviation bands shown in green and yellow, respectively. The solid black curve is the observed upper limit. The red curve represents the theoretical cross section for the signal process H → XX → 4 . The discontinuity at 12 GeV in uncertainty is due to the switch from experimental to theoretical uncertainty estimates of B(Z D → ee or µµ), as described in Ref. [7]. The symbol κ is the Higgs-mixing parameter. The grey band corresponds to the excluded region around the bb bound states of Υ.

Limits on the ALP model
Upper limits at 95% CL are calculated on the Wilson coefficients C ZH /Λ and C aH /Λ 2 , as shown in Fig. 7, where C ZH is the effective coupling parameter of the Higgs boson, Z boson, and the ALP, C aH is the effective coupling parameter of the Higgs boson and the ALP, and Λ is the new physics scale. In both interpretations, the ALP is assumed to decay promptly with B(a → ee or µµ) = 1, with equal fractions to muons and electrons. The last six mass hypotheses are omitted in the calculation of upper limits on C ZH /Λ to match the m a range adopted in Ref. [20]. Kinematic differences between the dark photon and ALP models are included as corrections on signal region yields, as detailed in Section 6.   No significant deviations from the standard model expectations are observed. The search imposes experimental constraints on products of model-independent branching fractions of B(H → ZX), B(H → XX) and B(X → ee or µµ), assuming flavor-symmetric decays of X to dimuons and dielectrons, exclusive decays of X to dimuons, and exclusive decays of X to dielectrons, for m X > 4 GeV. In addition, two well-motivated theoretical frameworks beyond the standard model are considered. Due to the presence of the Higgs boson production in LHC proton-proton collisions, the search provides unique constraints on the Higgs-mixing parameter κ < 4 × 10 −4 at 95% confidence level (CL) in a dark photon model with the XX selection, in Higgs-mixing-dominated scenarios, while searches for Z D in Drell-Yan processes [14,69] provide better exclusion limits on ε in kinetic-mixing-dominated scenarios. For the axion-like particle model, upper limits at 95% CL are placed on two relevant Wilson coefficients C ZH /Λ and C aH /Λ 2 . This is the first direct limit on decays of the observed Higgs boson to axion-like particles decaying to leptons.

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 centres and personnel of the Worldwide LHC Computing Grid and other centres 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, the CMS detector, and the supporting computing infrastructure provided by the follow- [25] ATLAS Collaboration, "Search for heavy particles in the b-tagged dijet mass distribution with additional b-tagged jets in proton-proton collisions at √ s = 13 TeV with the ATLAS experiment", Phys. Rev. D 105 (2022) 012006, doi:10.1103/PhysRevD.105.012001, arXiv:2110.00313.
[26] ATLAS Collaboration, "Search for Higgs bosons decaying into new spin-0 or spin-1 particles in four-lepton final states with the ATLAS detector with 139 fb −1 of pp collision data at √ s = 13 TeV", 2021. arXiv:2110.13673. Submitted to JHEP.