Search for an $L_{\mu}$ $-$ $L_{\tau}$ gauge boson using Z $\to$ 4$\mu$ events in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for a narrow Z$'$ gauge boson with a mass between 5 and 70 GeV resulting from an $L_{\mu}$ $-$ $L_{\tau}$ $U(1)$ local gauge symmetry is reported. Theories that predict such a particle have been proposed as an explanation of various experimental discrepancies, including the lack of a dark matter signal in direct-detection experiments, tension in the measurement of the anomalous magnetic moment of the muon, and reports of possible lepton flavor universality violation in B meson decays. A data sample of proton-proton collisions at a center-of-mass energy of 13 TeV is used, corresponding to an integrated luminosity of 77.3 fb$^{-1}$ recorded in 2016 and 2017 by the CMS detector at the LHC. Events containing four muons with an invariant mass near the standard model Z boson mass are analyzed, and the selection is further optimized to be sensitive to the events that may contain Z $\to$ Z$'\mu\mu$ $\to$ 4$\mu$ decays. The event yields are consistent with the standard model predictions. Upper limits of 10$^{-8}$-10$^{-7}$ at 95% confidence level are set on the product of branching fractions $\mathcal{B}$(Z $\to$ Z$'\mu\mu$) $\mathcal{B}$(Z$'$ $\to$ $\mu\mu$), depending on the Z$'$ mass, which excludes a Z$'$ boson coupling strength to muons above 0.004-0.3. These are the first dedicated limits on $L_{\mu}$ $-$ $L_{\tau}$ models at the LHC and result in a significant increase in the excluded model parameter space. The results of this search may also be used to constrain the coupling strength of any light Z$'$ gauge boson to muons.


Introduction
The standard model (SM) of particle physics [1][2][3] can not explain all experimental observations to date and is, therefore, generally believed to be an incomplete theory. Enlarging the SM gauge group to include an additional U(1) symmetry is a simple and well-motivated extension [4], which leads to a prediction of a new vector particle, a Z boson. In order for the extended gauge symmetry to be anomaly free, only certain generation-dependent couplings are allowed. The anomaly-free model we consider in this paper is the L µ − L τ gauge symmetry [5], where L µ and L τ are the µ and τ lepton numbers, respectively. The interaction between the Z and the second-and third-generation leptons can be described with the following Lagrangian [6]: where g is an arbitrary dimensionless coupling to the SM left-handed and right-handed µ and τ multiplets, respectively, are Additional U(1) gauge symmetries based on the difference in lepton family numbers are all anomaly free and require no new fermionic particle content. The model based on gauging L µ − L τ in particular is the least constrained experimentally, since it is coupled only to secondand third-generation leptons.This model has gained popularity in recent years [7][8][9][10][11][12][13] as an explanation to several anomalous experimental measurements in particle physics. These anomalies include the measurement of the anomalous muon magnetic moment by the Muon g−2 Collaboration [14], which can be explained for certain values of the Z mass and coupling strength (g) [7,9]. In addition, if the Z mediates an interaction between dark matter and ordinary matter, the bounds on the dark matter coupling strength from direct-detection experiments are less stringent [10,11] since the Z does not couple directly to quarks. Finally, if additional interactions beyond the minimal L µ − L τ model are assumed, abnormalities in kinematic angular distributions and lepton flavor universality tests observed in B → K * µ + µ − decays [15,16] can be explained by this model, given its flavor non-universal couplings [8,11].
The Z gauge boson associated with the putative L µ − L τ gauge symmetry can be sought at the CERN LHC. Since the Z couples only to second-and third-generation leptons (µ, ν µ , τ and ν τ ), it must be produced as a final state radiation product of a lepton originating from some other physics process. The Z → 4µ decay provides an extremely clean source of muons with excellent mass resolution. The resonant signal decay Z → µµ that may be present in Z → 4µ decays further reduces the background contamination. There are two types of irreducible background where the additional dimuon pair originates from annihilation or conversion topologies. The Feynman diagrams in Fig. 1 (left) for the signal and in Fig. 1 (right) for the background are examples of the annihilation topology, while the diagrams in Fig. 2 are examples of the conversion topology. The dominant background to the search comes from the annihilation diagram in Fig. 1 (right), while the background originating from the conversion diagrams in Fig. 2 are subdominant. The signal and background processes originating from the diagrams in Fig. 1 will hereafter be collectively referred to as the Z → 4µ process since they have very similar kinematic properties. The background processes originating from the diagrams in Fig. 1 (right) and Fig. 2 (left) will be collectively referred to as qq → 4µ, and the process originating from the diagram in Fig. 2 (right) will be referred to as gg → 4µ.
The Z → 4µ process has been studied by the ATLAS and CMS Collaborations [17-19] and constraints on the L µ − L τ parameter space have been derived. However, these measurements are not optimized for the presence of a Z particle. In particular, they do not utilize the fact that two of the four muons would form a resonant peak at the Z mass, providing a means to reduce the dominant background by several orders of magnitude. The subject of this paper is a dedicated counting experiment search with a final selection based on the reconstructed Z candidate mass. Figure 1: Leading order Feynman diagrams for the signal process (left) and the dominant background process (right), where in each diagram the additional dimuon pair originates from annihilation. Figure 2: Leading order Feynman diagrams for the subdominant quark-initiated (left) and gluon-initiated (right) background processes, where in each diagram the dimuon pairs originate from conversion.

The CMS detector
A 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. [20]. 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 detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. The silicon tracker measures charged particles with |η| < 2.5. For nonisolated particles with transverse momentum (p T ) between 1 and 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90 (45-150) µm in the transverse (longitudinal) impact parameter [21]. Muons are measured in the region |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive plate chambers. Matching muons to tracks measured in the silicon tracker results in a relative p T resolution for muons with 20 < p T < 100 GeV of 1.3-2.0% in the barrel (|η| < 0.9) and better than 6% in the endcaps (|η| > 0.9). The first level of the CMS trigger system [22], composed of custom hardware processors, uses information from the calorimeters and muon detectors to select the most interesting events in a fixed time interval of less than 4 µs. The high-level trigger (HLT) processor farm further decreases the event rate from around 100 kHz to less than 1 kHz before data storage.

Data and simulated samples
This analysis makes use of proton-proton (pp) collision data recorded by the CMS detector in 2016 and 2017, corresponding to an integrated luminosity of 77.3 fb −1 . Collision events are selected by HLT algorithms that require the presence of one, two, or three muons passing loose identification and isolation requirements. The overall trigger efficiency for simulated signal events that pass the full selection chain of this analysis (described in Section 4) is greater than 99%. The trigger efficiency is measured in data with a method based on the "tag-and-probe" technique [23] using a sample of 4µ events collected by the single-muon triggers. Events with four muons contain predominantly prompt muons and, therefore, background subtraction is not necessary. Muons matched to the single-muon triggers are used as tags and the other three muons are used as probes. The probe muons are then matched to the triggering muon objects from any of the one, two, or three muon triggers, and the combined efficiency is extracted. The efficiency in data is found to be in agreement with the expectation from simulation.
Monte Carlo simulation samples for the Z signal and for the background coming from the qq → 4µ and gg → 4µ processes are used to optimize the event selection, evaluate the signal acceptance, and estimate the background rate and systematic uncertainties. The signal is generated at leading order (LO) in perturbative quantum chromodynamics (pQCD) with MAD-GRAPH5 aMC@NLO (v2 4 2) [24] together with the Universal FEYNRULES Output (UFO) model from Ref. [6]. The signal samples are generated with m(Z ) ranging from 5 to 70 GeV in steps of 5 GeV. For m(Z ) below 5 GeV nonprompt muons become a challenging background, and for m(Z ) above 70 GeV the Z boson starts to be produced off mass-shell, requiring a dedicated event selection. The qq → 4µ process is generated at next-to-LO (NLO) in pQCD with POWHEG 2.0 [25][26][27], while the gg → 4µ process is generated at LO with MCFM 7.0 [28]. The default set of parton distribution functions (PDFs) used in all simulations is NNPDF30 nlo as 0118 [29]. The fully differential cross section for the qq → 4µ process has been computed at next-to-NLO (NNLO) [30], and the appropriate NNLO/NLO correction factor K of 1.03 at m(4µ) = m(Z) is used to correct the POWHEG sample. Because the signal production process is the same as that of the qq → 4µ background, an analogous NNLO/LO K factor of 1.29 is used to correct the signal process. The gg → 4µ process contributes at NNLO in pQCD and is corrected by a K factor of 2.4 [31][32][33][34][35][36][37].
After the final selection, described in Section 5, the gg → 4µ background contribution is typically less than 1% (and at most 7%) of the the qq → 4µ contribution. These simulations have been found to provide an accurate description of 4µ events in data by several previous

Object reconstruction
The techniques of the object reconstruction and event selection are based largely on Refs.
[19, [38][39][40]. Event reconstruction is based on the particle-flow (PF) algorithm [45], which exploits information from all the CMS subdetectors to identify and reconstruct individual particles in the event. Higher-level observables, such as muon isolation quantities, are built from the PF candidates classified as charged hadrons, neutral hadrons, photons, electrons, or muons.
Muons are reconstructed within the geometrical acceptance |η| < 2.4 by combining information from the silicon tracker and the muon system [46], and are required to satisfy p T > 5 GeV. The inner and outer tracks are matched using either an outside-in algorithm, starting from a track in the muon system, or an inside-out algorithm, starting from a track in the silicon tracker.
In the latter case, some very low-p T muons that may not have sufficient energy to penetrate the entire muon system are also collected by considering tracks that match track segments in only one or two planes of the muon system. Muons are identified from the reconstructed muon track candidates by applying minimal requirements on the inner and outer tracks, taking into account their compatibility with small energy deposits in the calorimeters.
Muons originating from nonprompt decays of hadrons are suppressed by requiring each muon track to have the ratio between its impact parameter in three dimensions, computed with respect to the chosen primary vertex position, and its uncertainty to be less than 4. The primary pp interaction vertex is taken to be the reconstructed vertex with the largest value of summed physics-object p 2 T . The physics objects are the jets, clustered using the jet finding algorithm [47,48] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
A relative isolation requirement of I µ < 0.35 is imposed to discriminate between prompt muons from Z boson decays and those arising from electroweak decays of hadrons within jets, where the relative isolation is defined as The isolation sums involved are all restricted to PF candidates within a volume bounded by a cone of angular radius ∆R = 0.3 around the muon direction at the primary vertex, where the angular distance between two particles i and j is ∆R where φ is the azimuthal angle in radians. The quantity ∑ p An algorithm is used to recover the final-state radiation (FSR) photons from muons. Photons reconstructed by the PF algorithm within |η γ | < 2.4 are required to satisfy p γ T > 2 GeV and I γ < 1.8. The photon relative isolation I γ is defined as for the muons in Eq. (3). Every FSR candidate is associated with the closest selected muon in the event, and we require FSR candidates to satisfy ∆R(γ, µ)/(p γ T ) 2 < 0.012 GeV −2 and ∆R(γ, µ) < 0.5. Finally, for every muon we retain the FSR candidate, if any, with the lowest ∆R(γ, µ)/(p γ T ) 2 . About 5% of signal events are found to have one FSR photon attached. Any selected FSR photons are excluded from the corresponding muon isolation computation.
The decay products of known dimuon resonances (J/ψ meson, Z boson) are used to calibrate the muon momentum scale and resolution in bins of p T and η. Muon momenta are calibrated using a Kalman filter approach [49]. A tag-and-probe technique [23, 50] is used to measure the efficiency of the reconstruction and selection for prompt muons in several bins of p T and η. The difference between the efficiencies measured in simulation and data, which on average is 1.2% per muon, is used to correct the selection efficiency in the simulated samples. The combined muon reconstruction and identification efficiency for signal events, including these corrections, is about 92% per muon.

Event selection
Events are required to contain at least four well-identified and isolated muons, with at least two muons required to have p T > 10 GeV and at least one to have p T > 20 GeV. The four selected muons must have zero net charge. Dimuon candidates are formed from muon pairs of opposite sign (µ + µ − ) and are required to pass 4 < m µ + µ − < 120 GeV. All recovered FSR photon candidates are included in the invariant mass computation. The dimuon candidates are then combined into Z → 4µ candidates, and we define Z 1 to be the dimuon candidate with an invariant mass closest to the nominal Z boson mass (m Z ) [51], and Z 2 as the other one. In events with more than four muons, the Z → 4µ candidate with m(Z 1 ) closest to m Z is selected. If several Z → 4µ candidates have the same m(Z 1 ), the Z 2 candidate formed from the two muons with the highest scalar sum of p T is chosen.
To be considered for the analysis, Z → 4µ candidates have to pass a set of kinematic requirements. The Z 1 invariant mass must be larger than 12 GeV and all muons must be separated in angular space by at least ∆R(µ i , µ j ) > 0.02. To further suppress events with muons originating from hadron decays in jet fragmentation or from the decay of low-mass hadronic resonances, all four opposite sign muon pairs that can be constructed with the four muons are required to satisfy m µ + µ − > 4 GeV, where selected FSR photons are disregarded in the invariant mass computation. Finally, the four-muon invariant mass m(4µ) must be between 80 and 100 GeV. The Z candidate is most often reconstructed as Z 2 for m(Z ) < 42.65 GeV and as Z 1 for m(Z ) > 42.65 GeV. The search is a counting experiment with a sliding mass window, and a final selection made on either m(Z 1 ) or m(Z 2 ) values, depending on the Z mass hypothesis. The exclusion limit for m(Z ) = 42.65 GeV using either m(Z 2 ) or m(Z 1 ) as an observable is about the same. For m(Z ) < 42.65 GeV, m(Z 2 ) is required to be within 2% of the m(Z ). While for m(Z ) > 42.65 GeV, the same requirement is applied on m(Z 1 ). The search window size of 2% was chosen to simultaneously optimize the expected significance and exclusion limit for different Z mass hypotheses. The efficiency of this requirement varies with Z mass and is 60% for m(Z ) = 5 GeV, 25% for m(Z ) = 40 GeV, and 75% for m(Z ) = 70 GeV. The low efficiency for m(Z ) ≈ m Z /2 is due the combinatoric ambiguity in selecting the correct Z candidate from the four possible dimuon candidates. The selection is, however, still extremely beneficial, as it eliminates approximately 99.8% of the SM Zγ * background for m(Z ) = 40 GeV. Additional backgrounds to the signal that can arise from processes in which heavy-flavor jets produce secondary muons, and from processes in which decays of heavy-flavor hadrons or nonprompt decays of light mesons within jets are misidentified as prompt muons, are found to be negligible after the final event selection.

Systematic uncertainties
Experimental uncertainties that equally affect the signal and background estimations include the uncertainty in the integrated luminosity measurement of 2.5% [52] and 2.3% [53] for the 2016 and 2017 data sets, respectively, and the uncertainty in the muon identification and reconstruction efficiency (4.9% on the overall event yield). An uncertainty in the signal and background yields due to the muon momentum scale is determined using Z → 4µ events in data and simulation and found to be negligible (0.1%). An uncertainty in the signal and background yields of 2% coming from the determination of the muon momentum resolution is obtained by smearing the dimuon mass resolution by 20% [40] with respect to the nominal resolution and recomputing the expected yields. The uncertainties due to the finite sizes of the simulated samples amount to 3% for the background estimation and 1.4% for the signal estimation.
Theoretical uncertainties that affect both the signal and background estimations include uncertainties in the finite-order perturbative calculations and the choice of the PDF set. The uncertainties arising from finite-order perturbative calculations are estimated by varying the renormalization and factorization scales between 0.5 and 2 times their nominal value, while keeping their ratios between 0.5 and 2. This uncertainty is found to be 3.5 (3.9)% for the qq → 4µ (gg → 4µ) process and is taken to be correlated between the signal and the dominant qq → 4µ background process. Following Ref. [31] and taking into account differences in selection, an additional uncertainty of 10% in the K factor used for the gg → 4µ prediction described in Section 3 is applied to account for the fact that the K factor was computed for the closely related gg → H process. The uncertainty from the PDF set is determined following the PDF4LHC recommendations [54] and is found to be 3.1 (3.5)% for the qq → 4µ (gg → 4µ) process. This uncertainty is also taken to be correlated between the signal and the dominant qq → 4µ background process.
To estimate the effect of the interference between the signal and background processes, three types of samples, pp → 4µ (inclusive), pp → Z µµ → 4µ (signal only), and pp → 4µ (background only), are generated using MADGRAPH5 aMC@NLO (v2 4 2), with g varied from 0.01 to 0.5. The inclusive sample contains background, signal, and interference contributions. The effect of the interference on the normalization of the signal is estimated by taking the difference between the inclusive sample's cross section and the sum of the signal and background sample's cross sections. This difference is at most 5% after the final event selection, and additional 5% uncertainty in the signal yield is applied to account for this effect.

Results
The number of candidates observed in data and the expected yields for the backgrounds and the different Z signals after the full event selection are reported in Table 1. The reconstructed four-muon invariant mass distributions are shown in Fig. 3 and compared with the expectations from signal and background processes. Fig. 4 shows the reconstructed m(Z 1 ) and m(Z 2 ) distributions. Table 1: The numbers of expected background and signal events and the numbers of observed candidate events after the full selection with 80 < m 4µ < 100 GeV. The signal and qq/gg → 4µ background rates are both estimated from simulation. The signal predictions are reported with systematic uncertainties only, while the background predictions are reported with statistical and systematic uncertainties, respectively. Also shown are the numbers of expected background and signal events and the numbers of observed candidate events in the relevant mass windows for three m(Z ) hypotheses. The values of the coupling strengths are chosen for the purpose of illustration.   In all cases, the observed distributions agree with the expectations within the assigned uncertainties. Upper limits at 95% confidence level (CL) are derived on the product of the Z µµ production cross section and the branching fraction B(Z → µµ) using the CL s method [55,56] with the test statistic described in Ref. [57], in the asymptotic approximation [58]. The asymptotic approximation was verified to be valid by computing limits with the full CL s method using pseudo-experiments for several m(Z ) hypotheses. A linear interpolation of the expected event yields between generated signal MC samples is assumed in the limit calculations. Systematic uncertainties are incorporated into the likelihood as nuisance parameters with lognormal prior constraints. Due to the low number of events passing the final selection, the statistical uncertainty is always larger than 22% within the entire m(Z ) search region, and dominates the sensitivity of this analysis. These limits are shown in Fig. 5. The upper limits on the B(Z → Z µµ)B(Z → µµ) are also shown. For the derivation of branching fraction limits, the Z boson production cross section prediction computed at NNLO in pQCD with the program FEWZ 2.1 [59][60][61] is used. Upper limits are also derived on the gauge coupling strength g and compared to other experimental constraints, shown in Fig. 6. These limits assume the B(Z → µµ) is equal to 1/3 as in the minimal L µ − L τ model with equal left-and right-handed coupling strengths, and the additional constraints are adapted from Ref. [11]. The mass of the dark matter candidate in the model from Ref.
[11] is assumed to be much larger than the largest Z mass considered and the gauge coupling strengths to other particles, such as b-and s-quarks, are taken to be much smaller than the coupling strength to leptons so that B(Z → µµ) is constant. The natural width of the Z is also assumed to be less than the detector resolution, which is a valid approximation in the minimal L µ − L τ model when g 2 /4π < 0.01. The shaded yellow region shows constraints derived in Ref.
[11] from the ATLAS B(Z → 4µ) measurement at √ s = 7 and 8 TeV [18]. The shaded red region is excluded by the measurement of the so-called neutrino trident cross section by the CCFR Collaboration [62,63]. The green region is excluded by a global analysis of B s mixing measurements performed in Ref. [11]. The region in between those two constraints and for m(Z ) > 10 GeV is a candidate region to explain the LHCb B decay anomalies. It is important to note that in order to explain these anomalies, additional assumptions on the couplings of the Z boson to b-and s-quarks are required, and the constraints from B s mixing measurements are therefore not generally applicable to the minimal L µ − L τ model. It can be seen that this search is able to exclude a significant portion of the previously allowed parameter space.

Summary
The first dedicated search at the LHC for a Z gauge boson resulting from an L µ − L τ U(1) local gauge symmetry has been presented. Events containing four muons with an invariant mass near the standard model Z boson mass are analyzed, and the selection is further optimized to be sensitive to the events that may contain Z → Z µµ → 4µ decays. A data sample of proton-proton collisions at a center-of-mass energy of 13 TeV is used, corresponding to an

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus  [19] CMS Collaboration, "Measurements of the pp→ZZ production cross section and the Z→ 4 branching fraction, and constraints on anomalous triple gauge couplings at √ s = 13 TeV", Eur. Phys. J C 78 (2018) 165, doi:10.1140/epjc/s10052-018-5567-9, arXiv:1709.08601.