Determination of the pseudoscalar decay constant $f_{D_s^+}$ via $D_s^+\to\mu^+\nu_\mu$

Using a $3.19~\mathrm{fb}^{-1}$ data sample collected at an $e^+e^-$ center-of-mass energy of $E_{\rm cm}=4.178$ GeV with the BESIII detector, we measure the branching fraction of the leptonic decay $D_s^+\to\mu^+\nu_\mu$ to be $\mathcal{B}_{D_s^+\to\mu^+\nu_\mu}=(5.49\pm0.16_{\rm stat.}\pm0.15_{\rm syst.})\times10^{-3}$. Combining our branching fraction with the masses of the $D_s^+$ and $\mu^+$ and the lifetime of the $D_s^+$, we determine $f_{D_s^+}|V_{cs}|=246.2\pm3.6_{\rm stat.}\pm3.5_{\rm syst.}~\mathrm{MeV}$. Using the $c\to s$ quark mixing matrix element $|V_{cs}|$ determined from a global standard model fit, we evaluate the $D_s^+$ decay constant $f_{D_s^+}=252.9\pm3.7_{\rm stat.}\pm3.6_{\rm syst.}$\,MeV. Alternatively, using the value of $f_{D_s^+}$ calculated by lattice quantum chromodynamics, we find $|V_{cs}| = 0.985\pm0.014_{\rm stat.}\pm0.014_{\rm syst.}$. These values of $\mathcal{B}_{D_s^+\to\mu^+\nu_\mu}$, $f_{D_s^+}|V_{cs}|$, $f_{D_s^+}$ and $|V_{cs}|$ are each the most precise results to date.

Using a 3.19 fb −1 data sample collected at an e + e − center-of-mass energy of Ecm = 4.178 GeV with the BESIII detector, we measure the branching fraction of the leptonic decay D + s → µ + νµ to be B D + s →µ + νµ = (5.49 ± 0.16stat. ± 0.15syst.) × 10 −3 . Combining our branching fraction with the masses of the D + s and µ + and the lifetime of the D + s , we determine f D + s |Vcs| = 246.2±3.6stat. ±3.5syst. MeV. Using the c → s quark mixing matrix element |Vcs| determined from a global standard model fit, we evaluate the D + s decay constant f D + s = 252.9 ± 3.7stat. ± 3.6syst. MeV. Alternatively, using the value of f D + s calculated by lattice quantum chromodynamics, we find |Vcs| = 0.985 ± 0.014stat. ± 0.014syst.. These values of B D + s →µ + νµ , f D + s |Vcs|, f D + s and |Vcs| are each the most precise results to date. The leptonic decay D + s → ℓ + ν ℓ (ℓ = e, µ or τ ) offers a unique window into both strong and weak effects in the charm quark sector. In the standard model (SM), the partial width of the decay D + s → ℓ + ν ℓ can be written where f D + s is the D + s decay constant, |V cs | is the c → s Cabibbo-Kobayashi-Maskawa (CKM) matrix element, G F is the Fermi coupling constant, m ℓ is the lepton mass, and m D + s is the D + s mass. In recent years, much progress has been achieved in the measurements of f D + s and |V cs | with D + s → ℓ + ν ℓ decays at the CLEO-c [2][3][4], BaBar [5], Belle [6] and BESIII [7] experiments. However, compared to the precision of the most accurate lattice quantum chromodynamics (LQCD) calculation of f D + s [8], the accuracy of the measurements is still limited. Improved measurements of f D + s and |V cs | are critical to calibrate various theoretical calculations of f D + s , such as those from quenched and unquenched LQCD, QCD sum rules, etc., and to test the unitarity of the quark mixing matrix with better precision.
In the SM, the ratio of the branching fraction (BF) of D + s → τ + ν τ over that of D + s → µ + ν µ is predicted to be 9.74 with negligible uncertainty and the BFs of D + s → µ + ν µ and D − s → µ −ν µ decays are expected to be the same. However, hints of lepton flavor universality (LFU) violation in semileptonic B decays were recently reported at BaBar, LHCb and Belle [37][38][39][40][41]. It has been argued that new physics mechanisms, such as a two-Higgs-doublet model with the mediation of charged Higgs bosons [42,43] or a Seesaw mechanism due to lepton mixing with Majorana neutrinos [44], may cause LFU or CP violation. Tests of LFU and searches for CP violation in D + s → ℓ + ν ℓ decays are therefore important tests of the SM.
In this Letter, we present an experimental study of the leptonic decay D + s → µ + ν µ [45] by analyzing a 3.19 fb −1 data sample collected with the BESIII detector at an e + e − center-of-mass energy of E cm = 4.178 GeV. At this energy, D + s mesons are produced mainly through the process e + e − → D + s D * − s + c.c. In an event where a D − s meson (called a single-tag (ST) D − s meson) is fully reconstructed, one can then search for a γ or π 0 and a D + s meson in the recoiling system (called a double-tag (DT) event).
Details about the design and performance of the BE-SIII detector are given in Ref. [46]. The endcap time-offlight (TOF) system was upgraded with multi-gap resistive plate chamber technology and now has a time resolution of 60 ps [47,48]. Monte Carlo (MC) events are generated with a geant4-based [49] detector simulation software package [50], which includes both the geometrical description of the detector and the detector's response. An inclusive MC sample is produced at √ s = 4.178 GeV, which includes all open charm processes, initial state radiation (ISR) production of the ψ(3770), ψ(3686) and J/ψ, and qq (q = u, d, s) continuum processes, along with Bhabha scattering, µ + µ − , τ + τ − and γγ events. The open charm processes are generated using ConExc [51]. The effects of ISR [52] and final state radiation (FSR) [53] are considered. The decay modes with known BF are generated using EvtGen [54] and the other modes are generated using LundCharm [55].
All charged tracks except for those from K 0 S decays must originate from the interaction point (IP) with a distance of closest approach less than 1 cm in the transverse plane and less than 10 cm along the z axis. The polar angle θ of each track defined with respect to the positron beam must satisfy | cos θ| < 0.93. Measurements of the specific ionization energy loss (dE/dx) in the main drift chamber and the TOF are combined and used for particle identification (PID) by forming confidence levels for pion and kaon hypotheses (CL π , CL K ). Kaon (pion) candidates are required to satisfy CL K(π) > CL π(K) .
To select K 0 S candidates, pairs of oppositely charged tracks with distances of closest approach to the IP less than 20 cm along the z axis are assigned as π + π − without PID requirements. These π + π − combinations are required to have an invariant mass within ±12 MeV of the nominal K 0 S mass [56] and have a decay length of the reconstructed K 0 S larger than 2σ of the vertex resolution away from the IP. The π 0 and η mesons are reconstructed via γγ decays. It is required that each electromagnetic shower starts within 700 ns of the event start time and its energy is greater than 25 (50) MeV in the barrel (endcap) region of the electromagnetic calorimeter (EMC) [46]. The opening angle between the shower and the nearest charged track has to be greater than 10 • . The γγ combinations with an invariant mass M γγ ∈ (0.115, 0.150) and (0.50, 0.57) GeV/c 2 are regarded as π 0 and η mesons, respectively. A kinematic fit is performed to constrain M γγ to the π 0 or η nominal mass [56]. The η candidates for the ηπ − ST channel are also reconstructed via π 0 π + π − candidates with an invariant mass within (0.53, 0.57) GeV/c 2 . The η ′ mesons are reconstructed via two decay modes, ηπ + π − and γρ 0 , whose invariant masses are required to be within (0.946, 0.970) and (0.940, 0.976) GeV/c 2 , respectively. In addition, the minimum energy of the γ from η ′ → γρ 0 decays must be greater than 0.1 GeV. The ρ 0 and ρ + mesons are reconstructed from π + π − and π + π 0 candidates, whose invariant masses are required to be larger than 0.5 GeV/c 2 and within (0.67, 0.87) GeV/c 2 , respectively.
The momentum of any pion not originating from a K 0 S , η, or η ′ decay is required to be greater than 0.1 GeV/c to reject soft pions from D * decays. For π + π − π − and K − π + π − combinations, the dominant peaking backgrounds from K 0 S π − and K 0 S K − events are rejected by requiring the invariant mass of any π + π − combination be more than ±0.03 GeV/c 2 away from the nominal K 0 S mass [56].
To suppress non-D + s D * − s events, the beam-constrained mass of the ST D − s candidate is required to be within (2.010, 2.073) GeV/c 2 , where p D − s is the momentum of the ST D − s candidate. In each event, we only keep the candidate with the D − s recoil mass closest to the nominal D * + s mass [56] per tag mode per charge. Figure 1 shows the invariant mass spectra of the accepted ST candidate events for 14 ST modes. The ST yield for each tag mode is obtained by a fit to the corresponding invariant mass spectrum. The signal is described by the MC-simulated shape convolved with a Gaussian function representing the resolution difference between data and the MC simulation. For the tag mode D − s → K 0 S K − , the peaking background from D − → K 0 S π − is described by the MC-simulated shape and then smeared with the same Gaussian function used in the signal shape with its size as a free parameter. The non-peaking background is modeled by a secondor third-order Chebychev polynomial function. Studies of the inclusive MC sample validate this parametrisation of the background shape. The fit results on these invariant mass spectra are shown in Fig. 1. The events in the signal regions are kept for further analysis. The total ST yield in data is N tot ST = 388660 ± 2592. The selected ST D − s meson may be directly produced in e + e − annihilation or produced by a D * − s → γ(π 0 )D − s decay. To separate D + s → µ + ν µ signals from combinatorial backgrounds, we define two kinematic variables and Here E miss ≡ | p miss | 2 + m 2 D + s and p miss ≡ − p tag − p γ(π 0 ) , where E i and p i (i = µ, γ(π 0 ) or tag) denote the energy and momentum of the muon, γ(π 0 ) or ST D − s , respectively. We loop over all remaining γ or π 0 candidates and choose the one that gives a minimum for |∆E|. The events with ∆E ∈ (−0.05, 0.10) GeV are accepted. The muon candidate is required to have an opposite charge to the ST D − s meson and a deposited energy in the EMC within (0.0, 0.3) GeV. It must also satisfy a two dimensional (| cos θ µ | and momentum p µ ) requirement on the hit depth (D) in the muon counter, as explained in Ref. [57]. To suppress the backgrounds with extra photon(s), the maximum energy of the unused showers in the DT selection (E max extra γ ) is required to be less than 0.4 GeV and no additional charged track that satisfies the charged track selection criteria is allowed. To improve the M 2 miss resolution, the candidate tracks, plus the missing neutrino, are subjected to a 4constraint kinematic fit requiring energy and momentum conservation. In addition, the invariant masses of the two D s mesons are constrained to the nominal D s mass, the invariant mass of the D − s γ(π 0 ) or D + s γ(π 0 ) combination is constrained to the nominal D * s mass, and the combination with the smaller χ 2 is kept. Figure 2 shows the M 2 miss distribution for the accepted DT candidate events.
To extract the DT yield, an unbinned fit is performed to the M 2 miss distribution. In the fit, the background events are classified into three categories: events with correctly reconstructed D − s and µ + but an unmatched γ(π 0 ) from the D * − s (BKGI), events with a correctly reconstructed ST D − s but misidentified µ + (BKGII), and other events with a misreconstructed ST D − s (BKGIII). The signal and BKGI shapes are modeled with MC simulation. The signal shape is convolved with a Gaussian function with its mean and width as free parameters. The ratio of the signal yield over the BKGI yield is constrained to the value determined with the signal MC events. The size and shape of the BKGII and BKGIII components are fixed by analyzing the inclusive MC sample. From the fit to the M 2 miss distribution, as shown in  The efficiencies for reconstructing the DT candidate events are determined with an exclusive MC sample of e + e − → D * + s D − s , where the D − s decays to each tag mode and the D + s decays to µ + ν µ . Dividing them by the ST efficiencies determined with the inclusive MC sample yields the corresponding efficiencies of the γ(π 0 )µ + ν µ reconstruction. The averaged efficiency of finding γ(π 0 )µ + ν µ is (52.67 ± 0.19)% as determined from where N i ST , ε i ST , and ε i DT are the ST yield, ST efficiency and DT efficiency in the i-th ST mode, respectively. The factor f PID µ = 0.897 accounts for the difference between data and MC simulation in the two dimensional (momentum and cos θ) muon PID efficiencies. This factor is estimated using e + e − → γµ + µ − events, and is further validated using e + e − → γ ISR ψ(3686), ψ(3686) → π + π − J/ψ, J/ψ → µ + µ − events. The absolute BF of D + s → µ + ν µ is then determined to be (5.49 ± 0.16 stat. ± 0.15 syst. ) × 10 −3 from where the radiative correction factor f rad cor = 0.99 is due to the contribution from D + s → γD * + s → γµ + ν µ [58], with D * + s as a virtual vector or axial-vector meson. This contribution is almost identical with our signal process for low energy radiated photons.
The systematic uncertainties in the BF measurement are estimated relative to the measured BF and are described below.
For uncertainties in the event selection criteria, the µ + tracking and PID efficiencies are studied with e + e − → γµ + µ − events. After considering the differences of the two dimensional weighted efficiencies between data and MC simulation and correcting the detection efficiency by f PID µ , we assign 0.5% and 0.8% as the uncertainties in µ + tracking and PID efficiencies, respectively. The photon reconstruction efficiency has been previously studied with J/ψ → π + π − π 0 decays [59]. The uncertainty of finding γ(π 0 ) is weighted according to the BFs of D * + s → γD + s and D * + s → π 0 D + s [56] and assigned to be 1.0%. The efficiencies for the requirements of E max extra γ and no extra good charged track are studied with a DT hadronic sample. The systematic uncertainties are taken to be 0.3% and 0.9% considering the efficiency differences between data and MC simulation, respectively. The uncertainty of the ∆E requirement is estimated by varying the signal region by ±0.01 GeV, and the maximum change of the BF, 0.5%, is taken as the systematic uncertainty.
To determine the uncertainty in the M 2 miss fit, we change the fit range by ±0.02 GeV 2 /c 4 , and the largest change of the BF is 0.6%. We change the signal shape by varying the γ(π 0 ) match requirement and the maximum change is 0.2%. Two sources of uncertainty in the background estimation are considered. The effect of the background shape is obtained to be 0.2% by shifting the number of the main components of BKGII by ±1σ of the uncertainties of the corresponding BFs [56], and varying the relative fraction of the main components of BKGII by 50%. The effect of the fixed number of the BKGII and BKGIII is estimated to be 0.5% by varying the nominal numbers by ±1σ of their uncertainties. To size up the uncertainty in the fixed ratio of signal and BKGI, we perform an alternative fit to the M 2 miss distribution of data without constraining the ratio of signal and BKGI. The change in the DT yield, 1.1%, is assigned as the relevant uncertainty.
The uncertainty in the number of ST D − s mesons is assigned to be 0.8% by examining the changes of the fit yields when varying the signal shape, background shape, bin size and fit range and considering the background fluctuation in the fit. The uncertainty due to the limited MC size is 0.4%. The uncertainty in the imperfect simulation of the FSR effect is estimated as 0.4% by varying the amount of FSR photons in signal MC events [53]. The uncertainty due to the quoted BFs of D * − s subdecays from the particle data group (PDG) [56] is examined by varying each subdecay BF by ±1σ. The efficiency change is found to be 0.4% and is taken as the associated uncertainty. The uncertainty in the radiative correction is assigned to be 1.0%, which is taken as 100% of its central value from theoretical calculation [58]. The ST efficiencies in the inclusive and signal MC samples are slightly different with each other due to different track multiplicities in these two environments. This may cause incomplete cancellation of the uncertainties of the ST efficiencies. The associated uncertainty is assigned as 0.6%, by taking into account the differences of the efficiencies of tracking/PID of K ± and π ± , as well as the selections of neutral particles between data and MC simulation in different environments. The total systematic uncertainty is determined to be 2.7% by adding all the uncertainties in quadrature.
Combining our BF with the world average values of G F , m µ , m D + s and the lifetime of D + s [56] in Eq. (1)  Combining the obtained f D + s |V cs | and its counterpart f D + |V cd | measured in our previous work [65], along with |V cs | and |V cd | from the SM global fit [56], yields f D + s /f D + = 1.24 ± 0.04 stat. ± 0.02 syst. . It is consistent with the CLEO measurement [2] within 1σ and the LQCD calculation within 2σ [8]. Alternatively, with the input of f D + s /f D + calculated by LQCD [8], we obtain |V cd /V cs | 2 = 0.048 ± 0.003 stat. ± 0.002 syst. , which agrees with the one expected by |V cs | and |V cd | given by the CKMfitter within 2σ. Here, no systematic uncertainty is canceled since the two data samples were taken in different years.
Based on our result for B D + s →µ + νµ and those measured at the CLEO-c [2], BaBar [5] and Belle [6] experiments, along with a previous measurement at BESIII [7], the inverse-uncertainty weighted BF is determined to bē B D + s →µ + νµ = (5.49 ± 0.17) × 10 −3 [66]. The ratio of B D + s →µ + νµ over the PDG value of B D + s →τ + ντ = (5.48 ± 0.23)% [56] is determined to be where the uncertainties in the tracking and PID efficiencies of the muon, the ST yields, the limited MC statistics, as well as the signal shape and fit range in M 2 miss fits for D + s and D − s have been studied separately and are not canceled.
In summary, by analyzing 3.19 fb −1 of e + e − collision data collected at E cm = 4.178 GeV with the BESIII detector, we have measured B(D + s → µ + ν µ ), the decay constant f D + s , and the CKM matrix element |V cs |. These are the most precise measurements to date, and are important to calibrate various theoretical calculations of f D + s and test the unitarity of the CKM matrix with better accuracy. We also search for LFU and CP violation in D + s → ℓ + ν ℓ decays, and no evidence is found.