Search for long-lived particles decaying into muon pairs in proton-proton collisions at $\sqrt{s} =$ 13 TeV collected with a dedicated high-rate data stream

A search for long-lived particles decaying into muon pairs is performed using proton-proton collisions at a center-of-mass energy of 13 TeV, collected by the CMS experiment at the LHC in 2017 and 2018, corresponding to an integrated luminosity of 101 fb$^{-1}$. The data sets used in this search were collected with a dedicated dimuon trigger stream with low transverse momentum thresholds, recorded at high rate by retaining a reduced amount of information, in order to explore otherwise inaccessible phase space at low dimuon mass and nonzero displacement from the primary interaction vertex. No significant excess of events beyond the standard model expectation is found. Upper limits on branching fractions at 95% confidence level are set on a wide range of mass and lifetime hypotheses in beyond the standard model frameworks with the Higgs boson decaying into a pair of long-lived dark photons, or with a long-lived scalar resonance arising from a decay of a b hadron. The limits are the most stringent to date for substantial regions of the parameter space. These results can be also used to constrain models of displaced dimuons that are not explicitly considered in this paper.

Other scenarios may produce a low-mass, long-lived resonance decaying into a muon pair. For instance, one of the minimal extensions to the SM adds a singlet scalar field φ, which mixes with the H boson and couples to all SM fermions [15,16]. In the hypothesis of weak coupling to SM fermions with a signal mixing angle θ 10 −3 , φ will be long-lived. Such a scalar resonance may be produced in the decay of a b hadron, h b → φX, as illustrated in Fig. 2. Constraints on this model have been previously placed by the CHARM [17] and LHCb [18,19] experiments. dedicated dimuon trigger stream (scouting) with low transverse momentum (p T ) thresholds, recorded at high rate by retaining a reduced amount of information. The rate of scouting triggers is higher than that of the standard triggers owing to less stringent requirements that enable dimuon resonance searches across mass and lifetime ranges that are otherwise inaccessible. The selected data correspond to a total integrated luminosity of 101 fb −1 (41.5 fb −1 collected in 2017 and 59.7 fb −1 collected in 2018). Additional details on the data sets and triggers are provided in Sections 2 and 3. The search targets narrow, low-mass, long-lived resonances decaying into a pair of oppositely charged muons, where the lifetime of the long-lived particle is such that the transverse displacement (l xy ) of its decay vertex is within 11 cm of the primary interaction vertex (PV). The signal is expected to appear as a narrow peak on the dimuon mass continuum, with a resonance width smaller than the experimental mass resolution.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The pixel tracker consists of four concentric barrel layers at radii of 29, 68, 109, and 160 mm, and three disks on each end at distances of 291, 396, and 516 mm from the center of the detector [20]. 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. [21].
A two-level trigger system is used to select events of potential physics interest. The first level of the CMS trigger system (L1) [22], composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events of interest in a time interval of less than 4 µs. The high-level trigger (HLT) processor farm further decreases the event rate from around 100 kHz to about 1 kHz, before data storage. A more detailed description of the CMS trigger system can be found in Ref. [23].
In addition to the standard HLT selection streams, a dedicated set of triggers exists that allows the exploration of otherwise inaccessible phase space. This approach is referred to as data scouting [14,24]. The scouting trigger algorithms used in this search select events containing muon pairs with a mass m µµ 2m µ at a rate of about 3 kHz. In order to compensate for the large rate, the event size is reduced by up to a factor of 1000 by only retaining limited information, as reconstructed at the HLT.

Data sets and triggers
The standard trigger streams typically select events containing at least one particle (e.g., a lepton) or jet with large p T . Thus, such triggers are not necessarily ideal in a generic search for light narrow resonances decaying into two muons, with m µµ 45 GeV [14]. The scouting triggers used in this search compensate for the large trigger rate by recording only limited event information reconstructed at the HLT. For muons, which are reconstructed by combining the information from the tracker and the muon detectors, the recorded information includes for each muon its four-momentum, the number of hits left by each muon track in the tracker, the normalized χ 2 of the muon track, and the isolation of the muon, computed as the scalar sum of the p T of all additional tracks in a cone of size ∆R = √ (∆η) 2 + (∆φ) 2 = 0.3 around the muon, where φ is the azimuthal angle with respect to the counterclockwise beam axis. Muon tracks are required to leave a hit in at least two layers or disks of the pixel tracker and are used in pairs to form dimuon vertices, considering all possible pairs. In the following, these vertices are termed secondary vertices (SVs), and they may be displaced from the PV or not. Candidate proton-proton (pp) interaction vertices are reconstructed using any track with a hit in each layer or disk of the pixel tracker, and the candidate vertex with the largest value of summed track p 2 T is taken to be the PV of the pp interaction. Events in the scouting data sets contain information on the normalized χ 2 of all reconstructed vertices, together with their position with respect to the nominal interaction point. The four-momenta of jets formed by clustering calorimeter deposits using the anti-k T algorithm [25,26] with a distance parameter of 0.4 are also available, for jets with p T > 20 GeV and |η| < 3.0.
Events are preselected at the L1 trigger by requiring data collected in 2017 and 2018 to pass at least one of three criteria: • the presence of two L1 muons with opposite charges (OS), p µ T > 4 GeV (for 2017 data) or 4.5 GeV (for 2018 data) and ∆R(µ 1 , µ 2 ) < 1.2; • the presence of two OS L1 muons with |η µ | < 1.4 and ∆R(µ 1 , µ 2 ) < 1.4; • the presence of two L1 muons with p µ 1 T > 15 and p µ 2 T > 7 GeV. Then, only the preselected events containing two OS muons with p µ T > 3 GeV and |η µ | < 2.4 at the HLT are retained.
The efficiency of the scouting triggers is measured in data, as a function of the dimuon l xy and of the minimum p µ T of the pair, using events with at least two muons selected with independent standard triggers. The measured efficiency and its uncertainty are accounted for in the normalization of simulated events, as discussed in Section 4.
The resolution in p µ T for muons reconstructed at the HLT is measured to be roughly 10% worse relative to muons undergoing offline reconstruction, because of differences in the calibration of the HLT and the offline reconstruction. As only a reduced amount of information is available in scouting data, no additional calibration is performed for the muons used in this search.

Monte Carlo simulation
The SM background is estimated directly from data as a continuum background in the m µµ spectrum parametrized by analytical functions. Thus, the search does not rely on Monte Carlo (MC) simulation in order to estimate the background.
On the other hand, MC simulations of beyond the SM (BSM) signals are used to interpret the results of the search. Two signal models are simulated, as described in Section 1: • A model of H boson production via the dominant mode of gluon-gluon fusion, with H → Z D Z D , where Z D is a long-lived resonance and at least one Z D decays to a pair of oppositely charged muons. This model is generated using the POWHEG v2.0 [27][28][29] generator at next-to-leading order (NLO) precision and the JHUGEN 7.4.0 [30][31][32][33][34] generator at leading order (LO) precision. The JHUGEN generator is used to model the H boson decay, setting the properties (mass and lifetime) of the Z D boson. The generators are interfaced with PYTHIA 8.230 [35] for the simulation of parton showering and hadronization. The simulated signal is normalized to a Higgs production cross section calculated at next-to-next-to-NLO (N 3 LO) in the strong coupling α S [36]. The decays H → ZZ D and H → Z D Z D through the hypercharge portal have negligible branching fractions for parameter values explored in this analysis.
• A model in which a scalar resonance φ is emitted in the decay of b hadrons, where φ is long lived and decays to a pair of oppositely charged muons with 100% probability. The production of b hadrons is initially generated at LO precision using PYTHIA 8.2. Subsequently, a fixed-order calculation at next-to-leading log (FONLL) in α S [37][38][39][40][41] is used to reweight the p T spectrum of b hadrons decaying into φ. The simulated signal is also normalized to the FONLL calculation.
The generators use the CP5 PYTHIA tune [42] and the NNPDF3.1 sets of parton distribution functions (PDFs) [43], while the detector response is simulated with a GEANT4 model [44] of the CMS detector. The products from multiple pp interactions are superimposed on those from the hard collision (pileup), and the simulated samples are reweighted such that the number of collisions per bunch crossing reflects the distribution observed in data.
The search explores a dimuon mass range from about 200 MeV (i.e., m µµ 2m µ ) to about 50 GeV, and a transverse displacement range 0 < l xy < 11 cm from the PV. The hard cutoff in l xy at 11 cm is due to the definition of the scouting trigger stream, which requires the presence of hits in at least two layers of the pixel tracker. Signal samples are generated to cover the full ranges of interest, in both the m µµ and the cτ 0 (mean proper decay length) dimensions. Only masses m φ 5 GeV can be probed in the model where a scalar resonance (φ) is produced on shell in the decay of a b hadron.

Event selection
The events collected using the scouting triggers (see Section 3) are further selected to contain at least one pair of OS muons associated with an SV that satisfies a number of quality criteria: • l xy < 11 cm; • x coordinate uncertainty < 0.05 cm; • y coordinate uncertainty < 0.05 cm; • z coordinate uncertainty < 0.10 cm; • χ 2 /dof SV < 5, where dof is the number of degrees of freedom.
Muons are also required to satisfy additional quality criteria: the χ 2 /dof track of each muon track is required to be less than 3, and the track itself is required to cross at least 6 layers (pixel or strip) of the tracker.
A large fraction of events contain only one OS muon pair and an associated SV. Among those events with more than one pair of OS muons, the pairs are ordered according to increasing χ 2 SV and only the first or the first two pairs are used in this search. A few of the selection criteria applied to the second muon pair are less stringent than those applied to the first pair in order to maximize the sensitivity to BSM signatures.
For the first (second) muon pair, a muon is considered isolated if the track isolation, as defined in Section 3, is less than 10 (20)% of p µ T , and if no jet is found within a cone of size ∆R < 0.3 around the muon itself.
As further described in Section 5.2, we explore isolated, partially isolated, and nonisolated dimuon topologies. Isolated topologies consist of events where both muons are isolated; partially isolated topologies consist of events where exactly one of the two muons fails the isolation criteria; finally, in nonisolated topologies both muons fail the isolation criteria. In events with two muon pairs, the only category considered requires all four muons to be isolated.
For a BSM signal, we expect the dimuon system transverse momentum vector, p µµ T , to be collinear with the SV vector, which is defined as the vector − → SV connecting the PV to the SV. This does not apply to a muon pair where the two muons originate from different vertices and erroneously form a common SV because of the accidental crossing of their trajectories. In order to suppress such background, we require the azimuthal angle between the two vectors (∆φ( p µµ T , − → SV)) to be less than 0.02 (0.10) for the first (second) muon pair. Similarly, the opening angle between two muons erroneously forming a common SV tends to be very large, independently of the boost of the dimuon system. Hence, in order to further reduce the background from events with the wrong muon SV association, we also require ∆φ(µ 1 , µ 2 ) < 2.8. The main goal of such selections is to suppress background where SVs are formed from cosmic ray muons, muons from pileup collisions, and muons from back-to-back quantum chromodynamics dijet events.
Muons from pileup vertices can overlap in the r-φ plane despite being apart in r-z at their origin, and, as a consequence, they can be wrongly associated with a common SV. These muon pairs are characterized by a large ∆η(µ 1 , µ 2 ) and a small ∆φ(µ 1 , µ 2 ). We exploit this feature to suppress this background contribution with the requirement log 10 (|∆η µµ |/|∆φ µµ |) < 1.25.
Additional background sources include "prompt" (i.e., nondisplaced) muons and material interaction vertices, i.e., vertices from particles produced in interactions with detector material, e.g., in the beampipe or the pixel layers. Muons from truly displaced vertices do not leave hits in the tracker layers between the interaction point and the SV. In order to suppress the background from prompt muons, we therefore require that the observed number of hits found in the pixel tracker for each muon be no larger than the number expected from the outward propagation of the muon itself from the associated SV. This requirement is applied only when l xy > 3.5 cm to ensure that the SV is located beyond the first pixel tracker layer. Additionally, material interaction vertices are suppressed by rejecting events with an SV located within 0.5 mm of a pixel tracker module.
In order to identify displaced vertices of interest, we also exploit the muon track transverse impact parameter, d xy . To this end, we select events based on the d xy significance (|d xy |/σ d xy ) and the lifetime-scaled d xy (|d xy |/(l xy m µµ /p µµ T )). The first expression normalizes d xy against its measurement uncertainty, while the second normalizes d xy against the lifetime corresponding to the observed SV transverse displacement. This allows the same selection requirements to be used for all lifetime hypotheses. For the first (second) muon pair, we require each muon to satisfy: The total selection efficiency after the trigger selections typically ranges from about 55 to 75% in simulated signal events, depending on the details of the models.

Event categorization
This search targets long-lived low-mass dimuon resonances that are characterized by displaced decays into pairs of muons. Therefore, the primary feature to be exploited in order to discriminate a potential signal from the SM background is the presence in each event of at least one pair of muons and an associated SV. Most of the physics processes contributing to the SM background are more highly suppressed at increased displacement of the SV from the PV. Because the vertex position resolution along the beam (z axis) is not as precise as the one achieved in the transverse (xy) plane, we only use the l xy of the SV to categorize dimuon events. The l xy categorization is intended to maximize the search sensitivity to a range of potential BSM signal lifetimes, and is based on the CMS pixel tracker geometry. We define the following six l xy categories: In each l xy category, we further divide dimuon events in two bins of p µµ T (< 25 and ≥ 25 GeV). This categorization improves the sensitivity to different signal topologies since p µµ T tends to be low in signal models where long-lived particles arise from the decay of a b hadron, while it is expected to be larger in other models. Finally, in each (l xy , p µµ T ) category, dimuon events are divided into bins of isolation (to distinguish isolated, partially isolated, and nonisolated topologies).
In each (l xy , p µµ T , isolation) bin, we define mass windows sliding along the dimuon invariant mass spectrum (see Section 6), and we perform a search for a resonant dimuon peak in each mass window. Mass regions corresponding to known resonances decaying either in dimuon pairs or in a pair of charged hadrons are not considered, i.e., they are "masked". As an example, in Fig. 3 (4) we show the m µµ distribution in bins of l xy as obtained from selected dimuon events where both muons are isolated and have p µµ T below (above) 25 GeV. In events with two muon pairs each associated with an SV, we further require the difference between the two dimuon masses to be within 5% of their mean, and the four-muon mass to be consistent with the mass of the SM Higgs boson (115 < m 4µ < 135 GeV). The selected fourmuon events are treated as an exclusive independent category, aimed at improving the search sensitivity to models of BSM physics where an SM Higgs boson decays to a pair of Z D bosons, each decaying to two muons. After applying all selection criteria, we observe zero events in  In each dimuon mass window, the SM background is modeled by fitting the dimuon invariant mass distribution in data, as further described in Section 6. Finally, the dimuon mass windows in all (l xy , p µµ T , isolation) bins and the single four-muon event category are used in combination as input to a final likelihood fit for the interpretation of the results, as described in Section 7. In the four-muon category, the two dimuon invariant mass distributions are not fitted separately.  = 100 mm, corresponding to = 1.2 × 10 −7 , is overlaid. The signal event yield corresponds to a value of the branching fraction product B(H → Z D Z D ) B(Z D → µµ) = 1.2 × 10 −5 , equal to the median expected exclusion limit at 95% CL set in this paper (see Section 7), and is multiplied by a factor of 100 for display purposes. The vertical gray bands indicate mass ranges containing known SM resonances, which are masked for the purpose of this search.

Parametrization of the invariant mass distribution
The search is performed by fitting mass distributions in the various categories to the sum of signal plus SM background models. We target low-mass dimuon resonances whose intrinsic width is assumed to be much narrower than the detector resolution. The signal dimuon mass resolution (σ µµ mass ) is estimated from simulated events to be about 1.1% of the mass, with an uncertainty of about 50%, and both are independent of the signal mass and lifetime hypotheses. The signal four-muon mass resolution is similarly estimated from simulated H → Z D Z D → 4µ events to be about 1.1% of the Higgs boson mass, with an uncertainty of about 50%. We use the sum of a Gaussian function and a double-sided Crystal Ball function [45,46] to model signal dimuon (four-muon) invariant mass distributions in dimuon (four-muon) events.
In each dimuon event category, the SM background is modeled in sliding windows of ±5σ µµ mass around the considered signal mass hypothesis. These windows are required to not overlap with the masked mass regions corresponding to known SM resonances (see Section 5.2). Instead, in the four-muon event category we consider a single window centered around the known Higgs boson mass (115 < m 4µ < 135 GeV). A number of functional forms are considered in order to model the background mass distribution. These include Bernstein polynomials, exponential functions, as well as combinations of the two. For each family of functional forms, the best order in each mass window and event category is selected by means of a Fisher test [47]. First, the lowest order (N) function is used to fit the data. Then, the next-order (N + 1) is used, and the difference 2∆NLL N+1 = 2 (NLL N+1 − NLL N ) (with NLL denoting the negative logarithm of the likelihood of the fit) is evaluated to determine whether the data support the need for a higher order function. This decision is based on the fact that 2∆NLL N+1 is asymptotically described by a χ 2 distribution with M degrees of freedom, where M is the difference in the number of free parameters in the (N + 1)th and Nth order functions. A p-value is calculated as P M χ 2 > 2∆NLL N+1 where P M χ 2 min is the χ 2 tail probability for M degrees of freedom. If the p-value is less than 0.05, the higher order function is retained, since it is determined to significantly improve the description of the data. Once the best order N for each family of functions has been determined, the corresponding functional forms are entered in an envelope then used in the discrete profiling method [48], where the choice of the background function is treated as a discrete nuisance parameter in the fit to account for the uncertainty associated with the arbitrary choice of the function. For Bernstein polynomials, we also include orders N − 1 (if the selected best order N > 0) and N + 1 in the list of suitable functions. Typically, the selected best order N is ≤ 4 and is smaller in event categories with a small number of observed events. In addition, we use a goodness-of-fit test to remove background models that do not describe the data appropriately. The test is based on a χ 2 test statistic, which is converted into a p-value through a set of pseudo-experiments. Models with p < 0.01 are not considered. The potential presence of a bias in the measurement of a signal due to the choice of the functional forms used to model the SM background is also assessed by means of pseudo-experiments. First, a varying amount of signal is injected on top of the background generated according to the selected functional form. Then, a background+signal fit is performed allowing the signal yield to float freely. The bias, quantified as the difference between the measured and injected signal yields relative to the statistical uncertainty in the measured signal yield, is found to be negligible.
The background-only fit results for the selected functional forms in one of the dimuon search bins are shown in Fig. 5 together with the dimuon invariant mass distribution expected for a representative signal.

Results and interpretation
Binned maximum likelihood fits to the data are performed simultaneously in all search bins, under either background-only or background+signal hypotheses, using background and signal models and uncertainties described in Section 6. Additional log-normal constraint terms are used to account for the uncertainties in the signal yields, when considered (Section 7.1).
No significant peak-like structures are observed in data. The background+signal fits are used to set 95% confidence level (CL) upper limits on the branching fractions for the signal models under consideration. Limits are set using a modified frequentist approach, employing the CL s criterion [49][50][51][52]. These limits are then used, in conjunction with the theoretical cross section calculations, to exclude ranges of masses and/or lifetimes for the BSM particles of the signal

Systematic uncertainties in the signal yields
A systematic uncertainty in the trigger efficiency (see Section 3) is assessed that ranges between 3 and 30%, depending on the p T and displacement of the muons. It is derived from studies of data collected with different triggers and accounts for differences between the measurement and the simulation. In order to account for potential mismodeling of the signal simulation, we also assess a systematic uncertainty in the signal yields equal to 100% of the inefficiency of the selection (see Section 5.1) applied after the trigger selections. The resulting systematic uncertainty in the signal yield is then between 25 and 45%, depending on the signal mass and lifetime. Additionally, we assess a systematic uncertainty in the expected signal yields arising from the uncertainty in the luminosity measurement of 2.3% [53] for 2017 and 2.5% [54] for 2018 data. Finally, we account for uncertainties resulting from the limited sizes of the simulated signal samples. Uncertainties arising from the choice of the PDF and of the renormalization and factorization scales used in the event generator are negligible compared to others. Uncertainties in the trigger efficiency and integrated luminosity are treated as correlated across search bins. Other uncertainties are taken as uncorrelated. Since we find no significant difference between 2017 and 2018 data and simulations, all uncertainties are treated as correlated across data taking periods.

Constraints on models of BSM physics
We set upper limits at 95% CL on the branching fractions for the h b → φX and H → Z D Z D signal models. For the h b → φX model, Fig. 6 shows the upper limits on the branching fraction product B(h b → φX) B(φ → µµ), as a function of the signal mass hypothesis, for a few representative lifetime hypotheses. For the H → Z D Z D model, the upper limits on the branching fraction product B(H → Z D Z D ) B(Z D → µµ), as a function of the signal mass, are shown in Fig. 7. In addition, Fig. 8 shows the upper limits on the branching fraction B(H → Z D Z D ), using values of B(Z D → µµ) from the model of Ref. [6]. The upper limits shown in Figs. 6 and 7 are obtained using only the combination of dimuon event categories. For the limits shown in Fig. 8, the four-muon event category is added to the combination, since in this case the relative acceptances for signal events in the dimuon and four-muon categories can be determined using the B(Z D → µµ) assumption.
Constraints for signals with mass above the b hadron mass are more stringent since the background from displaced muons from h b decays is completely eliminated. The constraints tend to get weaker at longer lifetime because of the loss of acceptance due to decays beyond the scouting trigger acceptance, which is determined by the radial coverage of the pixel detector. We note that our limits on the model with a long-lived scalar emitted in h b decays cannot be directly compared to similar limits from the LHCb experiment [18,19], since our limits are on the inclusive h b → φX branching fraction, while the LHCb limits are on the exclusive B 0 → φK * (892) 0 and B + → φK + branching fractions. However, we scale our inclusive limits based on the fraction of B 0 or B + production from MC and thus compare our exclusive branching fraction limits to those of LHCb. Such comparisons are provided as supplemental material [55].
Observed limits at 95% CL as contours in the parameter space of signal mass and signal lifetime for both the models considered in this search are shown in Fig. 9 and Fig. 10.
In Fig. 11, we show the excluded region at 95% CL for the inflaton model [15] in the parameter space containing the signal mass (m φ ) and the square of the signal mixing angle (θ 2 ) using the results from the h b → φX signal model. The signal mixing angle θ determines the inflaton's coupling to the SM fields via mixing with the Higgs boson. In Fig. 12, we show the observed limits at 95% CL as contours in the parameter space containing the signal mass (m Z D ) and the signal kinetic mixing coupling for the H → Z D Z D signal model.         : Region excluded at 95% CL using this search (shown in hatched gray) for the inflaton model [15] in the parameter space containing the signal mass (m φ ) and the square of the signal mixing angle (θ 2 ). The regions forbidden by theory and cosmological constraints are shown in green and yellow respectively. The vertical gray bands indicate mass ranges containing known SM resonances, which are masked for the purpose of this search. The limits are obtained using the combination of all dimuon event categories. Tabulated results are provided in the HEPData record for this analysis [55]. To facilitate the reinterpretation of these results in other BSM models, we also provide upper limits at 95% CL on the number of observed events as a function of the dimuon mass in a number of signal regions characterized by different l xy ranges, with and without an isolation requirement, inclusive in p µµ T as well as with p µµ T ≥ 25 GeV, together with parametrizations of the event selection efficiency.

Summary
A search for displaced dimuon resonances has been performed using proton-proton collisions at a center-of-mass energy of 13 TeV, collected by the CMS experiment at the LHC in 2017 and 2018, corresponding to an integrated luminosity of 101 fb −1 . The data sets used in this search are collected using a dedicated dimuon trigger stream with low transverse momentum thresholds, recorded at high rate by retaining a reduced amount of information, in order to explore otherwise inaccessible phase space at low dimuon mass and nonzero displacement from the primary interaction vertex. No significant excess beyond the standard model expectation is found, and the data are used to set constraints on a wide range of mass and lifetime hypotheses for models of physics beyond the standard model where a Higgs boson decays to a pair of long-lived dark photons, or where a long-lived scalar resonance arises from the decay of a b hadron. These constraints are the most stringent to date for substantial regions of the parameter space.