Search for $Z^{'} \rightarrow \mu^{+} \mu^{-}$ in the $L_\mu {-} L_\tau$ gauge-symmetric model at Belle

We search for a new gauge boson $Z'$ that couples only to heavy leptons and their corresponding neutrinos in the process $e^{+} e^{-} \rightarrow Z'(\rightarrow \mu^{+}\mu^{-}) \mu^{+}\mu^{-}$, using a 643 fb$^{-1}$ data sample collected by the Belle experiment at or near the $\Upsilon(1S,2S,3S,4S,5S)$ resonances at the KEKB collider. For the first time, effects due to initial state radiation are used in estimating the detection efficiency. No signal is observed in the mass range of 0.212 \ -- 10 GeV/${\it c}^2$ and we set an upper limit on the coupling strengh, $g'$, constraining $Z'$ as a possible contributor to the anomalous magnetic dipole moment of muon.

We search for a new gauge boson Z that couples only to heavy leptons and their corresponding neutrinos in the process e + e − → Z (→ µ + µ − )µ + µ − , using a 643 fb −1 data sample collected by the Belle experiment at or near the Υ(1S, 2S, 3S, 4S, 5S) resonances at the KEKB collider. While previous searches for Z did a data-based estimation of the initial state radiation effect, our search for the Z is the first to include effects due to initial state radiation in the signal simulated samples used in estimating the detection efficiency. No signal is observed in the Z mass range of 0.212 -10 GeV/c 2 and we set an upper limit on the coupling strength, g , constraining the possible Z contribution to the anomalous magnetic dipole moment of the muon.

I. INTRODUCTION
The lack of evidence for a Weakly Interacting Massive Particle by underground experiments [1,2], and the absence of supersymmetric particle signals at the LHC [3][4][5], suggest that dark matter might be composite and/or light. This gives rise to dark sector models [6][7][8][9][10][11][12][13][14][15][16] that introduce a zoology of dark particles which do not interact directly via Standard Model (SM) forces, but can interact by dark sector forces via new mediators and therefore only indirectly with (SM) particles, and could have masses between 1 MeV/c 2 and 10 GeV/c 2 .
Discrepancies observed at low-energy measurements [17,18] have fueled new precision studies. Within this context, the anomalous magnetic moment of the muon, (g − 2) µ , is one of the most precisely measured quantities in particle physics, where the difference between the experimental value and the SM prediction [19] is about 4.2σ [20]. This discrepancy might be a sign of new physics and has led to a variety of attempts to create physics models involving the leptonic sector of the SM [21][22][23][24].
These attempts include the set of SM extensions which add a new U (1) gauge boson (Z ) coupled to the difference between lepton family numbers, L i where i = e, µ and τ [22]. The electron number differences have been well constrained by measurements performed at e + e − colliders [25,26] and will not be discussed here. In this study we present a search for the gauge boson coupled to the L µ − L τ difference.
The partial widths for the Z decay to leptons [27,28] are given by: where g is the L µ − L τ coupling strength, and θ(M Z − 2M ) is a step function, and For M Z M the branching fraction to one neutrino flavor is half of that to a charged lepton. This is due to the fact that the Z boson only couples to left-handed neutrinos, but couples to both left-and right-handed charged leptons.
The visible branching fraction to muons is: which is identical to B(Z → τ τ ) except for the replacement of the decay width with the appropriate decay channel.
We search for the Z of an L µ − L τ model via the decay Z → µ + µ − . In this model, the Z only couples to the second and third generation of leptons (µ, τ ) and their neutrinos. We search for four-muon events in the reaction depicted in Fig. 1 Feynman diagram for the main production channel of the Z in e + e − colliders.
process is followed by Z radiation from a muon, and then, the Z decays to µ + µ − . In addition to its possible contribution to the (g − 2) µ anomaly, the effects of a Z have been searched for in other scenarios. It could be a source of an increase in neutrino trident production ν µ N → N ν µ µ + µ − [23]. No increase has been observed, and a limit was set for the Z parameter space. It could also work as an indirect channel to sterile neutrino dark matter [24], and could provide predictions for the neutrino mass-mixing matrix [29][30][31].
Recently, the Belle II collaboration published the search result with Z → νν decay using a 276 pb −1 luminosity data [32]. No Z signature was found so an upper limit of the parameter space of this decay mode was set. Previously BABAR searched for the Z with e + e − → Z (→ µ + µ − )µ + µ − using a 514 fb −1 luminosity data and since no Z signature was found the most stringent upper limits as a function of Z mass [33] was set. In this paper, we present a search for the same Z model in the full available Belle data sample.
The Belle detector surrounds the interaction point of KEKB. It is a large-solid-angle magnetic spectrometer consisting of a silicon vertex detector, a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters, a barrel-like arrangement of timeof-flight scintillation counters, and an electromagnetic calorimeter (ECL) comprised of CsI(TI) crystals located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux return located outside of the coil is instrumented with resistive plate chambers to detect K 0 L mesons and identify muons (KLM). The signal events are guaranteed to pass the trigger with nearly full efficiency because the muonic Z decay topology features more than three charged tracks. In addition, the large radius of the CDC (880 mm) [38] allows an excellent mass resolution and muon detection efficiency in Belle.

III. SELECTION CRITERIA
The selection is optimized based on the validation sample as well as a Monte Carlo (MC) simulation done in two steps. First, signal events are generated for different Z mass hypotheses using WHIZARD [39], which takes into account the Initial State Radiation (ISR) as well as the Final State Radiation (FSR) at the Υ(4S) centerof-mass energy. WHIZARD also has an option to generate events without radiative corrections. Then, the detector response to these events is simulated using GEANT3 [40]. There were 54 mass hypotheses generated for each of the Z MC samples from m Z = 212 MeV/c 2 to m Z = 1.015 GeV/c 2 in 100 MeV/c 2 steps, and subsequently in 200 MeV/c 2 steps up to m Z = 10.00 GeV/c 2 . The change in the steps is due to the behavior of the detection efficiency observed in the analysis.
The irreducible background, e + e − → µ + µ − µ + µ − is studied with an MC sample corresponding to a luminosity of 336 fb −1 generated with Diag36 [41] at Υ(4S) center-of-mass energy, Diag36 generates events without ISR corrections (non-ISR). There is no event generator available for the QED 4-lepton final state with radiative correction. In addition, other leptonic and hadronic background sources, such as e + e − → e + e − e + e − and e + e − → π + π − J/ψ(→ µ + µ − ) were studied through MC samples, and found to give negligible or no contributions after the application of the selection criteria.
We select events with two pairs of oppositely charged tracks in the final state. To ensure these tracks originate from the interaction point, their transverse and longitudinal impact parameters must be less than 0.2 and 1.5 cm, respectively. At least two tracks are required to have a muon likelihood ratio, Lµ Lµ+L K +Lπ , greater than 0.1. The value of L µ depends on the difference between the expected and actual muon penetration of the track in the KLM, and the distance between its KLM hits and the extrapolation of the track from the CDC. The efficiency for a track to be identified as a muon is about 95% for momenta between 1 to 3 GeV/c, and slightly lower momenta below 1 GeV/c. In addition, a hadron veto is applied. The muon candidate must not have a likelihood ratio corresponding to a pion, kaon, or a proton. This is implemented by comparing the likelihood ratio of two particles (proton and kaon, kaon and pion, and proton and pion) as P (i|j) = Li Li+Lj , where L i is the likelihood product from three detectors (ACC, TOF, and CDC). A pion is defined as P (K|π) < 0.4 and P (p|π) < 0.4. A kaon is defined as P (p|K) < 0.4 and P (K|i) > 0.6. A proton is defined as P (p|K) > 0.6 and P (p|π) > 0.6.
To suppress the background due to neutral particles, the sum of ECL clusters unrelated to any charged tracks with energy greater than 30 MeV is required to be less than 200 MeV. Additionally, the visible energy, E vis , calculated from the four muons must be consistent with the center-of-mass energy, E CMS , so that |E CMS − E vis | < 500 MeV. A kinematic fit based on the least square method for the final state is carried out under the constraint that the four-momentum of the final state be compatible with the initial e + e − system. The chi-squared is minimized by Lagrange multipliers, they issue a set of non-linear equations that are solved using the multi-dimensional Newton-Raphson method. As a result, the reconstructed Z mass resolution is improved.
As there are four possible combinations of oppositely charged muons in the final state, all four possible combinations correspond to four Z candidates counted per event. To improve the sensitivity in the low Z mass region, we introduce a reduced mass, defined as where m µµ is the invariant masses of Z candidate and m 2 µ,PDG is the muon nominal mass [42]. The m R distribution is smoother than the invariant mass distribution around the Z mass close to the dimuon threshold.
The Z reduced mass distributions for data and MC are compared in Fig. 2. Although the normalization of the data is almost 70% of the background level, determined by a fit to a constant probability density function (pdf) as shown in Fig. 2 (bottom). This difference arises due to the ISR effect, which is not simulated in the background MC sample.
We veto the reduced mass distribution around the J/ψ mass, 3.05 < m R < 3.13 GeV/c 2 , as its muonic decay can mimic a signal. This was not necessary around the ψ(2S) mass since the ψ(2S) decay into muons is negligible compared to the main background.

IV. RESULTS
We perform a binned maximum-likelihood fit to the reduced mass distribution with the range of m Z ± 25σ Z . The fit is repeated 9788 times with a different Z mass hypotheses in steps of 1 MeV/c 2 from 0 to 9787 MeV/c 2 . The Z resolution starts from less than 1 MeV/c 2 at the dimuon mass threshold increasing until 5.5 GeV/c 2 where it is valued at 6 MeV/c 2 then it starts decreasing until 9.5 GeV/c 2 where it is valued at ∼3 MeV/c 2 . The step is set around half of the width of the reduced mass distribution for MC generated signal.
The signal m R distribution is modeled as a sum of two Crystal Ball [43] functions with a common mean. The shape parameters as a function of the m R are determined with signal MC samples while the normalization is floated in the fit. The width is calibrated using J/ψ → µ + µ − events in the veto region. The background is modeled with a third-order polynomial which is the lowest order (Top) Reduced mass, mR, distributions. Red points represent the data after all selection criteria is applied. Black squares represent the non-ISR MC expectation for the e + e − → µ + µ − µ + µ − (Diag36) [41] scaled to the data luminosity. (Bottom) The ratio between data and the non-ISR e + e − → µ + µ − µ + µ − MC expectation. The red line represents a fit of a 1st order polynomial where its constant term is 0.700 ± 0.003 and the slope term is negligible. (Both) Black shaded region at 3.1 GeV/c 2 represents the J/ψ region which is not used in this analysis.
function that can fit the e + e − → µ + µ − µ + µ − background well. Background normalization and shape parameters are floated in the fit.
The efficiency is determined using a fit to the MC signal samples with different mass hypotheses. It is the result of the integration of the fit function over m Z ± 3σ Z , where σ Z is the Z mass resolution. This efficiency is interpolated between the different discrete mass hypotheses.
This procedure is done identically for non-ISR MC samples where the m R distribution is also modeled as a sum of two Crystal Ball functions with a common mean, however, for the non-ISR MC samples the parametrization of the pdf is different than for the ISR case. Comparing ISR and non-ISR detection efficiencies is key to understand the gap between data and MC background on Fig. 2. Fig. 3 shows efficiencies as a function of reduced mass. It is clear that the detection efficiency increases with increasing Z mass up to 6 GeV/c 2 and then it decreases. This behavior is due to the muon detection efficiency in the KLM, which has a threshold momentum of 600 MeV/c reaching maximum at 1 GeV/c then flattens for even higher values.
Systematic uncertainties arise from luminosity, track identification, muon identification and fitting bias. The luminosity uncertainty is 1.4% and is measured using Bhabha and two-photon events. The track identification uncertainty is 0.35% per charged track, or 1.4% in this analysis, and is determined by comparing the track finding efficiency of partially and fully reconstructed D * + → D 0 (→ K − π + )π + decays. A muon identification uncertainty of 1.15% is determined from the change in event yields while varying the muon likelihood ratio criterion from 0.1 to 0.2. With muon likelihood cuts there is also a systematic error to be considered on the detection efficiency calculated through MC signal samples. This error is calculated comparing γγ → µ + µ − data and MC samples. Due to the large number of these events it is possible to map the dependency between momentum, muon likelihood ratio and error rate. Comparing our MC signal calculated detection efficiency with and without this correction gives a 1% difference. Finally, a correction from the hadron veto is implemented on the MC samples. This correction factor is also of 1% and it is obtained by comparing MC samples with and without the hadron veto.
The effect of fitting bias is investigated using a bootstrap study to check whether allowing third-order polynomial components to float in the fit end up inducing a bias on the yield extracted. For each mass scan, this study is done by varying the data with a Poisson distribution, varying each individual bin of the histogram. This changed data set is then injected with a signal of yield corresponding to a Poisson distribution of the upper limit on the number of observed events and a distribution following its pdf. This reconstructed ensemble is then fitted in the same way as the data. The newly ex-tracted yield, N sig , is then compared to the true number of events injected, N true sig , divided by the uncertainty in the newly yield extracted, σ Nsig , as (N true sig − N sig )/σ Nsig . This procedure is repeated 1000 times for each mass scan. We find that the extracted yield and its uncertainty are systematically overestimated by 3% and 4%, respectively. These biases are accordingly taken into account in the Z scan and g upper limit calculation by correcting the yield extracted and the error on the yield extracted. They correspond to N cor sig = N sig (1 + b) and N errcor sig = N err sig × b err , where b stands for bias and the variables with a superscript err are related to the error on the yield.
The significance of each possible Z candidate is evaluated as where sign(N obs ) is the sign of the number of observed events and L S+B /L B is the ratio between the maximum likelihoods of the fits with signal plus background hypothesis (L S+B ) and background only hypothesis (L B ). The distribution of significances is shown Fig. 4. The largest local significance observed in an excess (deficit) is 3.7σ (3.5σ) around m Z = 3.3 GeV/c 2 (3.1 GeV/c 2 ), in Fig. 5. After incorporating the lookelsewhere-effect the global significance for the excess becomes 2.23σ.
Since no fit resulted in a global significance of at least 5σ, we set upper limits on the coupling strength g as a function of m Z . A Bayesian method [44] is used to estimate the 90% credibility level (C.L.) upper limit on the number of observed signal events, N obs . A flat prior is assumed for the signal yield and two nuisance parameters are added one for the signal yield and another for  The results are shown in Fig. 6. Using the calculated detection efficiency as shown in Fig. 3, the branching fraction from Eq.(3) and the Belle luminosity (L) of 643 fb −1 , the 90% upper limit on the Born e + e − → Z µ + µ − cross-section is obtained using: where N obs , ISR , B, (1 + δ) and |1 − Π(s)| 2 are the upper limit on the yield extracted from the data scan as shown in Fig. 6, the ISR signal MC sample based detection efficiency, the branching fraction from Eq.(3), the ISR correction factor, and the vacuum polarization factor, respectively. In order to test the ISR and the vaccum polarization effects, we check the ratio between the number of observed signal N S obs and the number of simulated signal events N S MC can be written as: where ISR ( non−ISR ) is the detection efficiency obtained by the ISR (non-ISR) signal MC. Since the cross-section with the ISR and vacuum polarization corrections (σ V ) is related to the Born cross section by 46] the ratio, Eq.(6), becomes: As the ISR and vacuum polarization corrections are common for the signal and the e + e − → µ + µ − µ + µ − (4µ) background process, one can expect that the ratio, Eq.(6), is the same for the signal and the 4µ background: Checking the consistency of the efficiency and ISR correction factors can be carried out by the 4µ MC background and data. From Fig. 2, we observe the ratio between data and the MC expectation for the 4µ process to be: The 90% C.L. upper limits on Born cross-section as a function of m Z are calculated and shown in Fig. 7.

A. Limits on the Coupling Strength g
With a Born theoretical cross-section σ th ( √ s), for a given m Z , at √ s and the coupling g , the expected number of signal events for data samples used in this analysis is given as: With Eq.(8), the 90% C.L. upper limit on g corresponding to N exp = N obs , is calculated and shown in Fig. 8. The result excludes most of the Z parameter space that could be related to the updated (g − 2) µ region, from the Muon (g − 2) experiment [17,20]. Also shown in Fig. 8 are comparisons with the CHARM-II experiment, the first measurement of the neutrino trident production [47], the reinterpretation of the Columbia-Chicago-Fermilab-Rochester (CCFR) results [23,48] and the first Z → µ + µ − search done by BABAR [33].

V. CONCLUSION
In summary, we report a search for a new gauge boson Z in the L µ − L τ model with the on-shell production of e + e − → Z µ + µ − , followed by Z → µ + µ − . This is the first search with the ISR effect directly included in the MC signal sample, while previous searches did a datadriven estimation of the ISR effect. Since no significant excess is observed, the upper limit on the coupling is set and the Z parameter space constraint is improved.
This result specifically improves the previous g upper limit between 2 and 8.4 GeV/c 2 .
The Z mass region lighter than the dimuon threshold, does not have any constraints but in the future, Belle II will be able to perform a more stringent test for the region [49][50][51].  [33], the light gray shaded area is the result of CCFR, and the light purple hashed area the result of CHARM-II over the Z parameter space [23,48], and the green region indicates the values of the Z coupling needed to explain (g−2)µ suggested by the Muon g−2 collaboration [17,20].

VI. ACKNOWLEDGMENTS
We thank B. Shuve for providing the models for MadGraph5 and the branching fractions for Z . Our gratitude goes to K. Mawatari for showing us the limitations of MadGraph5 when simulating ISR events and to J. Reuter for explaining how to use WHIZARD for generating ISR Z events. We also thank T. Shimomura for the enlightening discussions about Z .
T. C. is supported by the Japan Society for the Promotion of Science (JSPS) Grant No. 20H05858 and A. I. is supported by Grant No. 16H02176 and 22H00144.
We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for strong computing support; and the National Institute of Informatics, and Science Information NETwork 5 (SINET5) for valuable network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the JSPS including Grant No. 20H05850, and the Tau-Lepton Physics Research Center of Nagoya University; the Australian Research Council including grants DP180102629, DP170102389, DP170102204, DP150103061, FT130100303; Austrian Federal Min-