Search for long-lived particles decaying to a pair of muons in proton-proton collisions at $\sqrt{s}$ = 13 TeV

An inclusive search for long-lived exotic particles decaying to a pair of muons is presented. The search uses data collected by the CMS experiment at the CERN LHC in proton-proton collisions at $\sqrt{s}$ = 13 TeV in 2016 and 2018 and corresponding to an integrated luminosity of 97.6 fb$^{-1}$. The experimental signature is a pair of oppositely charged muons originating from a common secondary vertex spatially separated from the pp interaction point by distances ranging from several hundred $\mu$m to several meters. The results are interpreted in the frameworks of the hidden Abelian Higgs model, in which the Higgs boson decays to a pair of long-lived dark photons Z$_\mathrm{D}$, and of a simplified model, in which long-lived particles are produced in decays of an exotic heavy neutral scalar boson. For the hidden Abelian Higgs model with $m_\mathrm{Z_D}$ greater than 20 GeV and less than half the mass of the Higgs boson, they provide the best limits to date on the branching fraction of the Higgs boson to dark photons for $c\tau$(Z$_\mathrm{D}$) (varying with $m_\mathrm{Z_D}$) between 0.03 and ${\approx}$ 0.5 mm, and above ${\approx}$ 0.5 m. Our results also yield the best constraints on long-lived particles with masses larger than 10 GeV produced in decays of an exotic scalar boson heavier than the Higgs boson and decaying to a pair of muons.


Introduction
Long-lived particles (LLPs) are predicted by many extensions of the standard model (SM), in particular by various supersymmetric scenarios [1,2] and "hidden sector" models [3,4]. Such particles could manifest themselves through decays to SM particles at macroscopic distances from the proton-proton (pp) interaction point (IP). This paper describes an inclusive search for an exotic massive LLP decaying to a pair of oppositely charged muons, referred to as a "displaced dimuon", that originates from a common secondary vertex spatially separated from the IP. The search is based on an analysis of pp collisions corresponding to an integrated luminosity of 97.6 fb −1 collected with the CMS detector at √ s = 13 TeV during Run 2 of the CERN LHC. A minimal set of requirements and loose event selection criteria allow the search to be sensitive to a wide range of models predicting LLPs that decay to final states that include a pair of oppositely charged muons. We interpret the results of the search in the framework of two benchmark models: the hidden Abelian Higgs model (HAHM), in which displaced dimuons arise from decays of hypothetical dark photons [5], and a simplified model, in which a non-SM Higgs boson decays to a pair of long-lived exotic heavy neutral scalar bosons, at least one of which decays into a pair of muons [6].
The present search explores the LLP mass range above 10 GeV accessible with the Run 2 dimuon triggers and is sensitive to secondary vertex displacements ranging from several hundred µm to several meters. It is a continuation and extension of two CMS analyses performed using data taken at √ s = 8 TeV during Run 1 of the LHC. One analysis was dedicated to a search for LLPs decaying to two electrons or two muons in the tracker [7]; the other looked for LLP decays to final states containing two muons reconstructed only in the muon system [8]. The analysis of the Run 2 data described here contains numerous improvements over these Run 1 searches, notably in refined event selection and improved background evaluation procedures. It also benefits from an increase in integrated luminosity by almost a factor of five, collected at a higher √ s. A search for LLPs decaying to displaced dimuons has also been performed by the ATLAS Collaboration, using 2016 data corresponding to an integrated luminosity of 32.9 fb −1 [9]. This paper is organized as follows. Section 2 describes the CMS detector. Section 3 presents the signal models as well as the samples analyzed from data and from the Monte Carlo simulation. Section 4 describes the analysis strategy, the triggers, and the offline event selection. Estimation of backgrounds and the associated systematic uncertainties are described in Section 5. Section 6 summarizes the systematic uncertainties affecting signal efficiencies. Section 7 describes the results obtained in the individual dimuon categories and their combination. The analysis summary is presented in Section 8. Tabulated results and supplementary material for reinterpreting the results in the framework of models not explicitly considered in this paper are provided in HEPData [10].

The CMS detector
The central feature of the CMS detector 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 extending outwards to a radius of 1.1 m, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the coverage in pseudorapidity η provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers covering the range |η| < 2.4 and embedded in the steel flux-return yoke outside the solenoid. The muon system is composed of three types of chambers: drift tubes (DTs) in the barrel, cathode strip chambers (CSCs) in the endcaps, and resistive-plate chambers in both the barrel and the endcaps. The chambers are assembled into four "stations" at increasing distance from the IP; each station provides reconstructed hits in several detection planes, which are combined into track segments, forming the basis of muon reconstruction in the muon system [11]. 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. [12].
Events of interest are selected using a two-tiered trigger system. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of approximately 100 kHz within a fixed time interval of less than

Signal models, data and simulated samples
The search is performed using pp collision data collected at √ s = 13 TeV in 2016 and 2018 corresponding to integrated luminosities of 36.3 ± 0.4 and 61.3 ± 1.5 fb −1 , respectively [15,16]. The data collected in 2017 are not used because the triggers required for the analysis were not included when those data were recorded.
As mentioned above, two signal models with different final-state topologies and event kinematics are used in the optimization of event selection criteria and in the interpretation of results. The first belongs to a class of models featuring a "hidden" or "dark" sector of matter that does not interact directly with SM particles, but can manifest itself through mixing effects. This HAHM benchmark contains an extra dark U(1) D gauge group whose symmetry is broken by a new dark Higgs field [5,17]. The spin-1 mediator of the U(1) D group, known as the dark photon Z D , mixes kinetically with the hypercharge SM gauge boson ("vector portal"), whereas the dark Higgs boson H D mixes with the SM Higgs boson H ("Higgs portal") and gives mass m(Z D ) to the dark photon. If there are no hidden-sector states with masses smaller than m(Z D ), the mixing through the vector portal with the SM photon and Z boson causes the dark photon to decay exclusively to SM particles, with a sizable branching fraction to leptons. Pair production of the Z D via the Higgs portal with subsequent decays of dark photons via the vector portal is shown in Fig. 1 (left).
The present search probes the regime of m(Z D ) > 10 GeV with small values of the Z-Z D kinetic mixing parameter ϵ [5]. In this regime, the dark photon is long-lived, since its mean proper lifetime τ(Z D ) is proportional to ϵ −2 . In particular, the dark photon with 10 GeV ≲ m(Z D ) < m(H)/2 is expected to have macroscopically large mean proper decay lengths cτ(Z D ) ≳ O(100 µm) for ϵ < O(10 −6 ). The Z D production rate is governed by the branching fraction B(H → Z D Z D ), which does not depend on ϵ but is proportional to the square of κm 2 (H)/|m 2 (H) − m 2 (H D )|, where κ is the H-H D mixing parameter. Since κ and m(H D ) affect only the overall dark photon production rate, sampling of m(Z D ) and ϵ is sufficient to explore different kinematical and topological scenarios of the model. We generated a set of 24 HAHM samples with m(Z D ) between 10 and 60 GeV and ϵ between 10 −6 and 2 × 10 −9 . In this mass range, the model's prediction for B(Z D → µµ) varies between 15.4% at m(Z D ) = 10 GeV and 10.7% at m(Z D ) = 60 GeV. The dark Higgs boson is assumed to be heavy enough such that H → H D H D decays are kinematically forbidden. In the sample generation, we specify m(H D ) = 400 GeV and κ = 0.1. The production of dark photons is modeled at leading order by MADGRAPH5 aMC@NLO [18] H κμ μ The generation of the samples is done for the dominant gluon-fusion production mechanism. The Higgs boson production cross section is normalized to the most recent theoretical prediction for the sum of all production modes for m(H) = 125 GeV, 55.7 pb [19]. The decays of the dark photons are modeled by PYTHIA 8.212 and 8.230 [20] in samples corresponding to the 2016 and 2018 data sets, respectively.
At the LHC, another way that LLPs might arise is via production of mediators heavier than the Higgs boson that decay into LLPs. To explore ranges of kinematic variables and event topologies broader than those offered by HAHM, we also consider a simplified benchmark model [6], previously used in Run 1 searches for displaced dimuons [7,8], in which the LLP is an exotic spin-0 boson X. The scalar X has a non-zero branching fraction to dimuons and is pair produced in the decay of a new heavier scalar boson Φ, which is produced in gluon-gluon fusion: gg → Φ → XX, X → µ + µ − . The Feynman diagram for this process is shown in Fig. 1 (right).
The samples for Φ → XX are generated with PYTHIA. Two sets of samples are produced, depending on whether one or both X bosons are forced to decay to dimuons. Samples in each set are generated with different combinations of Φ boson masses m(Φ) (ranging from 125 GeV to 1 TeV) and X boson masses m(X) (ranging from 20 to 350 GeV). The width of the Φ boson is assumed to be small for the purpose of simulation, but the analysis has negligible dependence on this assumption. Each sample is furthermore produced with three different mean proper lifetimes τ(X) of the X bosons, corresponding to mean transverse decay lengths of approximately 3, 30, and 250 cm. Events generated at the selected values of m(Φ), m(X), and τ(X) allow us to study wide ranges of signal displacements, kinematical variables, and event topologies.
Since the optimization of the event selection criteria and the evaluation of the residual backgrounds are performed using data, the simulated background samples are used primarily to gain a better understanding of the nature and composition of surviving background events. Simulated background samples used in the analysis include Drell-Yan (DY) dilepton production; tt, tW, and tW events; W and Z boson pair production (dibosons); W+jets; and events comprised of jets produced through the strong interaction that are enriched in muons from semileptonic decays of hadrons containing b or c quarks.
The 2016 simulated signal and background samples are produced with either the NNPDF2.3 (leading order) or NNPDF3.0 (next-to-leading order) parton distribution functions (PDFs) [21], using the CUETP8M1 [22] tune to model the underlying event. All 2018 simulated samples are produced with the NNPDF3.1 PDFs [23] (next-to-next-to-leading order), using the CP5 [24] tune, which is optimized for the NNPDF3.1 PDFs. Simulation of the passage of particles through detector material is performed by GEANT4 [25]. Simulated minimum bias events are superimposed on a hard interaction in simulated events to describe the effect of additional inelastic pp interactions within the same or neighboring bunch crossings, known as pileup; the samples are weighted to match the pileup distribution observed in data. All simulated events are then reconstructed with the same algorithms as used for data.

Analysis strategy and event selection
An LLP produced in the hard interaction of the colliding protons may travel a significant distance in the detector before decaying into muons. While trajectories of the muons produced well within the silicon tracker can be reconstructed by both the tracker and the muon system, tracks of muons produced in the outer tracker layers or beyond can only be reconstructed by the muon system. Since the dimuon vertex resolution and the background composition differ dramatically depending on whether the muon is reconstructed in the tracker, we classify all reconstructed dimuon events into three mutually exclusive categories: a) both muons are reconstructed using both the tracker and the muon system (TMS-TMS category); b) both muons are reconstructed using only the muon system, as "standalone" muons (STA-STA category); and c) one muon is reconstructed only in the muon system, whereas the other muon is reconstructed using both the tracker and the muon system (STA-TMS category). These three categories of events are analyzed separately, each benefiting from dedicated event selection criteria and background evaluation. The results in each category are statistically combined to provide the final results.
The beamspot is identified with the mean position of the pp interaction vertices. The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4 of Ref. [26]. A pair of reconstructed muon tracks is fitted to a common vertex (CV), which is expected to be displaced with respect to the PV. The transverse decay length L xy is defined as the distance between the PV and the CV in the plane transverse to the beam direction. The transverse impact parameter d 0 is defined as the distance of closest approach (DCA) of the muon track in the transverse plane with respect to the PV.
Events were collected with dedicated triggers aimed at recording dimuons produced both within and outside of the tracker. Therefore, these triggers require two muons reconstructed in the muon system alone, without using any information from the tracker, and do not impose the beamspot constraint in the muon track fit at the HLT [27]. Each muon is required to be within the region |η| < 2.0 and to have transverse momentum magnitude p T > 28(23) GeV in 2016 (2018) data taking. To reduce the trigger rate caused by hadron punch-through and poorly measured muons, each muon track is required to be composed of segments found in two or more muon stations. To reduce the contribution to the trigger rate from cosmic ray muons and low-mass dimuon resonances, the trigger used to collect 2016 data also required that the 3D angle between the muons be less than 2.5 rad, and that the invariant mass of the two muons be larger than 10 GeV. The optimization of the online selection prior to the 2018 data taking made it possible to remove these two requirements from the 2018 trigger, thus providing additional validation regions for background evaluation and increasing the signal efficiency. The efficiency of triggering on signal events in 2018 was further improved by complementing the above trigger with one very similar to it, but using a modified version of the initial "seeding" stage of the muon trajectory building at the HLT. The seed generator used in the new trigger was specifically designed for muons not pointing to the beamspot and helped to increase the reconstruction efficiency for displaced muons.
The high-level triggers used in the analysis were seeded by L1 dimuon triggers that required the p T of the muons to be above certain thresholds. The values of the thresholds were varied during the data taking, depending on the instantaneous luminosity, from 11 and 4 GeV (for the leading and subleading L1 muons, respectively) during most of 2016, to 15 and 7 GeV at the end of Run 2. At L1, the p T assignment for muons was made under the assumption that the muons originated at the beamspot. As a result, the p T of the displaced muons not pointing to the beamspot were underestimated and could fall below the L1 trigger thresholds. The ensuing signal efficiency loss was larger when higher L1 trigger p T thresholds were used. Since this effect is decoupled from the collision environment (e.g., instantaneous luminosity), it can be studied using cosmic ray muons recorded with very loose triggers during periods with no beam. Figure 2 shows that the decrease in the L1 muon trigger efficiency as a function of the impact parameter d 0 for various L1 trigger p T thresholds used in 2016 (left) and 2018 (right) is well reproduced by the signal simulation in the barrel.  Figure 2: L1 muon trigger efficiency in cosmic ray muon data (blue) and signal simulation (red) as a function of d 0 , for the L1 trigger p T thresholds used in (left) 2016 and (right) 2018. The denominator in the efficiency calculations is the number of STA muons with |η| < 1.2 and p T > 33 (28) GeV in 2016 (2018).
As noted in Section 1, no single muon reconstructor provides optimal performance over the wide range of displacements of secondary vertices considered in the analysis. Muons produced near the IP can be accurately reconstructed by using commonly used algorithms developed for prompt muons and combining measurements in the tracker and the muon system. Among them are the global muon and tracker muon reconstruction algorithms [11,28]. The first algorithm builds "global muons" by using hits in the tracker and segments in the muon system in a common track fit. The second constructs "tracker muons" by propagating tracks in the inner tracker to the muon system and requiring loose geometrical matching to DT or CSC segments. The efficiency of these algorithms, however, rapidly decreases as the distance between the IP and muon origin increases, dropping to zero for muons produced in the outer tracker layers and beyond. On the other hand, such muons can still be efficiently reconstructed by algorithms that use only information from the muon system. These STA algorithms [11,28] can reconstruct muons with displacements of up to a few meters, but they have poorer spatial and momentum resolution than muons reconstructed using more precise information from the silicon tracker.
To benefit from the advantages offered by both types of algorithms and to follow what was done in the trigger, we begin the muon selection with the most efficient standalone muons and replace them with more accurately reconstructed global and tracker muons whenever such muons are found. We use muons reconstructed by an STA algorithm with the beamspot constraints removed from all stages of the muon reconstruction procedure, which yields the highest efficiency and the best resolution for displaced muons, out of all available STA algorithms. The event selection starts with the requirement that the event is selected by the triggers described above and has at least two STA muons, each containing more than 12 valid CSC or DT hits. The requirement of the minimal number of hits suppresses backgrounds from hadron punch-through and other sources, and ensures that the STA muons have acceptable p T resolution and charge assignment. The STA muons that satisfy this basic quality requirement form the initial list of the muon candidates retained for the analysis.
Next, we reject events in which no HLT muon pair that triggered the event matches two STA muons in the list. This requirement suppresses events that triggered on muons not related to the signal and facilitates application of trigger efficiency measurements in the analysis. We then attempt to match each STA muon in the list with a TMS muon, i.e., a global or a tracker muon. The STA and TMS muons are considered to be matched if they share at least two thirds of their segments or if ∆R STA−TMS < 0.1, where ∆R STA−TMS = √ (η hit − η pca ) 2 + (ϕ hit − ϕ pca ) 2 is the separation between η hit (ϕ hit ) of the position of the innermost hit of the STA muon and η pca (ϕ pca ) of the point of closest approach of the TMS muon to this hit. If an associated TMS muon is found, it replaces the corresponding STA muon in the list of the muon candidates used for further analysis; otherwise, the original STA muon is kept. The matching procedure was optimized using events in the simulated signal and background samples as well as data in the signal-free control regions discussed in Section 5. It eliminates most of the pp collision background to LLP decays outside of the tracker and greatly increases sensitivity to LLP decays in the tracker, thanks to a far superior resolution of TMS muons compared to that of STA muons.
The impact of the STA-to-TMS muon association procedure on the event selection is further illustrated in Fig. 3, which shows the fraction of simulated Φ → XX → µµ+anything signal events with zero, one, and two STA muons matched to TMS muons as a function of L true xy , defined as the transverse distance between the simulated positions of the hard-interaction and LLP decay vertices. While almost all dimuons produced close to the IP have both STA muons matched to TMS muons, the fraction of these events rapidly decreases with L true xy , reflecting the dependence on L true xy of the tracker reconstruction efficiency. Events with one STA muon matched to a TMS muon start to dominate at L true xy = 25 cm and remain the dominant component up to ≈50 cm, where events with no STA-to-TMS matches take over. When LLPs decay in the outer tracker layers or beyond the tracker, all STA-to-TMS associations are purely accidental and occur for fewer than 5% of the simulated signal muons. Therefore, the association procedure gives rise to three categories of dimuons, each dominating in a certain L true xy range: TMS-TMS at small L true xy , STA-TMS at intermediate L true xy , and STA-STA at large L true xy . The number of STA-STA dimuons beyond the solenoid, at L true xy > 3.2 m, is low because of the low trigger efficiency.
The STA and TMS muons are then subjected to additional selection criteria optimized using simulated signal and background samples, and samples of dimuons misreconstructed as displaced in the signal-free regions in data. The STA muons are required to have p T > 10 GeV and to satisfy the following criteria: • relative p T uncertainty σ p T /p T < 1.0, where σ p T is the internal uncertainty from the muon track fit;  Figure 3: Fractions of signal events with zero (green), one (blue), and two (red) STA muons matched to TMS muons by the STA-to-TMS muon association procedure, as a function of true L xy , in all simulated Φ → XX → µµ+anything signal samples combined. The fractions are computed relative to the number of signal events passing the trigger and containing two STA muons with more than 12 muon detector hits and p T > 10 GeV matched to generated muons from X → µµ decays.
• χ 2 /dof of the muon track fit less than 2.5; • more than 18 DT hits for muons reconstructed only in the barrel; • time difference |∆t| with respect to the current bunch crossing less than 12 ns; • muon identified as traveling outwards based on the timing measurements [11].
The σ p T /p T , χ 2 /dof, and DT hits requirements suppress background events arising from poorly measured prompt muons and ensure acceptable p T resolution and charge assignment for signal muons. The muon timing requirement serves several purposes: it rejects events in which the trigger timing is early, causing the tracker hits that are read out to be from unrelated pp collisions 25 ns earlier and not to include any tracker hits of the triggering STA muons; it rejects out-of-time collision muons; and it helps suppress background arising from cosmic ray muons crossing the detector outside of the tracker acceptance, which are often reconstructed as two STA muons with no associated TMS muons. The |∆t| < 12 ns requirement does not impact signal efficiency in the studied mass range. The muon direction, determined using the muon timing measurements, offers another handle for suppressing cosmic ray muon background by identifying cosmic ray muons in the upper part of the detector, which travel inwards. We do not impose any explicit requirements on the isolation of the STA muons, thus making the STA-STA analysis sensitive to models predicting highly displaced b quarks, such as LLPs decaying to bb.
The initial list of TMS muon candidates consists of the global and tracker muons with p T > 10 GeV that are matched to STA muons by the STA-to-TMS muon association procedure. These TMS muon candidates are further required to satisfy the following criteria: • σ p T /p T < 1.0; • be identified as tracker muons with two or more matched CSC or DT segments; • be isolated from other activity in the event.
The first two requirements reject poorly measured muons and misidentified hadrons. The isolation requirement suppresses copious background from dijet and multijet events. Specifically, we calculate the sum of the p T of other tracks reconstructed within a cone of radius ∆R = (∆η) 2 + (∆ϕ) 2 centered on the direction of the TMS muon track. Both the p T of the muon track itself and the p T of the track of any other muon forming a dimuon (if it lies within the cone) are not included in the sum. For the muon to be considered isolated, the relative tracker isolation I rel trk , defined as the ratio of the p T sum in a cone of ∆R < 0.3 to the muon track p T , is required to be smaller than 0.075. The cone size, isolation threshold, and other parameters of the algorithm such as the minimum |∆z| between the reference point of the track and that of the muon, are the result of the optimization procedure that uses displaced muons in the simulated signal samples and background muons in the signal-free control regions in data.
Dimuon candidates are then formed from pairs of STA and TMS muons passing the above muon selection criteria. All possible combinations of muons (STA-STA, STA-TMS, and TMS-TMS) are considered. Each pair of selected muons with the distance of closest approach of the muon tracks smaller than 15 cm is fit to a common vertex by means of the Kalman vertex fitter algorithm [29,30] and forms a dimuon. If more than two muons pass the selection criteria, multiple dimuons can be formed, including dimuons sharing one of the two muons; this situation would be particularly common in signal events with two LLPs decaying to two muons each. Using simulated signal samples, we developed a set of muon pairing criteria that efficiently choose up to two dimuons among the formed dimuons. The selection is first done separately in each of the three dimuon categories. In events with three selected muons, we choose a dimuon with the lowest value of χ 2 of the vertex fit, χ 2 vtx ; in events with four or more selected muons, we choose the pair of dimuons whose χ 2 vtx sum is the smallest among all distinct pairs of dimuons formed from the four highest p T muons in the event. In rare cases, when more than two dimuons are found in all categories combined, up to two dimuons with the smallest χ 2 vtx are kept for further analysis; we also ensure that there are no STA or TMS muons included in more than one dimuon. The efficiency of choosing the right dimuon is 85-95% for 4µ signal samples and close to 100% for 2µ samples.
To suppress background from dimuons formed from unrelated muons, the dimuons chosen by the muon pairing criteria are required to have χ 2 of the common vertex fit below a certain threshold. Since χ 2 distributions in data are difficult to reproduce accurately in the simulation, the requirements on χ 2 vtx are chosen to be loose, leaving the signal efficiency almost intact while rejecting a sizable fraction of background events in the control regions in data. The values of the thresholds are determined separately in each of the dimuon categories and are χ 2 vtx < 10 in the STA-STA and TMS-TMS categories, and χ 2 vtx < 20 in the STA-TMS category. To further ensure that the muon tracks are consistent with the vertex, we require that the number of tracker hits assigned to the track of the TMS muon upstream of the vertex position does not exceed 2 in the TMS-TMS category and 5 in the STA-TMS category. In the TMS-TMS category, background from dimuons formed from unrelated or mismeasured muons is further suppressed by requiring TMS muons forming the dimuon to have a similar number of pixel hits, namely that the difference between the number of pixel hits on two TMS muons, ∆N(pixel hits), is smaller than 3 hits.
A cosmic ray muon crossing the detector within the acceptance of the muon system is often reconstructed as two back-to-back muons, one in the upper half and one in the lower half of the detector. This background is very efficiently suppressed by requiring that the 3D opening angle α between the two muons be less than a certain threshold value. In the analysis of 2016 data, we require α < 2.5 rad (or, equivalently, cos α > −0.8), consistent with the presence of the cos α > −0.8 requirement in the 2016 version of the trigger. Since the 2018 version of the trigger has no cos α requirement, the offline requirements on α in the analysis of 2018 data are re-optimized in each of the dimuon categories separately, taking into account the angular resolution of STA and TMS muons. As the result of this optimization, the 2018 selection criteria are loosened to α < 2.7 rad (cos α > −0.9) for STA-STA and STA-TMS dimuons and to α < 3.0 rad (cos α > −0.99) for TMS-TMS dimuons.
The cos α requirements fail to reject one particular type of cosmic ray muon background occurring when a cosmic ray muon with a large incoming angle crosses the muon system diagonally in the r-z plane. A cosmic ray muon of this type can lead to multiple STA muons, each reconstructed from a small number of segments. These poorly measured STA muons are not necessarily back-to-back, and can give rise to a mistakenly formed displaced dimuon passing the full event selection. To suppress such events, we reject STA-STA dimuons if the sum of the numbers of segments belonging to the two muons, N(dimuon segments), is smaller than 5.
Furthermore, there are events with multiple cosmic ray muons produced in an atmospheric shower. Such events typically result in a large number of nearly parallel STA muons reconstructed in both upper and lower halves of the detector. While dimuons formed from muons in different hemispheres are rejected by the cos α requirement, dimuons formed from two muons in the same hemisphere are not. To reject the background from cosmic ray muon showers, we require that the event contain fewer than 4 nearly parallel STA muons and have at least one reconstructed pp interaction vertex with more than 3 associated tracks and with transverse (longitudinal) coordinates within 2 (24) cm of the IP. In the STA-STA dimuon category, in which the contribution from the cosmic ray muon showers is the largest, we also require that neither of the two muons forming a dimuon be back-to-back (cos α < −0.9) with another STA muon with p T > 10 GeV. The dimuon is rejected only if the time difference |∆t b2b | between the muon in the dimuon and its back-to-back muon is larger than 20 ns, i.e., consistent with the time difference between the two reconstructed parts of a cosmic ray muon. This requirement also helps to suppress dimuons reconstructed from a cosmic ray muon and a muon from an overlapping pp collision.
Displaced dimuons produced in decays of nonprompt J/ψ mesons and other low-mass SM resonances, or formed from the products of the b hadron cascade decays (b → cµ 1 X followed by c → µ 2 X), are suppressed by requiring that the reconstructed dimuon invariant mass m µµ be larger than 10 GeV. This requirement, applied in all three dimuon categories, also suppresses dimuons from decays of promptly produced low-mass SM resonances. However, lowp T muons can appear as muons with higher p T , with straighter tracks, when reconstructed from a small number of measurements. This gives rise to dimuons with an overestimated m µµ (above the 10 GeV threshold) and a mistakenly formed displaced vertex. To suppress such events, we reject STA-STA dimuons whose separation in η is small (|∆η µµ | < 0.1) if one of the muons is reconstructed in the barrel from fewer than 25 DT hits or if N(dimuon segments) < 6. In the other two dimuon categories, this background is suppressed by requiring that the hits associated with each TMS muon populate a minimum number of tracker layers, N(tracker layers). This minimum number depends on L xy and decreases from 6 for L xy < 15 cm to 3 for L xy > 45 cm.
Another source of SM background events is prompt high-mass dimuons that are reconstructed as displaced due to instrumental or reconstruction failures. Such dimuons mostly arise from DY dimuon production; contributions from processes such as tt and diboson production are relatively small. Events from DY ττ production with both τ leptons decaying to muons lead to a background with characteristics similar to those of the mismeasured DY µµ events. An important discriminating variable between signal and these backgrounds is the transverse collinearity angle ∆Φ between ⃗ p T µµ of the dimuon system and ⃗ L xy . When a pair of muons is produced in the decay of an LLP originating at the PV, both the resulting ⃗ p T µµ and ⃗ L xy point away from the PV, and |∆Φ| is small. On the other hand, mismeasured prompt dimuons, in particular those arising from the DY process, are expected to have a |∆Φ| distribution symmetric around π/2, because the directions of ⃗ p T µµ and ⃗ L xy are independent. We require |∆Φ| < π/4 in all three dimuon categories. The requirement is kept loose in order to preserve sensitivity to signal models featuring LLPs that decay into a dimuon and a particle (or particles) escaping detection, such as neutrinos, dark matter particles, or lightest supersymmetric particles [31,32]. We use the symmetric region |∆Φ| > 3π/4 as a control region for evaluating the contribution from DY and other prompt backgrounds, and the regions with π/4 < |∆Φ| < π/2 and π/2 < |∆Φ| < 3π/4 for validating background predictions.
The last important source of SM backgrounds is dijet and multijet events yielding dimuons that are formed from particles arising from different jets, either genuine muons or particles misidentified as muons. Such events contribute mainly to the TMS-TMS and STA-TMS categories and are strongly suppressed by the isolation criteria applied to TMS muons, as well as by the requirement that at least one of the muons that form the TMS-TMS dimuon has p T > 25 GeV. Furthermore, such events are expected to result in dimuons with both opposite-sign (OS) and same-sign (SS) electric charges, each with a symmetric |∆Φ| distribution. Therefore, in addition to binning in |∆Φ|, we classify selected dimuons as OS or SS, based on the observed muon charges. The signal selection requires that dimuons be OS, while SS dimuons constitute a control region used to evaluate backgrounds arising from dijet and multijet events, and from b hadron cascade decays.
A class of background events peculiar to the STA-TMS category consists of dimuons with their common vertex reconstructed on the wrong side of the PV relative to the direction of the TMS momentum vector. Since the TMS muon is usually reconstructed much more precisely than the STA muon, the reconstructed CV in STA-TMS events is usually located at or near the trajectory of the TMS muon. With ϕ TMS µ defined as the angle between the TMS muon ⃗ p T and ⃗ L xy , the location of the CV along the TMS muon trajectory depends on the STA muon and can be situated on either the correct side (|ϕ TMS µ | < π/2) or the wrong side (|ϕ TMS µ | > π/2) of the PV. We reject events on the extreme wrong side (|ϕ TMS µ | ∼ π) by requiring |ϕ TMS µ | < 2.9.
Next, we require that the dimuons be displaced with respect to the primary vertex. This is achieved by imposing requirements on the L xy significance L xy /σ L xy and the muon d 0 significance d 0 /σ d 0 . The L xy uncertainty σ L xy is calculated by combining uncertainties in the transverse positions of the CV and the PV; the d 0 uncertainty σ d 0 includes both track and PV uncertainties. The optimal values of requirements on L xy /σ L xy and d 0 /σ d 0 are chosen separately, for each of the dimuon categories, by maximizing an approximate figure of merit for the expected statistical significance of a discovery [33]. The expected number of signal events used for the optimization is obtained from simulation, while the number of background events is estimated from events in control regions in data using the procedures described in the next section. The criteria applied to L xy /σ L xy and d 0 /σ d 0 are: • L xy /σ L xy > 6 in the STA-STA dimuon category; • L xy /σ L xy > 6 and d 0 /σ d 0 > 6 in the TMS-TMS dimuon category; and • L xy /σ L xy > 3 and d 0 /σ d 0 > 6 in the STA-TMS dimuon category.
In the STA-STA category, the analysis sensitivity does not increase appreciably for any requirement on the d 0 significance once a requirement on L xy /σ L xy is in place; therefore, we do not apply any selection on d 0 /σ d 0 . In the TMS-TMS category, a requirement on d 0 /σ d 0 is applied to both muons and, since the residual background has a falling d 0 /σ d 0 distribution, the signal region is divided into three bins in the minimum of the two d 0 /σ d 0 values, min(d 0 /σ d 0 ): 6-10, 10-20, and >20. In the STA-TMS category, a requirement on d 0 /σ d 0 is applied only to the TMS muon, with no selection on d 0 /σ d 0 of the STA muon. Since the L xy resolution in the STA-TMS category is of the order of a few cm, a requirement on L xy /σ L xy is kept loose to preserve good efficiency for events with 10 ≲ L gen xy ≲ 60 cm, many of which are found in this category. Table 1 summarizes the event, muon, and dimuon selection criteria used in the analysis.
Finally, to test for the existence of an LLP with a given mass, dimuons satisfying the selection criteria are required to have m µµ within a specified interval containing the probed LLP mass. The width of each interval is chosen according to the mass resolution and the expected background. The resulting intervals typically contain a large fraction (90-99%) of putative signal with the probed mass. Since the mass resolution in the TMS-TMS category is far superior to that in the other two categories (1-3% compared to 10-25%, for LLP masses between 20 and 350 GeV), the minimum width of mass intervals varies from 3 GeV in the TMS-TMS category to ≈20 GeV in the STA-STA category.
The signal efficiency is defined as the ratio of the number of simulated signal events in which at least one dimuon candidate of any type passes all selection criteria (including the trigger) to the total number of simulated signal events. The efficiencies are estimated for LLP lifetimes corresponding to mean proper decay lengths in the range of 50 µm-1 km by reweighting events in the available simulated signal samples. The maximum efficiency, which is attained in the heavy-scalar model with m(Φ) = 1 TeV, m(X) = 150 GeV, and cτ = 1 cm, is approximately 50% for events that have only one LLP (Z D or X) decaying to muons. The efficiency becomes significantly smaller at low masses or at longer and shorter lifetimes, mostly because of a lower trigger and geometric acceptance and insufficient displacement. For example, the efficiency at the same set of masses decreases to approximately 15% at cτ = 100 cm, while the efficiency at cτ = 100 cm drops to below 5% when m(Φ) = 125 GeV and m(X) = 20 GeV. In accordance with Fig. 3, the vast majority of the signal with the shortest (longest) lifetimes is found in the TMS-TMS (STA-STA) dimuon category, whereas all three categories contain important fractions of the signal with the intermediate lifetimes. The efficiency in 2018 is higher than that in 2016, thanks mostly to the improved trigger and its relaxed requirements. The gain in efficiency is especially large, about a factor of two, at m(Φ) = 125 GeV and long lifetimes.

Background estimation and associated systematic uncertainties
Since the background events passing the event selection criteria arise mostly from misreconstructed prompt muons and muons in jets, their yield cannot be reliably ascertained from simulation. Therefore, we evaluate the expected background using events in data. The control regions used to estimate contributions from various types of background processes are chosen by inverting one or more selection criteria in order to obtain a region populated mostly by a given type of background and containing a negligible contribution from the signal processes. The definitions of the control regions and the details of the background evaluation procedure <π/4 <π/4 <π/4 opposite-sign muons yes yes yes differ for different dimuon categories and are described in the rest of this section. To avoid potential bias in the event selection, the events passing the full selection (i.e., those in the signal region) were "blinded" until the last steps of the analysis.
The contribution from cosmic ray muons is evaluated separately for each dimuon category from the number of dimuons satisfying all selection criteria but failing the cos α requirements. The evaluation procedure makes use of the efficiency of the cos α requirements measured from a sample of cosmic ray muons collected during periods with no beam. In all dimuon categories, the residual background arising from cosmic ray muons is estimated to be smaller than 0.1 events in all mass intervals combined.

Estimation of Drell-Yan and other prompt backgrounds
In all three dimuon categories, the contribution from prompt misreconstructed dimuons, collectively referred to as DY-like events, is evaluated from events in the signal-free |∆Φ|-symmetric control region, |∆Φ| > 3π/4: where N i DY (OS; |∆Φ| < π/4) and N i DY (OS; |∆Φ| > 3π/4) are, respectively, the numbers of DY background events in the signal and its |∆Φ|-symmetric control region; R i DY is the transfer factor accounting for the residual asymmetry in the population of events in the two |∆Φ| regions and obtained from auxiliary measurements; and the index i denotes the dimuon category (STA-STA, STA-TMS, or TMS-TMS). The number of DY dimuons in the |∆Φ| > 3π/4 region is taken to be the total number of events in that region minus the expected contribution from other types of background events estimated as discussed in Section 5.2.
The symmetry of the |∆Φ| distributions in this class of background events is studied using data and simulated events. In the STA-STA and STA-TMS categories, we use events in the control regions obtained by reversing the STA-to-TMS association. Specifically, we select events that consist of STA-STA or, alternatively, STA-TMS dimuons passing all selection criteria, but in which each of the constituent STA muons is associated with a TMS muon. To ensure that such STA-STA and STA-TMS dimuons are promptly produced (and thus are not signal), we require that the associated TMS-TMS dimuons, which have a far superior spatial resolution, are prompt. This is achieved by requiring L xy /σ L xy < 1.0 for the associated TMS-TMS dimuon in the STA-STA category, and d 0 /σ d 0 < 1.5 for the TMS muon associated with the STA muon in the STA-TMS category. To minimize contamination from muons from jets, which we discuss separately in what follows, each TMS muon in the associated TMS-TMS dimuon is required to satisfy the isolation requirement I rel trk < 0.05. Since the TMS and STA muons are predominantly reconstructed from information in different detectors (the tracker and the muon system, respectively), a genuine prompt muon giving rise to a displaced STA muon is usually accurately reconstructed as prompt by the TMS muon reconstruction. As a result, the aforementioned control regions contain genuine prompt dimuons that are reconstructed as displaced STA-STA or STA-TMS dimuons because of reconstruction failures or vertex fit anomalies in these categories, i.e., exactly the type of background events that we wish to study. The |∆Φ| distributions of STA-STA and STA-TMS dimuons in these control regions, in 2018 data and simulated background samples, are shown in Fig. 4. (The distributions in 2016 data are very similar.) The observed distributions, sculpted by the interplay between the geometric effects and the trigger and offline selection requirements, are well reproduced by the simulation. The distributions are still fairly symmetric around π/2, but there is a small asymmetry caused by the event selection criteria. Corrections accounting for this asymmetry are obtained from the ratio of events with |∆Φ| < π/4 and |∆Φ| > 3π/4 in the aforementioned control regions,  In the TMS-TMS dimuon category, the symmetry of the |∆Φ| distribution in DY-like backgrounds is assessed from events in the control region obtained by reversal of the requirement on χ 2 vtx . We observe a strong correlation between vertex χ 2 and DCA in both data and simulated DY events, which suggests that χ 2 vtx is effectively a measure of the distance between the TMS muons forming the dimuon. Further studies of TMS-TMS dimuons in DY events passing and failing the χ 2 vtx requirement confirm that they differ only in how far away the two muons are reconstructed from each other, and have very similar properties otherwise.
We use events in the inverted vertex χ 2 control region to evaluate the transfer factor R TMS-TMS

Estimation of nonprompt backgrounds
The background evaluation method described above is based on the symmetry of the |∆Φ| distribution and does not account for the contributions from background sources that yield dimuons exclusively or predominantly at small |∆Φ|. Such background sources include: dimuon decays of nonprompt low-mass resonances such as a J/ψ meson from b hadron decay; cascade decays of b hadrons; and dimuons formed from a pair of unrelated nonprompt muons in the same jet. If well reconstructed, most such background events have m µµ not exceeding a few GeV and are rejected by the m µµ > 10 GeV requirement. However, a small fraction of them with mismeasured m µµ can satisfy this requirement and pass the event selection. Such dimuons mostly have small |∆Φ| values (because the p µµ T and L xy vectors are collinear) and may have large, signal-like L xy /σ L xy and d 0 /σ d 0 values. They are also likely to have invariant masses close to the 10 GeV threshold. The other source of nonprompt background consists of dimuons formed from muons embedded in different jets. Since all these nonprompt background events arise from jets produced through the strong interaction, we collectively refer to them as quantum chromodynamics (QCD) events.
To gain insight into a contribution from this class of background events to the STA-STA and STA-TMS categories, we study events in control regions similar to those described above, but tailored to select events with muons embedded in jets. Once again, we invert the STA-to-TMS association and select STA-STA or, alternatively, STA-TMS dimuons passing all selection criteria, except that at least one of the constituent STA muons is associated with a TMS muon. To suppress |∆Φ|-symmetric background events such as DY as well as potential contributions from signal processes, we require each TMS muon to be nonisolated, defined as I rel trk > 0.1 in the STA-STA and I rel trk > 0.125 in the STA-TMS category. According to the simulation, this requirement selects a subset of events almost entirely composed of QCD events. Figure 5 (left) shows the |∆Φ| distribution of OS STA-STA dimuons in 2018 data in the samples thus obtained. Unlike the DY events in Fig. 4, which are approximately symmetric around π/2, the QCD events have a signal-like peak at |∆Φ| = 0. Most of these events are genuine low-mass dimuons that are reconstructed at higher m µµ because of poor p T resolution of STA muons, and hence pass the m µµ > 10 GeV requirement. This is demonstrated by Fig. 5 (right), which shows the distribution of well-measured m µµ of TMS-TMS dimuons associated with STA-STA dimuons with m µµ > 10 GeV.
Many of the background processes yielding small-|∆Φ| OS dimuons also give rise to small-|∆Φ| SS dimuons, either because these processes are charge symmetric or via the muon charge misassignment. Thus, we evaluate the contribution from the QCD background to the signal region, N i QCD (OS; |∆Φ| < π/4), from the number of small-|∆Φ| SS dimuons, N i (SS; |∆Φ| < π/4): The transfer factor R i QCD between the numbers of QCD events in these two regions is obtained from the ratio of OS to SS dimuons in the aforementioned control region with the STA-to-TMS association reversed and TMS muons not isolated: Because the charge misassignment and the probability of finding the STA-to-TMS association were found to be anticorrelated, the measurement of R i QCD in the STA-STA category is performed using dimuons with exactly one of the constituent STA muons associated with a TMS  A priori, we do not expect a large contribution from |∆Φ|-asymmetric low-mass dimuons in the TMS-TMS category, because of a far superior dimuon invariant mass resolution. Indeed, the study of simulated QCD events shows that a vast majority of both OS and SS TMS-TMS dimuons passing all selection criteria arise from a pair of unrelated nonprompt muons in two different jets. Such events do not contain genuine displaced dimuons and are expected to have a symmetric |∆Φ| distribution. Nevertheless, since some contribution from |∆Φ|-asymmetric dimuons may still be present in the background events in data, we prefer not to rely on the |∆Φ| symmetry in the evaluation of nonprompt backgrounds. Instead, similarly to the STA-STA and STA-TMS categories, we use the fact that dijet and multijet events give rise to both OS and SS dimuons, and base our estimate of the QCD background on the number of SS dimuons following Eq. (3).
The transfer factor R TMS-TMS QCD is obtained from the ratio of OS to SS dimuons in the control region with the muon isolation requirement reversed, which comprises dimuons passing the nominal event selection but with at least one muon with I rel trk > 0.075 and both with I rel trk < 0.5. We have verified that these events, as well as SS dimuons passing isolation requirements, contain negligible contributions from signal and DY events. As the signal region is divided into several min(d To avoid overestimating the DY background, the same QCD background evaluation method is applied to dimuons in the |∆Φ| > 3π/4 control region. The obtained estimate of the QCD background is then subtracted from the total observed number of OS dimuons with |∆Φ| > 3π/4 to obtain the number of DY dimuons in this |∆Φ| region, N i DY (OS; |∆Φ| > 3π/4), used for the evaluation of the DY backgrounds in the signal |∆Φ| < π/4 region according to Eq. (1). This procedure is not applied in the STA-STA category, where the |∆Φ|-symmetric QCD background is negligible. The sum of the QCD and DY background estimates constitute the total predicted background in the signal region. According to the background evaluation method, the DY backgrounds are expected to dominate at small d 0 /σ d 0 and L xy /σ L xy values, whereas the relative QCD contribution becomes larger as d 0 /σ d 0 and L xy /σ L xy increase. The uncertainty in the background predictions is dominated by the statistical uncertainty in the numbers of events in the |∆Φ| > 3π/4 and SS control regions.

Validation of background predictions
The background evaluation method described in Sections 5.1 and 5.2 is tested in several validation regions (VRs) that are expected to contain negligible contribution from signal. The evaluation of DY backgrounds is examined in the VRs obtained by inverting the L xy /σ L xy and d 0 /σ d 0 requirements and thereby enriched in this class of events. One example of such studies is shown in Fig. 6, which compares the background predictions to the observed distributions in the L xy /σ L xy < 6 VR in the STA-STA category. The yields in data are consistent with predictions of the method, which also correctly predicts a larger STA-STA background in 2016 compared to 2018 due to a lower tracking efficiency in a part of 2016 data [34].
In another check, we apply the background evaluation procedure to the TMS-TMS dimuons in the 2 < min(d 0 /σ d 0 ) < 6 sideband. The comparison of the predicted background and data in bins of L xy /σ L xy is shown in Fig. 7. The expected and observed numbers of events are in agreement in the entire probed L xy /σ L xy range. There are more background events in 2018 data than in 2016 data because of looser trigger requirements and larger integrated luminosity.
The evaluation of the |∆Φ|-asymmetric component of QCD backgrounds, which is particularly important in the STA-STA category, is tested in the low-mass (6 < m µµ < 10 GeV) VR, as well as in the region obtained by inverting the requirement on the minimum number of DT hits and muon segments applied to dimuons with |∆η µµ | < 0.1, referred to as the small-|∆η µµ | VR. Using dimuons with STA muons associated with TMS muons and taking well-measured m µµ and |∆Φ| values of corresponding TMS-TMS dimuons as proxies for true values of these quantities, we have verified that the samples of events in these VRs predominantly consist of small-|∆Φ| dimuons with mismeasured m µµ . Figure 8 shows the comparison of the predicted background yields and data in these two VRs. The low-mass VR is only available in 2018 data because the trigger used to collect 2016 data for this analysis included the m µµ > 10 GeV   requirement. The small-|∆η µµ | VR is available in both data sets, but since the number of events in this VR in 2016 data is small, an additional test is performed on a subset of 2018 data collected using the 2016 trigger, and therefore enriched in events similar to those recorded in 2016. The m µµ intervals of 10-32, 15-60, and 20-80 GeV shown for the small-|∆η µµ | VRs are the intervals chosen to probe LLP masses of 20, 30, and 50 GeV, respectively. The yields in data are found to be consistent with background predictions in all tests and m µµ intervals.  Finally, to ensure the validity of the method at different values of the main discriminating variable in the TMS-TMS and STA-TMS categories, the validation checks are performed in bins of d 0 /σ d 0 of the TMS muon. Such checks include comparisons in the d 0 /σ d 0 sideband (d 0 /σ d 0 < 6) in the signal |∆Φ| region, as well as those in the entire d 0 /σ d 0 range in the |∆Φ| sideband, π/4 < |∆Φ| < π/2. In the latter, the region with π/4 < |∆Φ| < π/2 is used as a signal-free proxy for the |∆Φ| < π/4 signal region. The background evaluation procedure is applied to the OS and SS dimuons in the |∆Φ|-symmetric region, π/2 < |∆Φ| < 3π/4, as well as SS dimuons with π/4 < |∆Φ| < π/2. The comparisons of the predicted background and 2018 data in the TMS-TMS and STA-TMS categories in this VR are shown in Fig. 9. The observed and expected numbers of events are consistent within statistical uncertainties.  Figure 9: Distributions in the π/4 < |∆Φ| < π/2 VR in 2018 data: (left) the smaller of the two d 0 /σ d 0 values for the TMS-TMS dimuon; (right) d 0 /σ d 0 of the TMS muon in the STA-TMS dimuon. The observed distributions (black points with error bars) are compared to the results of the background prediction method applied to events with π/2 < |∆Φ| < 3π/4. The stacked histograms show the expected numbers of DY (green) and QCD (yellow) background events. The last bin includes events in the overflow. The lower panels show the ratios of the observed to predicted numbers of events. The shaded area shows the statistical uncertainty in the background prediction.

Systematic uncertainties affecting the signal
Most of the systematic uncertainties affecting the signal efficiencies are evaluated separately in each dimuon category and for each data-taking year. Unless stated otherwise, we consider sources of uncertainties to be uncorrelated among different dimuon categories and years.
In the STA-STA and STA-TMS categories, the dominant systematic uncertainties come from the STA muon identification and trigger efficiencies. At small displacements, both efficiencies are accurately measured as a function of muon p T and η by applying the "tag-and-probe method" [28] to muons from J/ψ meson and Z boson decays. The differences in the identification and trigger efficiencies between data and simulation are used to correct the signal simulation yields. In the STA-STA category, these corrections range from 0.78 to 1.13, depending on the signal sample. The evolution of efficiencies with displacement is studied using a sample of cosmic ray muons collected during periods with no beam, and additional d 0 -dependent corrections and systematic uncertainties are derived. At d 0 = 10 (100) cm, the correction amounts to 0.99 (0.96) per muon, whereas the uncertainty is on the order of 10 (35)%. Since the d 0dependent uncertainty is dominated by the accuracy of the L1 trigger efficiency measurements, it is taken to be correlated among all dimuon categories.
The dominant systematic uncertainties in the TMS-TMS category come from the d 0 dependence of the L1 trigger efficiency discussed above and the efficiency to reconstruct displaced muons in the tracker. The evolution of the tracking efficiency with d 0 is measured using a sample of cosmic ray muons and compared to the tracking efficiency predicted by simulation. Based on the results of the comparison, we assign a 5% systematic uncertainty per muon for muons with d 0 > 1 cm. The overall efficiency corrections applied to the simulated signal yields range from 0.74 to 1.08, depending on the signal sample, and arise mostly from imperfect modeling of the HLT efficiencies at small displacements.
The remaining systematic uncertainties related to the signal efficiency are much smaller. The impact of mismodeling of the muon p T resolution on the signal yield is evaluated by smearing the muon p T in simulated signal events according to the measurements performed using cosmic ray muons and muons from Z boson decays. This leads to variations that are less than 2% at all signal masses except for m(Z D ) = 10 GeV. Corrections of up to 2% are applied to the TMS muon efficiency to account for the difference in efficiency of isolation requirements measured using muons from Z boson decays, and an additional systematic uncertainty of 2% is assigned. A systematic uncertainty ranging from 1 to 8%, depending on the signal sample, is assigned to account for mismodeling of the DCA requirement in the STA-STA category. The efficiency of the vertex χ 2 requirement as a function of displacement is studied using cosmic ray muons and muons from Z boson decays in the STA-STA category, and muons from decays of nonprompt J/ψ mesons in the TMS-TMS category. The differences between data and simulation contribute an uncertainty of 2% in each category. The efficiencies of several other selection criteria, such as requirements on the number of tracker hits upstream of the vertex position and the difference between the number of pixel hits on two TMS muons, are found to be well modeled by simulation, and no additional uncertainty is assigned.
The uncertainty in the integrated luminosity, partially correlated between the years, is 1.2% in 2016 [15] and 2.5% in 2018 [16]. The uncertainty in the signal efficiency due to pileup is 2%. Both uncertainties are correlated among dimuon categories.

Results
The predicted background yields in the representative m µµ intervals and the corresponding numbers of observed events are shown in Fig. 10 for the STA-STA category and Fig. 11 for the STA-TMS category. For illustration, signals at the level of the median expected exclusion limits at 95% confidence level (CL) in the absence of signal are also shown. As expected, events observed in the STA-STA and STA-TMS categories are predominantly at low masses-14 out of 18 STA-STA and 9 out of 13 STA-TMS events have m µµ < 20 GeV-and have characteristics typical of those for QCD background events. The numbers of observed events and the predicted background and signal yields in the TMS-TMS category are shown in Fig. 12 as functions of m µµ in each of the three min(d 0 /σ d 0 ) bins and in Fig. 13 as a function of min(d 0 /σ d 0 ). The observed TMS-TMS events have a steeply falling min(d 0 /σ d 0 ) distribution and cluster at m µµ values of a few tens of GeV, which are both consistent with the characteristics of the expected background. The numbers of observed events are consistent with the predicted background yields in all dimuon categories and m µµ intervals, in both data sets. No significant excess of events above the SM background is observed.
For each of the two benchmark models, we compute upper limits on the product of the signal production cross section σ and the branching fraction B to two muons as a function of mass and mean proper decay length. The limit extraction is based on a modified frequentist approach [35,36] and uses the CMS COMBINE package developed for statistically combining the results of Higgs boson searches [37]. The method yielding background predictions in the signal region is implemented using a multibin likelihood, which is a product of Poisson distributions corresponding to the signal region and the control regions. The systematic uncertainties affecting the signal yield are incorporated as nuisance parameters using log-normal distributions. The expected and observed upper limits are evaluated through the use of simulated pseudoexperiments. For each signal model, the limits are first computed separately in each dimuon category and for each data-taking year. The individual likelihoods are then combined to obtain the combined limits.
The signal efficiencies are obtained from simulation and further corrected by the data-to-simulation scale factors described in Section 6; they are computed separately for each year, signal    model, dimuon category, and mass interval. A reweighting procedure is employed to calculate an estimated number of signal events for lifetimes other than the lifetimes used to produce the available simulated signal samples. Given the smallness of the expected background and taking into account the selection efficiencies discussed in Section 4, an introduction of a separate category for events with two dimuons does not increase the sensitivity of the analysis significantly even in the most favorable case for the 4µ signal events, namely B(X → µµ) = 1. The gain is negligible for smaller B(X → µµ) values such as those predicted by the HAHM model. Therefore, no distinction is made between events with one and two dimuons of the same type. Events with two TMS-TMS dimuons are assigned to the min(d 0 /σ d 0 ) bin encompassing the larger of the two min(d 0 /σ d 0 ) values.
The observed and expected 95% CL upper limits on the product σ(Φ → XX)B(X → µµ) as a function of cτ(X) in the simplified heavy-scalar model for m(Φ) = 125, 200, 400, and 1000 GeV are shown in Figs. 14, 15, 16, and 17, respectively. Each figure shows the limits obtained with the ensemble of 2016 and 2018 data in the individual dimuon categories, as well as their combination. The search is sensitive to a broad range of cτ from 30 µm to more than 1 km. As expected, the three dimuon categories reach maximum sensitivity at different cτ values and all give relevant contributions to the overall sensitivity of the search. The limits are most restrictive for cτ between 0.1 mm and 10-100 m, excluding σ(Φ → XX)B(X → µµ) smaller than 1 fb, and become more stringent at high LLP masses, owing to a higher signal efficiency and lower background. The smallest σ(Φ → XX)B(X → µµ) value excluded is 0.05 fb, at m(X) = 350 GeV and cτ between 0.3 and 30 cm. compared to the theoretical predictions for a set of representative B(H → Z D Z D ) values ranging from 1% to 0.001%. In the m(Z D ) range of 20-60 GeV, the branching fraction B(H → Z D Z D ) of 1% is excluded in the cτ(Z D ) range from a few tens of µm to approximately 100 m, whereas B(H → Z D Z D ) as low as 0.01% is excluded in the range of 1 mm to 1 m. These constraints on rare SM Higgs boson decays are tighter than those derived from searches for invisible Higgs boson decays [38] and from indirect constraints from measurements of the SM Higgs boson couplings [39]. At m(Z D ) > 20 GeV, the limits obtained are the best to date for all cτ(Z D ) values except those between ≈0.5 and 500 mm (depending on m(Z D )), where our search is complemented by the CMS search using data collected with a dedicated high-rate data stream [40].
The observed 95% CL exclusion contours in the (m(Z D ), cτ(Z D )) plane, determined from the observed 2016-18 upper limits on cτ(Z D ) for several representative values of B(H → Z D Z D ), are shown in Fig. 19 (left). These results can be translated into limits on the kinetic mixing ϵ following the relationship between the dark photon mass, lifetime, and ϵ. The resulting 95% CL exclusion contours in the (m(Z D ), ϵ) plane are shown in Fig. 19 (right). Our analysis excludes a wide range of ϵ values, between 9 × 10 −9 and 6 × 10 −6 at m(Z D ) = 10 GeV and between

Summary
Data collected by the CMS experiment in proton-proton collisions at √ s = 13 TeV in 2016 and 2018 and corresponding to an integrated luminosity of 97.6 fb −1 have been used to conduct an inclusive search for long-lived exotic neutral particles (LLPs) decaying to a pair of oppositely charged muons. The search is largely model-independent and is sensitive to a broad range of LLP lifetimes and masses. No significant excess of events above the standard model background is observed. The results are interpreted as limits on the parameters of the hidden Abelian Higgs model, in which the Higgs boson decays to a pair of long-lived dark photons Z D , and of a simplified model, in which LLPs are produced in decays of an exotic heavy neutral scalar boson. In the mass range 20 < m(Z D ) < 60 GeV, a branching fraction of the Higgs boson to dark photons of 1% is excluded at 95% confidence level in the range of proper decay length cτ(Z D ) from a few tens of µm to ≈100 m. The results of this search significantly extend the previously excluded range of model parameters. For the hidden Abelian Higgs model with m(Z D ) greater than 20 GeV and less than half the mass of the Higgs boson, they provide the best limits to date on the branching fraction of the Higgs boson to dark photons for cτ(Z D ) (varying with m(Z D )) between 0.03 and ≈0.5 mm, and above ≈0.5 m. At exotic scalar boson masses larger than the Higgs boson mass, our results represent the best current constraints for all considered LLP masses and lifetimes.   [8] CMS Collaboration, "Search for long-lived particles that decay into final states containing two muons, reconstructed using only the CMS muon chambers", CMS Physics Analysis Summary CMS-PAS-EXO-14-012, 2015.