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 $\sqrt{s}=13$ TeV

Searches are conducted for new spin-0 or spin-1 bosons using events where a Higgs boson with mass $125$ GeV decays into four leptons ($\ell =$ $e$,$\mu$). This decay is presumed to occur via an intermediate state which contains two on-shell, promptly decaying bosons: $H \rightarrow XX/ZX \rightarrow 4\ell$, where the new boson $X$ has a mass between 1 and 60 GeV. The search uses $pp$ collision data collected with the ATLAS detector at the LHC with an integrated luminosity of 139 fb$^{-1}$ at a centre-of-mass energy $\sqrt{s}=13$ TeV. The data are found to be consistent with Standard Model expectations. Limits are set on fiducial cross sections and on the branching ratio of the Higgs boson to decay into $XX/ZX$, improving those from previous publications by a factor between two and four. Limits are also set on mixing parameters relevant in extensions of the Standard Model containing a dark sector where $X$ is interpreted to be a dark boson.


Introduction
Although the Higgs boson was discovered at the Large Hadron Collider (LHC) in 2012 [1,2], there is good reason to believe that the description of the Higgs sector of the Standard Model (SM) is still incomplete. Besides the well-known issues of naturalness and baryon asymmetry, astrophysical observations implying the existence of dark matter motivate extensions to the Higgs sector of the SM, particularly those that propose the existence of a 'dark' (i.e., hidden) sector, with its own hidden-sector particles [3,4].
An attractive way to search for new physics in the Higgs sector is through non-standard ('exotic') decays of the Higgs boson. Existing precision measurements of the properties of the Higgs boson still allow a branching ratio of up to about 30% to non-standard decays (assuming that the couplings of the Higgs boson to the and bosons are not larger than their SM values) [5][6][7]. Further, since the SM predicts a very narrow decay width for the Higgs boson, even a small coupling to a new light state could result in a significant branching ratio to that state. In addition, new hidden-sector particles may preferentially couple to the Higgs boson, making it a 'portal' to explore this new physics [8][9][10][11]. Such exotic decays of the Higgs boson are predicted by many proposed extensions to the SM, including models with a first-order electroweak phase transition [12,13], models with neutral naturalness [14][15][16], and models with a hidden sector [17][18][19][20][21][22][23][24][25][26][27], as well as by several models of dark matter [28][29][30][31][32][33], including some posited to explain observed excesses of astrophysical positrons [34][35][36]. They are also predicted by the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [37][38][39][40][41][42].
This paper reports three related searches, each of which looks for a SM Higgs boson decaying via a new boson into a final state consisting of four charged leptons (ℓ ≡ , ). All use the full LHC Run 2 data set of about 139 fb −1 that the ATLAS detector collected from proton-proton collisions at a centre-of-mass energy of √ = 13 TeV. Following the models motivating these analyses, the new boson could be either a dark-sector vector boson or a scalar boson, denoted by . The three searches considered are: • High-mass (HM): → → 4ℓ (15 GeV < < 60 GeV).
The LM analysis uses only the 4 final state because the selection efficiency for isolated muons is significantly larger than that for isolated electrons in this mass range (see Section 5.4). These searches are sensitive to any intermediate bosons within the considered mass ranges that are narrow, on-shell, and decay promptly. This paper provides model-independent fiducial cross-section limits, as well as limits based on the specific models described in Section 2.
This work extends previous searches performed by ATLAS with 20 fb −1 of data collected at √ = 8 TeV [43] and with 36 fb −1 of data collected at √ = 13 TeV [44]. In addition to a larger data sample and improved lepton identification, the signal region selection of the HM analysis has been re-optimized. Other similar searches, including searches for pairs of light bosons decaying into muons, -leptons, photons, and/or jets, as well as searches for a single light boson decaying into a pair of muons, using both √ = 8 TeV and 13 TeV data, have been performed by ATLAS [45][46][47][48][49], CMS [50][51][52][53], and LHCb [54]. Further searches for a SM Higgs boson decaying into undetected particles are reported in Refs. [55,56]. This section is followed by a summary of the theoretical models used in the interpretation of the results (Section 2). Next, the detector is described (Section 3), followed by discussions of features that are common to all three analyses, including the samples of data and simulated events (Section 4), the reconstruction of lepton candidates and of their combinations (Sections 5.1 and 5.2), the event selections (Section 5.3), and the common systematic uncertainties (Section 6). Next, aspects specific to each analysis are described (Sections 7 to 9). Finally, the ways the analyses are combined to extract limits, and the interpretations of the results in terms of the theoretical models, are presented in Section 10, and a summary is given in Section 11.
The coupling strength of the d boson to SM particles, and hence its lifetime (assuming no significant decays to non-SM particles), is determined by the mixing parameter . The decays of the d boson, on the other hand, are determined by the gauge couplings, and the decay branching ratios are largely independent of for 1. Over the d mass range 1 GeV < d < 60 GeV, the branching ratio for decays into electron or muon pairs can be 10%-15% [21]. Over the same mass range, the decay is prompt for 10 −5 [21]. For smaller values of , the decay vertex would be significantly displaced from the interaction point, while for 10 −8 the lifetime of the d boson becomes long enough for it to likely escape the detector. Also, the decay width of the d boson is very small ( 1 GeV) for 1 and d < 60 GeV. ATLAS and CMS have searched for these long-lifetime signatures in collisions at energies of both 8 TeV [61-64] and 13 TeV [65][66][67][68][69].
If the U(1) d symmetry is broken by an additional dark Higgs boson , then there could be mixing with strength between the SM Higgs boson and the dark Higgs boson [21][22][23][24][25][26]. The observed Higgs boson would be one of the mass eigenstates and could also decay into dark-sector particles, including dark Higgs bosons that subsequentially decay into SM fermions. The dark Higgs boson would inherit the Yukawa couplings from the SM Higgs boson and decay preferentially into high-mass fermion pairs. A further possibility is mass mixing between the d boson and the SM boson [23,24]. If the mass term for this mixing is written as 2 d , with = d / , then is the model parameter describing the mixing.
The processes probed in this paper that involve a SM Higgs boson decaying into d bosons are depicted in Figures 1(a) and 1(b) and are included in the Hidden Abelian Higgs Model (HAHM) [21]. The decay → d is sensitive to the parameters and d , but does not depend on . However, the presence of an irreducible background from the SM → * process means that this signal can be observed only as a peak in the dilepton mass spectrum over the background. The process → d d , in contrast, is much more easily separated from SM backgrounds and hence is potentially sensitive to smaller values of the kinetic mixing , where it is only required that the mixing be large enough for the d boson to decay promptly. However, this process does require mixing between the SM and dark-sector Higgs bosons and thus depends on .
Limits on the kinetic mixing of 0.03 have been set from precision electroweak measurements [21,70,71] over the range 1 GeV < d < 200 GeV. Searches for dilepton resonances, → d → ℓℓ, at the LHC for d < imply that 0.005-0.020 for 20 GeV < d < 80 GeV [72]. Other searches rule out 10 −3 for 10 MeV < d < 10 GeV [73][74][75][76][77][78]. The → → 4ℓ analyses constrain the Higgs mixing parameter , while the → d → 4ℓ analysis provides information about the kinetic mixing parameter .  [21] (to which the HM and LM analyses are sensitive). The d gauge boson decays into SM particles through kinetic mixing with the hypercharge field (with branching ratios that are nearly independent of ). The d vertex factor is proportional to whereas the d d vertex factor is proportional to . (c) illustrates the decay of a Higgs boson into dark Higgs scalars or pseudoscalars that couple to SM particles through mixing with the SM Higgs field in models with an extended Higgs sector (Section 2.2).

Extended Higgs sectors
Models containing two Higgs doublets and an additional scalar field (2HDM+S) [22,79] are also relevant for the search for → → 4 . Two-Higgs-doublet models (2HDMs) generically contain two neutral scalars 1,2 , two charged scalars ± , and one neutral pseudoscalar . The lighter of the neutral scalars 1 is identified as the observed Higgs boson , while the other states are constrained to be heavy by existing data [80,81]. Adding a complex scalar singlet that mixes weakly with 1,2 gives two additional states, a scalar and a pseudoscalar . If these are lighter than /2, then → and → decays are allowed (Figure 1(c)). This paper probes the process → → 4 , but limits on → → 4 also apply to → → 4 .
The decays of the scalar and pseudoscalar into fermions are determined by their Yukawa couplings [22], implying that the branching ratio to electrons is very small, and that the branching ratio to muons is smaller than that of the d vector bosons described previously. Branching ratios for → and → can be significant in the range 2 < < 2 , ranging from 10 −2 to 10 −1 in some regions of parameter space [22]. In 2HDMs, there are several possible ways for the Higgs sector to couple to fermions. Of these, type-III models (in which leptons and quarks couple to different Higgs doublets) at large tan (where tan is the ratio of the vacuum expectation values of the two Higgs doublets) are particularly interesting for these analyses. A light pseudoscalar can correspond to the -symmetry limit of the NMSSM [82,83], which reduces the need for fine-tuning and addresses the -problem [84]. Searches for exotic decays of the

Data and simulated event samples
The results in this paper are based on 139 fb −1 of √ = 13 TeV proton-proton ( ) collision data collected with the ATLAS detector over the period 2015-2018.
Monte Carlo (MC) simulation is used to determine expected contributions from both the signal processes and most background processes. For most samples, detector effects were included using a G 4 [97] simulation of the ATLAS detector [98]. The → , → d , and → signal samples, as well as a portion of the → * and triboson background samples, instead used a fast simulation [98] which relies on a parameterization of the calorimeter response [99]. The effects of pile-up (additional collisions in the same or a neighbouring bunch crossing) are included in the simulation. Weights are applied to the simulated events to correct for small differences between data and simulation in the reconstruction, identification, isolation, and impact parameter efficiencies for electrons and muons [100-102]. Further, the lepton momentum scales and resolutions in the simulation are adjusted to match the data [100, 102, 103]. Signal samples involving a d vector boson were generated according to the HAHM [21,22,25,26] implementation in M G 5_aMC@NLO, with the Higgs bosons being produced via gluon-gluon fusion (ggF) and the Higgs boson mass set to = 125 GeV. For the → d d process, and were both set to 10 −4 and samples were generated with d = 0.5 GeV, 1 GeV, 2 GeV, and every 5 GeV in the range 5 GeV ≤ d ≤ 60 GeV. For the → d process, was changed to 10 −10 , and samples were generated every 5 GeV in the range 15 GeV ≤ d ≤ 55 GeV. Final states with -leptons were not included; the change in signal region yield due to the omission of these decays was below 1% and thus neglected. The much smaller production of signal events by vector-boson fusion (VBF), , and¯was also omitted.
Prompt-lepton backgrounds are estimated directly from MC simulations. They arise primarily from the SM → * → 4ℓ process along with the non-resonant * → 4ℓ process. Smaller leptonic backgrounds arise from triboson production as well as¯+ decays. Decays involving → were found to contribute negligibly to the background yields and are thus not included in the simulation. Backgrounds with jets misidentified as leptons are estimated with data-driven methods, detailed below in the individual analysis sections.
The non-resonant¯→ * → 4ℓ background process [149] was simulated using S 2.2.2 at NLO for up to one additional parton and at LO for up to three additional partons. Matrix element calculations were matched and merged with the S parton shower based on the Catani-Seymour dipole factorization [150,151] using the MEPS@NLO prescription [152][153][154][155]. The virtual QCD corrections are provided by the O L library [156][157][158]. The gluon-initiated process ( → * → 4ℓ) was simulated in the same manner, except that it was at LO, and the s-channel diagrams were omitted to avoid double-counting. The gluon-initiated process has a large QCD correction at NLO, so the cross section was scaled by a NLO/LO -factor of 1.70 ± 0.15 [159]. Interference between the → * → 4ℓ and → → * → 4ℓ processes is neglected.
Higher-order electroweak processes include triboson production ( ) and vector-boson scattering (VBS). These processes can yield final states with four leptons along with two additional particles. They were generated with S 2.2.2 at NLO for the inclusive processes and at LO for up to two additional partons, with the same treatment as for¯→ * → 4ℓ. Higgs boson production via VBF was subtracted from these samples in order to avoid double-counting.
The process¯+ ( → ℓℓ) was generated with S 2.2.0 at LO with up to one additional parton emission.
Other, reducible, backgrounds have fewer than four prompt leptons in the final state, but can be accepted by the signal selection if there are additional leptons from heavy-flavour decay or jets misidentified as leptons. The + jets process was generated with S 2.2.1 using NLO matrix elements for up to two partons and LO matrix elements for up to four partons. The¯process was generated with P B at NLO with the ℎ damp parameter, which regulates the high-T radiation against which the¯system recoils, set to 1.5 top = 258 GeV [160]. The process was also generated with P B at NLO with the CT10 PDF.

Lepton reconstruction
For the analyses considered in this paper, the final-state objects of interest are electrons and muons.
Electrons are reconstructed and identified from charged-particle tracks in the ID that match energy deposits in the calorimeters [ To avoid identifying the same detector signature as multiple particles, an electron candidate that has the same ID track as a muon candidate is ignored, unless the muon is only calorimeter-tagged, in which case the muon is ignored instead. Electrons that have the same track or cluster as a higher-T electron are also ignored.

Invariant kinematic mass variables
All three analyses considered in this paper involve looking for mass resonances in final states consisting of a quadruplet of two same-flavour opposite-sign (SFOS) lepton pairs: ( + − + + − ), ( + − + + − ), or ( + − + + − ). The invariant masses of the two pairs are denoted by 12 and 34 , where 12 is taken to be the one closer in mass to the boson, If all four leptons have the same flavour, then for a given 12 and 34 labelling, alternative SFOS pairings can also be defined. An invariant mass 14 is constructed from the positively charged lepton of the 12 pair and the negatively charged lepton of the 34 pair. The other alternative pairing 23 is constructed analogously.

Common event selection
The analyses all involve a Higgs boson decaying into a pair of new bosons , or into a new boson along with a boson, which in turn decay into pairs of leptons. The bosons are presumed to be on-shell, so the strategy is to search for resonances in the relevant dilepton mass distributions. Each analysis defines a signal region (SR) via a series of selections on measured quantities which maximizes the sensitivity to the signal.
All three analyses share a common preselection, but differ in the subsequent steps of selecting the candidate final-state leptons, forming them into quadruplets, selecting one of those quadruplets, and applying further requirements to the selected quadruplet. Table 2 shows the event selections of the different analyses.
The common preselection requires that events were recorded with the detector in good operating condition [162] and without excess calorimeter noise [163]. Each event must have an identified primary vertex with at least two tracks [164] and at least four lepton candidates. Events were triggered by requiring either one or two lepton candidates, where the candidates could be either electrons or muons [165][166][167]. The lepton candidates identified offline must match candidates identified by the trigger. The trigger T requirements range from T > 7 GeV to T > 60 GeV, depending on lepton multiplicity and flavour. In either case, the trigger efficiency is above 95% (relative to signal region events surviving all other event selections).
Electron and muon candidates are reconstructed as described in Section 5.1. Electrons must be within the central region of the detector (| | < 2.47), have T > 7 GeV, have a longitudinal impact parameter 0 that satisfies | 0 sin | < 0.5 mm with respect to the primary vertex, and have an additional associated hit in the insertable B-layer. Muons must be within the acceptance of the muon spectrometer, | | < 2.7. All muons must have T > 5 GeV, while CT muons must pass the stronger requirement T > 15 GeV. Lastly, all muon candidates that are associated with a vertex, i.e. all except SA muons, must have a longitudinal impact parameter with respect to the reconstructed primary vertex satisfying | 0 sin | < 0.5 mm and a transverse impact parameter with respect to the position of the beam satisfying 0 < 1 mm.
All possible quadruplets (Section 5.2) are formed from the selected leptons. A quadruplet may contain no more than one SA or CT muon, and at least one lepton in the quadruplet must correspond to a lepton found by one of the triggers satisfied by the event. The three highest-T leptons must satisfy, respectively, T > 20 GeV, T > 15 GeV, and T > 10 GeV. Except for the LM analysis, for which the angular separation between leptons can be very small, all pairs of same-flavour leptons in the quadruplet must satisfy Δ (ℓ, ℓ ) > 0.1, while different-flavour pairs must satisfy Δ (ℓ, ℓ ) > 0.2. At least one quadruplet per event is required. For the HM and LM analyses, if there is more than one quadruplet passing these requirements, the one with the smallest mass difference between the two pairs, Δ ℓℓ = | 12 − 34 |, is chosen. The analogous procedure for the ZX analysis is described in Section 5.6.
The leptons in the quadruplet must be isolated from other deposits in the calorimeter or ID tracks [101,102]. This rejects backgrounds in which leptons arise from the decay of heavy-flavour jets, or in which hadronic jets are misidentified as leptons. For each lepton, the sum of the transverse energies of topological clusters [168] within a cone of Δ = 0.2 around it (excluding energy attributed to the lepton itself) must be less than 20% of its T for electrons, and less than 30% of its T for muons. The transverse momenta of tracks in a cone around the lepton are also summed, and must be less than 15% of its T . The -radius of the cone depends on the momentum of the lepton. For electrons, the radius is Δ = min(0.2, 10 GeV/ T ), while for muons, it is Δ = min(0.3, 10 GeV/ T ). In both cases, tracks and energy clusters attributed to other leptons in the quadruplet are also excluded from the sums. This is particularly important for the LM analysis (Section 5.5), where the angular separation between leptons may be very small.
In addition to the impact parameter requirements discussed earlier, each lepton in the quadruplet must have transverse impact parameter significance 0 / 0 < 5 for electrons and 0 / 0 < 3 for muons (with the exception of SA muons, which do not have an ID track), where 0 is the estimated error in the reconstructed transverse impact parameter 0 .

HM event selection
The high-mass analysis applies a set of kinematic requirements to select events consistent with → → 4ℓ decays. The invariant mass of the four-lepton system must be consistent with the SM Higgs boson: 115 GeV < 4ℓ < 130 GeV. Also, the lepton pairs must not be consistent with the decays of bosons ( -veto): 10 GeV < 12,34 < 64 GeV. For the 4 and 4 channels, it is possible that the leptons from a single or decay are not paired together, but rather a lepton from one / decay may be paired with a lepton from the other / decay. Therefore, there are also requirements on the alternative lepton pairings, 5 GeV < 14,23 < 75 GeV, in order to suppress * background events in which the leptons are mispaired.
Events with pairs consistent with / or Υ decay are also rejected with requirements on the lepton pair masses (see Table 2).

LM event selection
The LM analysis is designed to be sensitive to the mass range 1 GeV < < 15 GeV. For these low masses, the angular separation between the two leptons in the → ℓℓ decay can become very small (Δ (ℓ, ℓ) < 0.1 for = 1 GeV). In this case, the efficiency to select isolated electrons is significantly smaller than that for muons, so this analysis uses only the 4 final state. Otherwise, the event selection is very similar to that of the HM analysis (Section 5.4), except that a few kinematic criteria differ. The Δ requirements between final-state leptons are removed, and the -veto requirement is not relevant. In addition to the HM heavy-flavour veto, the two lepton pair masses 12 and 34 must not be in the ranges 2-4.4 GeV or 8-12 GeV. The 4ℓ requirement is narrowed to 120 GeV < 4ℓ < 130 GeV, because muons have smaller radiative losses than electrons, and both lepton pairs must satisfy 1.2 GeV < 12,34 < 20 GeV. Also, the final requirement for the signal region is simplified to 34 / 12 > 0.85.

ZX event selection
The selection for the ZX analysis differs from those of the HM and LM analyses as it is selecting a boson along with a new boson. It is, however, very similar to the selection used for the ATLAS SM → * → 4ℓ analysis [169]. In addition to the common criteria described in Section 5.3, each quadruplet must satisfy 50 GeV < 12 < 106 GeV, 12 GeV < 34 < 115 GeV, and, for the 4 and 4 channels, the alternative pairings must satisfy 14,23 > 5 GeV. The latter requirement suffices to remove mispaired / events. Backgrounds from Υ decays were found to be negligible after all selections. If there is more than one such quadruplet, quadruplets are ranked by the following criteria, applied in sequence: • Rank by the flavours of the two lepton pairs according to the reconstruction efficiencies of the leptons.
The reconstruction efficiency for muons is higher than that for electrons, especially at lower lepton momenta. Therefore, the final state lepton pairs in order of decreasing reconstruction efficiency are 4 , 2 2 , 2 2 , and 4 .
This is strictly applied to all quadruplets, even same-flavour ones, and thus differs slightly from the prescription used in the analysis of Ref. [169], where the lower-ranked alternative pairing of a same-flavour quadruplet is prevented from being considered in the quadruplet selection. For this analysis, the alternative pairing is treated as a separate quadruplet and participates in the ranking and quadruplet selection.
Following the selection of the quadruplet, the tracks associated with all four leptons are required to be consistent with originating from a common vertex: 2 / dof < 9 (tightened to < 6 for the 4 channel), where these upper bounds were chosen to give an efficiency of 99.5%. This removes additional reducible backgrounds, mainly + jets and¯. (These backgrounds are already very small for the HM and LM analyses, so this requirement is not applied in those cases.) Finally, the total invariant mass is required to be consistent with the decay of a Higgs boson, in the same manner as in the HM analysis: 115 GeV < 4ℓ < 130 GeV. None of the other requirements of the HM analysis ( boson/heavy-flavour veto and signal region requirements) are applied here.

Systematic uncertainties
Many systematic uncertainties are common to all the analyses considered here. The dominant ones include: • Luminosity and pile-up: The uncertainty in the integrated luminosity is 1.7% [170], obtained using the LUCID-2 detector [171] for the primary luminosity measurements. Uncertainty due to pile-up arises from differences between the predicted and measured inelastic cross sections, as well as from the reweighting procedure described in Section 4. This uncertainty is approximately 1%.
• Lepton-related uncertainties: The efficiency for events to pass the selection depends on the reconstruction and identification efficiencies for leptons, as well as the determination of their momentum scale. Tag-and-probe techniques are applied to the dilepton resonances → ℓ + ℓ − , / → ℓ + ℓ − , and Υ → + − in order to measure the efficiencies and momentum scales and resolutions for electrons and muons. This leads to corrections, usually of the order of up to a percent, to account for differences observed between data and simulation, as well as an estimate of the residual uncertainty [100, 102]. As there are four leptons in the final state, small single-lepton uncertainties can result in larger uncertainties in the final yields, which range up to 15%, dominated by the uncertainty in the electron reconstruction and identification efficiency.
• Theoretical uncertainties: Uncertainties in the modelling of the simulated signal and background processes are estimated by varying the parton distribution functions, the factorization, renormalization, and QCD scales, and the modelling of hadronization and the underlying event. The total uncertainty in the acceptance of the signal is around 3%, and the uncertainty in the background yield is 3%-9% for the → * → 4ℓ process [172] and about 5% for * → 4ℓ [149- 151,153,157,159].
Uncertainties related to data-driven background estimates are discussed in the analysis-specific sections below.
Each source of systematic uncertainty is considered to be uncorrelated with others: in the statistical description of the data, each source of systematic uncertainty is parameterized by several nuisance parameters that are constrained by Gaussian probability density distributions. The luminosity and lepton-related uncertainties are completely correlated among all Monte Carlo samples.

HM analysis: → → 4ℓ (15 GeV < < 60 GeV)
The high-mass analysis searches for decays of a SM Higgs boson into a pair of new bosons , where could be d , , or , which in turn decay into pairs of electrons or muons (see Figure 1). The event selection (detailed in Section 5.4) seeks two same-flavour opposite-sign pairs of leptons of similar invariant mass that are consistent with the decay of a SM Higgs boson and inconsistent with the subsequent decay of bosons.

Background estimate
Backgrounds with four prompt leptons are estimated from simulation (see Section 4) and validated using data from background-dominated control samples. The dominant backgrounds are → * → 4ℓ (about 72% of the total background) and * → 4ℓ (about 24% of the total background). Other such processes include¯→ 4ℓ and processes with three gauge bosons. These are found to be negligible.
Reducible backgrounds include those from processes with leptons originating from the decay of heavyflavour jets, or with jets misidentified as leptons. The background from the + jets process is estimated using data. Control regions enriched in misidentified leptons are defined by selecting quadruplets with one or two of its subleading leptons satisfying 'inverted' criteria but otherwise passing the signal region selection. 'Inverted' electrons fail either the Loose electron selection or the isolation requirement, but not both. 'Inverted' muons fail either the isolation requirement or the transverse impact parameter significance requirement ( 0 / 0 < 3), or both. Two samples are defined, both requiring two leptons consistent with the decay of a boson. The 'good' sample requires at least one extra lepton passing the nominal selection, while the 'inverted' sample requires at least one extra lepton passing the 'inverted' selection. Since both samples are highly enriched in → ℓℓ decays, the extra leptons originate mostly from jets misidentified as leptons. Transfer factors are defined as the ratio of the number of extra leptons passing the 'good' selection to the number passing the 'inverted' selection. These transfer factors are applied to events in the 'inverted' control regions in order to extrapolate to the signal region. The systematic uncertainties in this procedure are estimated by propagating the statistical uncertainties in the transfer factors as well as comparing the results from several different definitions of 'good' and 'inverted' leptons. This yields an estimate of the background due to the + jets process in the signal region compatible with zero.
Other reducible backgrounds are estimated from simulation. The dominant contribution is from¯, with about 3% of the total background. Other such backgrounds, including those from diboson production and heavy-flavour processes, are found to be negligible.

Background validation
The background estimates are validated using four dedicated background-enriched validation regions, defined so that they do not overlap with the HM signal region: • VR1: The -veto requirement on the alternative pairings is inverted, requiring 14,23 ≥ 75 GeV, and the compatibility requirement on 34 / 12 is removed. This produces a sample enriched in the → * → 4ℓ process as well as the non-resonant * → 4ℓ process. Only the 4 and 4 final states contribute to this region.
• VR2: The requirements on the four invariant mass pairings are removed and replaced with 12 ≥ 64 GeV, and the compatibility requirement on 34 / 12 is removed. This sample is also enriched in both the → * → 4ℓ and * → 4ℓ processes. All final states contribute to this region.
• VR4: The final 34 / 12 compatibility requirement is inverted, and all four dilepton mass requirements are changed to ℓℓ < 55 GeV. This sample mainly consists of → * → 4ℓ, but has a significant contribution from * → 4ℓ.
Although these validation regions are constructed so that they do not overlap with the signal region for the HM analysis, there is some overlap of VR1 and VR2 with the signal region of the ZX analysis. However, given the cross-section limits found for the → d → 4ℓ process (Figure 17(a)), the contribution of ZX signal events to either of these regions is less than 5% of the SM background expectation. Figure 2 compares the predicted backgrounds in these regions with the data for the variable ℓℓ = 1 2 ( 12 + 34 ). Good agreement is found in all cases. In these validation regions, the + jets background is estimated from MC simulations, while for the signal region it is estimated from data.

Results
The resulting ℓℓ distribution for this analysis is shown in Figure 3(a), while Table 3 summarizes the final yields and uncertainties in the signal region as defined in Table 2. A total of 20 events are observed, with a total predicted background of 15.6 ± 1.3 events. The -values for the background-only hypothesis as a function of are shown in Figure 4. The profile-likelihood ratio (−2 log[ ( = 0,ˆ)/ (ˆ,ˆ)]) is used as the test statistic, and the likelihood used is described in Section 10. Different final states are not distinguished in the fit; distributions used are summed over all channels. The largest deviation from SM expectations occurs around d = 28 GeV, corresponding to the two events with ℓℓ ≈ 28 GeV, with a local significance of 2.5 . Following procedures fixed before the data in the signal region were examined, the one event with ℓℓ < 15 GeV and the two with ℓℓ > 60 GeV are not considered when setting limits and do not affect Figure 4. The distribution of 34 versus 12 for the selected events is shown in Figure 3(b).      The LM analysis extends the HM analysis to the region 1 GeV < < 15 GeV, where = d , , or . Only the 4 final state is considered for this analysis. The event selection is detailed in Section 5.5 and is similar to that of the HM analysis, with some adjustments for the different kinematic region.

Background estimate
Backgrounds involving four prompt leptons are estimated directly from MC simulations (see Section 4). The → * → 4 and * → 4 processes together comprise about two thirds of the total background estimate. Higher-order electroweak processes, including triboson production and vector-boson scattering, are found to be negligible.
The remaining backgrounds involve non-prompt leptons, primarily from decays of heavy-flavour hadrons in events with multiple -quarks such as¯. A leading part of this contribution comes from double semileptonic decays, where a -hadron decays into a muon and a -hadron, which further decays into another muon and light hadrons. Resonances produced in the -hadron decay chain (i.e., , , , / ) are also an important background but are almost completely suppressed by the heavy-flavour vetoes on dilepton masses required as part of the LM event selection. There is also a small contribution from¯¯, where each muon originates from an independent -quark. As the muons selected here are all isolated, -jet tagging is not useful for reducing these backgrounds. The backgrounds from these processes are estimated together using a data-driven method [44,50].
The first step is to find the shape of the background in the 12 -34 plane. The invariant mass distribution of each muon pair is modelled separately to account for the different kinematic selections imposed on the leading, subleading, and remaining muons. Two distinct control samples are used, each of which contains an opposite-sign muon pair plus a third muon. The first sample, used to model 12 , requires a muon pair with T 1 > 20 GeV and T 2 > 10 GeV satisfying a dimuon trigger, and a third muon with T 3 > 5 GeV. The second sample, used to model 34 , requires a muon pair with T 1,2 > 5 GeV and a third muon with T 3 > 27 GeV, satisfying a single-muon trigger. In both cases, the muons must pass the same isolation and quality requirements as for the signal region. Ninety-seven percent of signal events pass both these selections, with the 12 and 34 pairs passing the requirements of the muon pair in the first and second samples, respectively. The invariant masses of the muon pairs are taken from the two control samples and used to form a 2D template in the 12 -34 plane as the direct product of the two distributions.
A correction to the 12 -34 template is made to account for a correlation between the kinematics of the two muon pairs, which is introduced by the Higgs boson mass requirement. Another control sample is defined by inverting the isolation and vertex requirements on the muons in the signal event selection, defining a sample enriched in events with muons from heavy-flavour quark decays. Comparing the distributions of muon pair invariant masses before and after the Higgs boson mass requirement yields the correction to the background shape as a function of 12 and 34 .
Finally, the overall normalization for the background from non-prompt leptons is determined from data in regions defined by inverting several selection criteria. As shown in Figure 5, region B is defined by inverting the compatibility requirement 34 / 12 > 0.85. In order to improve the statistical precision of the background prediction, additional regions are defined by inverting the Higgs boson mass requirement (region C in Figure 5) and also the isolation and vertex requirements (regions D and E in Figure 5). The regions with 81 GeV < 4ℓ < 101 GeV are excluded in order to reduce contributions from bosons. The contribution with prompt muons, mostly * in regions B and C, is subtracted from the data. The background with non-prompt leptons in region B is then estimated using = · / . The 2D template is then used to scale from the estimate in region B to the signal region satisfying the compatibility requirement, region A. The uncertainty in the heavy-flavour background estimate is found by varying each parameter of the background shape model up and down by 2 and taking the largest change in yield for each bin, giving an uncertainty of 38%. Statistical uncertainties in the normalization of the signal region are also propagated to the heavy-flavour background yield, giving an uncertainty of 33%. Adding these uncertainties in quadrature gives a total systematic uncertainty of 50% in the heavy-flavour background yield.

Results
The ℓℓ distribution in the LM signal region is shown in Figure 6(a). The distribution of 12 vs 34 is shown in Figure 6(b), while Table 4   invariant mass distribution of the other pair. For d > 55 GeV, the invariant mass distribution for d → ℓℓ starts to overlap significantly with that for → ℓℓ. Since this analysis accepts events with the leading lepton pair consistent with the decay of a boson, it relies much more on the invariant mass distribution of the other lepton pair to reject the → 4ℓ background than does the HM analysis. Therefore, the upper search range for this analysis is limited to 55 GeV, rather than 60 GeV as for the HM analysis.

Background estimate
The dominant backgrounds in this analysis are → * → 4ℓ (about 65% of the total) and non-resonant * → 4ℓ (about 33% of the total). Additional prompt backgrounds include the triboson processes , , and . These are estimated from simulation (see Section 4), but the * → 4ℓ background estimate is checked using background-enriched validation samples.
Other, reducible, backgrounds, such as those from + jets,¯, and processes, contain either additional non-isolated leptons from heavy-flavour decay or objects misidentified as leptons and constitute only a few percent of the background. The procedure used to estimate the total yield of these backgrounds is identical to that of the ATLAS SM → * → 4ℓ analysis [169,173].
The reducible background is estimated separately for the cases where the second lepton pair ( 34 ) decays into muons (ℓℓ ) and those where it decays into electrons (ℓℓ ). For the ℓℓ case, a number of mutually exclusive control regions are defined by inverting or relaxing some of the lepton identification requirements, including the isolation and impact parameter requirements for the subleading muon pair. A fit to the 12 distribution is then performed to estimate the amount of background due to each of¯, + heavy-flavour (having -or -quark content), and + light-flavour. Transfer factors derived from simulation are then used to extrapolate the fitted yield of each background in the control regions to the signal region. The contribution from production is estimated using simulation.
The ℓℓ background from + jets,¯, and production is classified into processes with jets being misidentified as electrons ( ), electrons from photon conversions ( ), and electrons from semileptonic decay of heavy-flavour hadrons ( ). The component is estimated from simulation. The other two components are estimated from a control region in which the identification requirements of the lowest-T electron are relaxed. Further, to suppress the * contribution, the two subleading electrons must have the same sign. The expectations for the and components are obtained by fitting to the distribution of the number of inner pixel detector hits associated with the track of the lowest-T electron. The estimated yields of all three components are then extrapolated to the signal region using transfer factors derived from simulation.
Finally, the shape of the 34 distribution for the reducible background is taken from simulation. An additional 10% systematic uncertainty is assigned to the reducible background estimate to account for differences in the lepton isolation requirements between this analysis and that of Refs. [169,173].
These requirements are illustrated in Figure 7. Distributions of 34 for the two validation regions are shown in Figure 8. Good agreement is found with background expectations.

Results
The final 34 distribution for this analysis is shown in Figure 9, while Table 5 summarizes the final yields and uncertainties. The dominant systematic uncertainty in final states that contain electrons arises from the modelling of the electron identification efficiency. For the 4 channel, the dominant systematic uncertainty arises from the modelling of muon isolation. A total of 356 events are observed with an expected background of 320 ± 17. Figure 10

Limits and interpretation
No significant excess is observed above SM background predictions for any of the analyses considered. Therefore, the results are interpreted in terms of exclusion limits. Firstly, model-independent limits are placed on fiducial cross sections. Model-dependent exclusion limits are then set for the benchmark models described in Section 2.
For the HM and LM → → 4ℓ analyses, evaluating the limits entails parameterizing the signal distribution as a function of both ℓℓ and , while the → → 4ℓ analysis requires the parameterization to be a function of 34 and . Since simulated events are generated only at discrete values of , the signal templates are interpolated between values. For the HM and ZX analyses, this is done using moment morphing [175]. The distributions at the generated values of are used as templates, and the normalization is determined from interpolation of the simulated signal yields. For the LM analysis, Gaussian distributions are fit to the ℓℓ distributions at each generated , and the fit parameters are interpolated in .
The data are described statistically by a likelihood function consisting of a Poisson factor for each histogram bin, summed over each channel, along with a Gaussian constraint for each nuisance parameter [176]: where is the number of observed events observed in bin for channel , is the set of nuisance parameters, ( ) and ( ) are the predicted numbers of signal and background events for each bin and channel, is the signal strength, and and are mean and width of the Gaussian constraint for nuisance parameter . Systematic uncertainties are modelled via nuisance parameters which are profiled in the calculation of the test statistic; the effect of systematic uncertainties on the limits is small.

Limits on fiducial and total cross sections
Model-independent cross-section limits for the HM, LM, and ZX analyses are derived in fiducial regions defined using generator-level quantities. These fiducial selections, shown in Table 6, are designed to mimic the signal region selection requirements. In order to account for the effects of quasi-collinear electromagnetic radiation from the leptons within the detector resolution, the four-momenta of prompt photons close to a lepton (Δ < 0.1) are added to the four-momentum of that lepton [177].
The fiducial selections are used to factorize the effects of the event selection into a largely model-independent 'efficiency' and a model-dependent 'acceptance'. The efficiency for a given channel is defined as the fraction of events passing the fiducial selection (using generator-level quantities) that also pass the full event selection (using reconstructed quantities). This mostly depends on the lepton reconstruction, but not on the model used. Systematic uncertainties relevant to the reconstruction of leptons are propagated to the efficiency. For a given theoretical signal model, the acceptance for a channel is defined as where fid is the yield for channel within the fiducial region (at generator level) and tot is the total yield for channel (simply the total generator-level yield for the channel). The efficiency may thus be used to find a model-independent fiducial cross-section limit, which may then be converted to a model-dependent total cross-section limit using the acceptance.

HM and LM limits
The efficiencies within the fiducial regions for the HM and LM analyses are shown in Figure 11(a). These were calculated using the benchmark → d d model, but the efficiencies are mostly model-independent: for → → 4 over the range 1 GeV < < 15 GeV the efficiencies are the same as for → d d → 4 to within a relative difference of 3%. The difference in efficiency between different final states is mainly due to the fact that the efficiencies for reconstruction, identification, and selection are lower for electrons than for muons. These efficiencies are used to compute 95% confidence level (CL) upper limits on the cross section within the fiducial region, using the CL s frequentist formalism [178] with the profile-likelihood-ratio test statistic [179]. The resulting limits are shown in Figure 12. These limits should be applicable to any models of the SM Higgs boson decaying into four leptons via two intermediate bosons that are narrow, on-shell, and that decay promptly. The model-dependent acceptances for the HM and LM analyses are shown in Figure 11(b) for the → d d and → → 4 models. The resulting upper limit on the product of the total cross section and decay branching ratio for the benchmark model, , for the HM analysis is shown in Figure 13, while Figure 14 shows upper limits on ( → → d d → 4 ) and ( → → → 4 ) for both the HM and LM analyses. These results are independent of assumptions about the decay branching ratios of the d and bosons. In particular, Figure 14(b) also applies to the scalar case ( → → → 4 ).

ZX limits
For limits involving ZX processes, the normalization of the non-resonant * → 4ℓ background is validated using control samples, but the normalization of the remaining significant background, → * → 4ℓ, is allowed to float in the limit determination as an unconstrained nuisance parameter. The model-independent efficiency within the fiducial region is shown in Figure 15(a), and the resulting 95% CL upper limit on the fiducial region cross section is shown in Figure 16. The fiducial region acceptance for the → d → 4ℓ process is shown in Figure 15(b), and the upper limits on the product of the total cross section and decay branching ratio for the benchmark models, ( → → d → 4ℓ) and ( → → → 2ℓ2 ), are shown in Figure 17.

Limits on branching ratios
A (model-dependent) cross-section limit may be converted to a branching ratio limit using the relations: where → →4ℓ is the model-dependent total cross section, is the SM Higgs boson production cross section for the ggF process (48.58 pb for = 125 GeV [124]), and B ( → 2ℓ) is the model-dependent branching ratio for each decay to one lepton flavour. The branching ratios for d → ℓℓ and → are taken from the benchmark models [21,22], where for the d → ℓℓ case, the branching ratios for the two lepton flavours are taken to be equal. For the → case, the branching ratio varies considerably in a model-dependent way over the range of considered here. The resulting branching ratio limits are shown in Figure 18.

Limits on Higgs mixing
The branching ratio limit can also be interpreted as a limit on the effective Higgs mixing parameter , defined as where is the Higgs portal coupling and is the mass of the dark Higgs boson. Using rather than combines the dependencies on and into a single parameter. Then, according to Eq. (2.33) of Ref. [21] and assuming > /2: where Γ SM is the SM width of the 125 GeV Higgs boson, and ≈ 246 GeV is the vacuum expectation value of the Higgs field. The resulting limit is shown in Figure 19.
The → d analysis can also be used to set limits on the d mixing parameter and on thed mass mixing parameter , as described in Refs. [21,43]. These are shown in Figure 20, assuming the SM Higgs boson production cross section.  Expected   The limits presented in this paper improve on those from the previous ATLAS search by factors between two and four due to a larger data sample, improved lepton reconstruction and identification, and a better optimized event selection. In addition to the improvements on the results from the previous search, this paper also presents limits on total cross sections and on the dark Higgs boson mixing parameters.

Appendix
The signal region for the HM analysis in Ref. [44] was defined by 34 / 12 > 0.85. This was the result of an optimization that assumes the d width is narrow, as expected in the HAHM, so that the observed width will be dominated by the detector performance. A ±2 change in the energy measurement gives about a 15% change in the lepton-pair invariant mass, motivating the coefficient of 0.85 in the selection. However, the background is low, especially in the lower mass region, so it is possible to widen the signal region somewhat relative to the background without significant loss of sensitivity. This is further motivated by models such as the one discussed in Chapter 7 of Ref. [181]. Accordingly, the signal region selection was re-optimized to allow for a larger d width.
The selection widens up to 3.5 at the low end of the dilepton mass spectrum, where the background is the lowest, but narrows to the previous value of 2 for higher masses. So the form of the selection is 34 where ( 12 ) is a modulating function that is 1 at the lowest mass and goes to 0 at higher masses.
The form of the modulating function is taken from the shape of the average dilepton mass spectrum in the background. The function fit to this shape consists of an exponential tail matched with half of a Gaussian: The fit parameters, with given in GeV, are ℎ = 3.73, = 51.6, = 16.6, 1 = −2.62, 2 = −0.0266, and = 6.39, with the maximum of ( ) occurring at max = 49.64 GeV. The modulating function is then defined by The final shape of the re-optimized signal region is shown in Figure 3(b).