A search for pair production of new light bosons decaying into muons in proton-proton collisions at 13 TeV

A search for new light bosons decaying into muon pairs is presented using a data sample corresponding to an integrated luminosity of 35.9fb − 1 of proton-proton collisions at a center-of-mass energy √ s = 13 TeV, collected with the CMS detector at the CERN LHC. The search is model independent, only requiring the pair production of a new light boson and its subsequent decay to a pair of muons. No signiﬁcant deviation from the predicted background is observed. A model independent limit is set on the product of the production cross section times branching fraction to dimuons squared times acceptance as a function of new light boson mass. This limit varies between 0.15 and 0.39fb over a range of new light boson masses from 0.25 to 8.5GeV. It is then interpreted in the context of the next-to-minimal supersymmetric standard model and a dark supersymmetry model that allows for nonnegligible light boson lifetimes. In both cases, there is signiﬁcant improvement over previously published limits.


Introduction
The standard model (SM) is known to give an incomplete description of particle physics and a number of extensions of the SM predict the existence of new light bosons [1][2][3]. In this Letter, we present a model independent search for the pair production of a light boson that decays into a pair of muons. A simple example of pair production in proton-proton (pp) collisions is pp → h → 2a + X → 4μ + X, where h is a Higgs boson (either SM or non-SM), a is the new light neutral boson, and X are spectator particles that are predicted in several models [4]. While production via the h boson is possible, it is not required in the search presented here: the only requirement is that a pair of identical light bosons are created at a common vertex and each light boson subsequently decays to a pair of muons. These muon pairs are referred to as "dimuons"; the dimuon and new light boson production vertices are allowed to be displaced. The generic nature of this signature means that any limit set on the product of the cross section, branching fraction to dimuons squared, and acceptance is E-mail address: cms -publication -committee -chair @cern .ch. model independent; it can thus be reinterpreted in the context of specific models.
We develop a set of search criteria intended to minimize background events while remaining model independent. Two different classes of benchmark models are used to design the analysis and to verify that the results are actually model independent: the next-tominimal supersymmetric standard model (NMSSM) [1,[5][6][7][8][9][10][11][12] and supersymmetry (SUSY) models with hidden sectors (dark SUSY) [3,13,14]. In the NMSSM benchmark models, two of the three charge parity (CP) even neutral Higgs bosons h 1 or h 2 can decay to one of the two CP odd neutral Higgs bosons via h 1,2 → 2a 1 . The light boson a 1 subsequently decays to a pair of oppositely charged muons; this is equivalent to B(a 1 → 2μ). In the dark SUSY benchmark models, the breaking of a new U (1) D symmetry (where the subscript "D" means "Dark") gives rise to a massive dark photon γ D .
This dark photon can couple to SM particles via a small kinetic mixing parameter ε with SM photons. The lifetime, and thus the displacement, of the dark photon is dependent upon ε and the mass of the dark photon m γ D . The signal topologies investigated feature an SM-like Higgs boson h that decays via h → 2n 1 , where n 1 is the lightest non-dark neutralino. Both of the n 1 then decay via n 1 → n D + γ D , where n D is a dark neutralino that is undetected. The dark photon γ D decays to a pair of oppositely charged muons.

Data selection
The data were collected with a trigger that uses muon reconstruction algorithms that have an efficiency greater than 80% up to the maximum vertex displacement (98 mm) studied in this analysis [34]. This maximum vertex displacement is motivated in Section 4. The HLT is seeded by requiring the presence of two muons selected by the L1 trigger in an event, the leading muon with p T > 12 GeV, the subleading muon with p T > 5 GeV, and both satisfying |η| < 2.4. Events that later pass the HLT are required to have at least three reconstructed muons: one with p T > 15 GeV and |η| < 2.4, the other two with p T > 5 GeV and |η| < 2.4. The final state particles in the events are reconstructed using the particleflow (PF) algorithm which performs a global fit that combines information from each subdetector [35].
The offline event selection in this analysis requires events to have a primary vertex reconstructed using a Kalman filtering (KF) technique [36]. In addition, each event contains at least four muons, reconstructed with the PF algorithm, and identified as muons either by the PF algorithm itself or by using additional information from the calorimeter and muon systems. Each muon is required to have p T > 8 GeV and |η| < 2.4. At least one muon must be a "high-p T " muon, i.e., it must be found in the barrel region (|η| < 0.9) and must have p T > 17 GeV in order to ensure that the trigger reconstruction has high efficiency and has no dependence on η.
Dimuons are constructed from pairs of oppositely charged muons that share a common vertex, reconstructed using a KF technique, and must have an invariant mass m (μμ) less than 9 GeV. This restriction ensures that there is no contribution to the SM background from the Z boson decays nor the Y meson system. These muons pairs must not have any muons in common with one another. Exactly two dimuons must be present in each event. A dimuon that contains a high-p T muon is called a "high-p T dimuon". When only one high-p T muon is present in the event, the high-p T dimuon is denoted as (μμ) 1 , while the other is denoted as (μμ) 2 . When both dimuons have at least one high-p T muon, the dimuons are labeled randomly to prevent a bias in kinematic distributions. Single muons not included in dimuons are called "orphan" muons. No requirement is applied on the number of orphan muons. Each reconstructed dimuon must contain at least one muon that has at least one hit that is recorded by a layer of the pixel system. This requirement preserves the high reconstruction efficiency for our signal benchmark models. The dimuons are required to originate from the same primary ver- is the z position of the secondary vertex associated with the dimuon propagated back to the beamline along the dimuon direction vector. Furthermore, each dimuon must be sufficiently isolated. The dimuon isolation I (μμ) is calculated as the p T sum of charged-particle tracks with p T > 0.5 GeV in the vicinity of the dimuon within R < 0.4 and |z track − z (μμ) | < 0.1 cm. Here, R is defined in terms of the track separation in η and azimuthal angle (φ, in radians) as 2 , while z track is defined as the z coordinate of the point of closest approach to the primary vertex along the beam axis. Tracks included in the dimuon reconstruction are excluded from the isolation calculation. The total isolation sum must be less than 2 GeV. Since the dimuons are expected to originate from the same type of light bosons, the dimuon masses should be consistent with each other to within five times the detector resolution. This requirement carves out a signal region (SR) in the two-dimensional plane of the dimuon invariant masses m (μμ) 1 and m (μμ) 2 . The signal region is illustrated in Fig. 1 (left).

Signal modeling
The pp collisions at √ s = 13 TeV are simulated for samples in each of the two benchmark models, NMSSM and dark SUSY. The parton distribution functions (PDFs) are modeled using NNPDF2.3LO [37]. The underlying event activity at the LHC and jet fragmentation is modeled with the Monte Carlo (MC) event generator pythia [38] using the "CUETP8M1" tune [39]. Specifically, pythia 8.212 is used for NMSSM and pythia 8.205 for the dark SUSY models. In each model, only Higgs boson production via gluon-gluon (gg) fusion is considered. A single mass point is also generated through vector boson fusion (VBF) and associated vector boson production (VH) to determine their contribution to the h 2 → 2a 1 rate; this is included in a simplified reference scenario discussed later.
In the case of the NMSSM, a simulated Higgs boson, either h 1 or h 2 (generically denoted by h 1,2 ), is forced to decay to a pair of light bosons a 1 . Each a 1 subsequently decays to a pair of oppositely charged muons. Since the h 1,2 in h 1,2 → 2a 1 might not be the observed SM Higgs boson [40][41][42], mass values of m h 1,2 between 90 and 150 GeV are simulated. This range is motivated by constraints set by the relic density measurements from WMAP [43] and Planck [44], as well as searches at LEP [45][46][47][48][49][50]. The light boson mass is simulated to vary between 0.25 and 3.55 GeV, or approximately 2m μ and 2m τ , as motivated in Ref. [51].
In the case of dark SUSY, production of SM Higgs bosons is simulated with the MC matrix-element generator MadGraph 4.5.2 [52] at leading order. The non-SM decay of the Higgs bosons is modeled using the BRIDGE 2.24 program [53]. Higgs bosons are forced to decay to a pair of SUSY neutralinos n 1 via h → 2n 1 . Each SUSY neutralino in turn decays to a dark photon and a dark neutralino via n 1 → n D + γ D . The dark neutralino mass m n D is set to 1 GeV; they are considered stable and thus escape detection. We set the dark photons to decay to a pair of oppositely charged muons 100% of the time, γ D → μ − μ + . Only signal events are generated because these MC generated events are used to determine the effect of the selection criteria on the signal. The Higgs boson and n 1 masses are fixed to 125 and 10 GeV, respectively. Dark photon masses m γ D are simulated between 0.25 and 8.5 GeV. The upper value was chosen such that any observed peak will be fully below the 9 GeV limit described in Section 3. Since dark photons interact weakly with SM particles, their decay width is negligible compared to the resolution in the dimuon mass spectrum. Muon displacement is modeled with an exponential distribution with cτ γ D between 0 and 100 mm. All MC generated events are run through the full CMS simulation based on Geant4 [54] and reconstructed with the same algorithms that are used for data.
One of the key features of this analysis is the model independence of the results. This is confirmed by verifying that the ratio of the full reconstruction efficiency full over the generator level acceptance α gen is independent of the signal model. The signal acceptance is defined as the fraction of MC-generated events that pass the generator level selection criteria. The criteria are as follows: at least four muons in each event with p T > 8 GeV and |η| < 2.4, at least one muon with p T > 17 GeV and |η| < 0.9, and both light bosons must have a transverse decay length L xy < 9.8 cm and longitudinal decay length |L z | < 46.5 cm. The upper limits on L xy and |L z | correspond to the dimensions of the outer layer of the CMS pixel system and define the volume in which a new light boson decay can be observed in this analysis. The parameter full is defined as the fraction of MC-generated events that pass the trigger and full offline selection described above. The insensitivity to the model used is displayed in Table 1.
Scale factors are determined to correct for the differences between observed data and simulated samples. Corrections for the identification and isolation of muons and isolation of dimuons are measured using Z → μ − μ + and J/ψ → μ − μ + samples using a "tag-and-probe" technique [55]; the samples used are events from simulated data and from observed data control regions enriched in events from the aforementioned SM processes. All muons in these samples are required to have p T > 8 GeV, the "tag" muon is required to be a loose muon as described in Ref. [30], while the "probe" muon criteria vary according to the variable under study.
Corrections for the trigger efficiency are calculated using WZ → 3μ and ttZ → 3μ events in simulated samples and in control data samples enriched with those processes. The control data samples are selected using a missing transverse energy requirement such that the control data sample is primarily composed of events that are different from those in the data sample used in this analysis.
A scale factor per event obtained from the efficiency seen in data, data , compared to the efficiency seen in MC generated data, sim , is determined to be data / sim = 0.93 ± 0.06 (stat).

Background estimation
The selection criteria described in Section 3 are effective at reducing and eliminating most SM backgrounds with similar topology to our signal. As a result, this analysis is expected to have a very small background contribution in the SR. Three SM backgrounds are found to be nonnegligible and are presented here: bottom quark pair production (bb), prompt double J/ψ meson decays, and electroweak production of four muons. Contributions from Y mesons are also considered; they are found to be negligible below the 8.5 GeV upper bound on the mass of the new light boson. Cosmic ray backgrounds are negligible. The total background contribution in the SR is estimated to be 7.95 ± 1.12 (stat) ± 1.45 (syst) events; the contributions from each process are described below.

The bb background
The largest background, bb production, is dominated by events in which both b quarks decay to μ − μ + + X or decay through lowmass meson resonances such as ω, ρ, φ, J/ψ, and ψ(2S). The J/ψ meson decay contribution considered in this background is nonprompt; the prompt J/ψ meson decay contribution is discussed in Section 5.2. A minor contribution comes from events with charged particle tracks misidentified as muons. A two-dimensional tem- is constructed in the plane of the two dimuon invariant masses and used to estimate the contribution to the SM background from bb decays. The template is constructed as follows.
First, a bb-enriched control sample is selected from events with similar kinematic properties as the signal events, but not included in the SR. Events are required to pass the signal trigger and have exactly three muons. One of these muons must have p T > 17 GeV within |η| < 0.9, while the other two have p T > 8 GeV within |η| < 2.4. In addition, the control sample selection requires a good primary vertex, exactly one dimuon, and one orphan muon. The longitudinal distance between the projections of the dimuon trajectory starting from its vertex and the orphan muon track back to the beam axis, z((μμ), μ orphan ) must have an absolute value of less than 0.1 cm. The dimuon is required to have at least one hit Table 1 The full reconstruction efficiency over signal acceptance full /α gen in % for several representative signal NMSSM (upper) and dark SUSY benchmark models (lower). All uncertainties are statistical. in the pixel system as explained in Section 3. Finally, the dimuon isolation value cannot be higher than 2 GeV.
Next, two one-dimensional templates, S I (m (μμ) ) and S II (m (μμ) ), are obtained from the bb-enriched events. In the case of S I (m (μμ) ), at least one high-p T muon is contained in the dimuon. In the case of S II (m (μμ) ), the high-p T muon is the orphan muon and the dimuon may or may not contain another high-p T muon. This procedure ensures that kinematic differences between signal events that have exactly two high-p T dimuons or just one high-p T dimuon are taken into account. Each distribution is fitted with a shape comprised of a Gaussian distribution for each light meson resonance, a double-sided Crystal Ball function [56] for the J/ψ meson signal peak, and a set of sixth-degree Bernstein polynomials for the bulk background shape. The template S(m (μμ) 1 Finally, the two-dimensional template is normalized in the dimuon-dimuon mass space from 0.25 to 8.5 GeV. The template is represented as a function of m (μμ) 1 and m (μμ) 2 in Fig. 1 (left) by a gray scale. The SR defined in Section 3 is outlined by dashed lines. The region of the mass space outside the SR represent the control region for the bb background. The ratio between the integral of the template in the SR A SR and the control region A CR is calculated to be R = A SR /A CR = 0.1444/0.8556. The same figure also shows the 43 events found in the data that pass all selection criteria except for the m (μμ) 1 m (μμ) 2 requirement and thus fall outside the SR. The number of bb events in the SR is then estimated to be This method of estimating the bb contribution to background events is further validated by repeating the procedure for different dimuon isolation values (5, 10, 50 GeV) and without any isolation. The bb event yield is stable in the SR within 20%, which is assigned as a systematic uncertainty.

Prompt double J/ψ meson background
Two mechanisms contribute to prompt double J/ψ meson production: single parton scattering (SPS) and double parton scattering (DPS); these processes have been measured by CMS and ATLAS [57,58]. They can mimic the signal process when each J/ψ meson decays to a pair of muons. The prompt double J/ψ meson decay background is estimated with a method that uses both experimental and simulated data. In a control sample of experimental data, the prompt and nonprompt double J/ψ meson decay contributions are separated using the matrix method (also called the "ABCD" method [59]). The prompt contribution is then extrapolated into the SR. Double J/ψ meson events are selected with a trigger dedicated to bottom quark physics. Each event is required to have at least four muons with p T > 3.5 GeV within |η| < 2.4. No high-p T muon is required. Events must have exactly two dimuons, with labels (μμ) 1 or (μμ) 2 assigned randomly. The dimuon isolation follows the same definition as in Section 3. The kinematic properties of SPS and DPS events are studied using MC simulation. These events are generated using pythia 8.212 and herwig 2.7.1 [60]. The variable with the best SPS-DPS separation power is found to be the absolute difference in rapidity between the two dimuons, | y|. To remove nonresonant muon pairs from the sample, the dimuon masses are required to be within 2.8 and 3.3 GeV. The ABCD method is then employed using the dimuon isolation values as uncorrelated variables in the plane (I (μμ) 1 , I (μμ) 2 ). The maximum isolation on (μμ) 1 and (μμ) 2 is set to 12 GeV. Here, region "A" is the region bounded by I (μμ) 1,2 < 2 GeV. Conversely, "B", "C", and "D" are nonisolated sideband regions used to extrapolate the nonprompt contribution into region "A". The nonprompt | y| distribution is determined from the sideband regions; this distribution is scaled to match the nonprompt contribution in region "A". This is then subtracted from the | y| distribution, leaving the prompt | y| distribution in region "A". To separate the prompt SPS from prompt DPS in data, a template distribution f SPS | y SPS | + (1 − f SPS )| y DPS | is fitted to the corresponding | y| distribution in data, where f SPS and 1 − f SPS are the fractions of prompt SPS and DPS events, respectively. Finally, this result is used to determine the number of events that are expected in the SR of our experimental data sample. The contribution of the prompt double J/ψ meson decay events in data passing the signal selections in Section 3 is calculated to be N data (SR) = 0.33 ± 0.08 (stat) ± 0.05 (syst).

Electroweak background
Electroweak production of four muons, pp → 4μ, is estimated using MC events generated with CalcHEP 3.6.25 [61]. The processes studied include qq → ZZ * → 2μ − 2μ + and qq → Z → μ − μ + , where one of the muons radiates a second Z boson that decays to a μ − μ + pair. Other electroweak processes, such as pp → h(125) → ZZ * → 2μ − 2μ + , are determined to be negligible a priori and thus are not included. Based on the simulation, the electroweak background is found to be 0.36 ± 0.09 (stat). Unlike the prompt double J/ψ meson decay background, the electroweak background is not concentrated at any particular mass value; its contribution to any mass bin is negligible compared to the bb background. Consequently, these background events are neglected in any limit setting computation.

Systematic uncertainties
Both instrumental and theoretical sources of uncertainty are considered in this section. The leading source of instrumental uncertainty is the triple-muon trigger scale factor (6%). It is dominated by the statistical uncertainty in events in the control region used to measure the scale factor. Other sources of instrumental uncertainty include the uncertainty in the measurement of the integrated luminosity recorded by the CMS detector (2.5%) [62], the muon identification data-to-simulation scale factor (0.6% per muon for all simulated muons), the reconstruction of the dimuon in the tracker (1.2% per dimuon) and in the muon system (1.3% per dimuon) from spatially close muons, and the effect on the acceptance of the dimuon mass shape used to determine the width of the SR (1.5%). The uncertainty in the dimuon isolation and the contributions of extraneous pp collisions are determined to be negligible.
The theoretical uncertainties are dominated by the uncertainty in the PDFs, knowledge of the strong coupling constant α S , and the renormalization (μ R ) and factorization (μ F ) scales. The PDF and α S uncertainties are estimated using a technique that follows the PDF4LHC recommendations [63,64]. The uncertainty in the scale factors is determined by simultaneously varying μ R and μ F up and down by a factor of two using MCFM 8.0 [65]. The effect of PDF choice and PDF parameter variation upon the central values is also studied. When all previously described theoretical uncertainties are added in quadrature, the sum is 8%. The uncertainty in the branching fraction B(h → 2a + X → 4μ + X) is taken to be 2% [42].

Results
After applying all selection criteria to the data sample, 9 events are found in the SR. Their distribution in m (μμ) 1 and m (μμ) 2 is shown in Fig. 1 (left). This result is consistent with the sum of all background estimates described in Section 5, which is found to be 7.95 ± 1.12 (stat) ± 1.45 (syst) events. A model independent 95% confidence level (CL) upper limit is set on the product of the production cross section times branching fraction to dimuons squared times acceptance. Limits are set using the CL s method [66,67]. The test statistic used is based on the logarithm of the likelihood ratio [68]. The systematic uncertainties and their correlations have been accounted for by profiling the likelihood with respect to the nuisance parameters for each value of the signal strength s; this results in the profile likelihood being a function only of s. The limit is shown as a function of m a in Fig. 1 (right) over the range 0.25 < m a < 8.5 GeV; the limit varies between 0.15 and 0.39 fb. Neglecting the large peak in the upper limit at the J/ψ meson mass, the largest upper limit is 0.25 fb. This result can be interpreted in the context of specific models.
For the NMSSM scenario, the 95% CL upper limit is derived for to simplify the expression. This choice is conservative because 1 , for any m h 1 . In this simplified scenario, B(a 1 → 2μ) is a function of m h 1 as calculated in Ref. [51].
To facilitate comparison between the upper limits derived from this analysis and upper limits following from setting parameters in theoretical models, we include reference curves (solid line) in both Fig. 2 left and right. For both reference curves, the ratio of the vacuum expectation values of the Higgs doublets tan β is set to 20. We also set σ (pp → h i ) = σ SM (m h i ) [69] and B(m h i → 2a 1 ) = 0.3% so that the resulting reference curves are similar to the upper limits that are determined from the yield of dimuon pair events observed in the data. In Fig. 2 (left), the reference curve is constructed with the assumption that B(a 1 → 2μ) = 7.7% and m a 1 ≈ 2 GeV. In the region where m h i < 125 GeV, m h 1 is the independent variable and it is assumed that m h 2 is the mass of the observed 125 GeV Higgs boson. In the region where m h i > 125 GeV, m h 2 is the independent variable and it is assumed that m h 1 is the observed Higgs boson mass. Compared to the upper limits shown in Refs. [15], Fig. 2 (left) represents an improvement of a factor of ≈1.5 for m a 1 = 3.55 GeV (dotted curve) and a factor of ≈3 for m a 1 = 0.25 GeV (dashed curve). In Fig. 2 (right), we present 95% CL upper limits as functions of m a 1 in the NMSSM scenario on m h 1 = 125 GeV (dash-dotted curve), and m h 2 = 150 GeV (dotted curve). It is assumed that all contributions come from either h 1 or h 2 ; there is no case in which both h 1 and h 2 decay to the a 1 . The sharp inflections in the reference curve are due to the fact that B(a 1 → 2μ) is affected by the a 1 → ss and a 1 → gg channels [51]. As m h 1 crosses the internal quark loop thresholds, B(a 1 → gg) changes rapidly, giving rise to structures in B(a 1 → 2μ) at these values of m h 1 .
For the dark SUSY scenario, a 90% CL upper limit is set on the product of the Higgs boson production cross section and the branching fractions of the Higgs boson (cascade) decay to a pair of dark photons. The limit set by this experimental search is presented in Fig. 3 as areas excluded in a two-dimensional plane of ε and m γ D . Also included in Fig. 3    namely, τ γ D (ε, m γ D ) = ε −2 f (m γ D ). The lifetime of the dark photon is allowed to vary from 0 to 100 mm and m γ D can range from 0.25 to 8.5 GeV. Because of the extensions in the ranges of these parameters, this search constrains a large and previously unexplored area in the ε and m γ D parameter space. The limits on ε presented in this Letter improve on those in Ref. [15] by a factor of approximately 2.5.

Summary
A search for pairs of new light bosons that subsequently decay to pairs of oppositely charged muons is presented. This search is developed in the context of a Higgs boson decay, h → 2a + X → 4μ + X and is performed on a data sample collected by the Compact Muon Solenoid experiment in 2016 that corresponds to an integrated luminosity of 35.9 fb −1 proton-proton collisions at 13 TeV. This data set is larger and collected at a higher centerof-mass energy than the previous CMS search [15]. Additionally, both the mass range of the light boson a and the maximum possible displacement of its decay vertex are extended compared to the previous version of this analysis. Nine events are observed in the signal region (SR), with 7.95 ± 1.12 (stat) ± 1.45 (syst) events expected from the standard model (SM) backgrounds. The distribution of events in the SR is consistent with SM expectations. A model independent 95% confidence level upper limit on the product of the production cross section times branching fraction to dimuons squared times acceptance is set over the mass range 0.25 < m a < 8.5 GeV and is found to vary between 0.15 and 0.39 fb.
This model independent limit is then interpreted in the context of dark supersymmetry (dark SUSY) with nonnegligible light boson lifetimes of up to cτ γ D = 100 mm and in the context of the next-tominimal supersymmetric standard model (NMSSM). For the dark SUSY interpretation, the upper bound of m γ D was increased from 2 to 8.5 GeV and the excluded ε was improved by a factor of approximately 2.5. In the NMSSM, the 95% CL upper limit was improved by a factor of ≈1.5 (3) for m a 1 = 3.55 (0.25) GeV over previously published limits.

Acknowledgements
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.