Observation of $D^+_s\to \eta^\prime \mu^+\nu_\mu$, Precision Test of Lepton Flavor Universality with $D^+_s\to \eta^{(\prime)} \ell^+\nu_\ell$, and First Measurements of $D^+_s\to \eta^{(\prime)}\mu^+\nu_\mu$ Decay Dynamics

By analyzing 7.33 fb$^{-1}$ of $e^+e^-$ annihilation data collected at center-of-mass energies between 4.128 and 4.226 GeV with the BESIII detector, we report the observation of the semileptonic decay $D^+_s\to \eta^\prime \mu^+\nu_\mu$, with a statistical significance larger than 10$\sigma$, and the measurements of the $D_s^+ \to \eta\mu^+\nu_\mu$ and $D_s^+ \to \eta^\prime\mu^+\nu_\mu$ decay dynamics for the first time. The branching fractions of $D_s^+ \to \eta\mu^+\nu_\mu$ and $D_s^+ \to \eta^\prime\mu^+\nu_\mu$ are determined to be $(2.235\pm0.051_{\rm stat}\pm0.052_{\rm syst})\%$ and $(0.801\pm0.055_{\rm stat}\pm0.028_{\rm syst})\%$, respectively, with precision improved by factors of 6.0 and 6.6 compared to the previous best measurements. Combined with the results for the decays $D_s^+ \to \eta e^+\nu_e$ and $D_s^+ \to \eta^\prime e^+\nu_e$, the ratios of the decay widths are examined both inclusively and in several $\ell^+\nu_\ell$ four-momentum transfer ranges. No evidence for lepton flavor universality violation is found within the current statistics. The products of the hadronic form factors $f_{+,0}^{\eta^{(\prime)}}(0)$ and the $c\to s$ Cabibbo-Kobayashi-Maskawa matrix element $|V_{cs}|$ are determined. The results based on the two-parameter series expansion are $f^{\eta}_{+,0}(0)|V_{cs}| = 0.452\pm0.010_{\rm stat}\pm0.007_{\rm syst}$ and $f^{\eta^{\prime}}_{+,0}(0)|V_{cs}| = 0.504\pm0.037_{\rm stat}\pm0.012_{\rm syst}$, which help to constrain present models on $f_{+,0}^{\eta^{(\prime)}}(0)$. The forward-backward asymmetries are determined to be $\langle A_{\rm FB}^\eta\rangle=-0.059\pm0.031_{\rm stat}\pm0.005_{\rm syst}$ and $\langle A_{\rm FB}^{\eta^\prime}\rangle=-0.064\pm0.079_{\rm stat}\pm0.006_{\rm syst}$ for the first time, which are consistent with the theoretical calculation.

a Also at the Moscow Institute of Physics and Technology, Moscow 141700, Russia b Also at the Novosibirsk State University, Novosibirsk, 630090, Russia c Also at the NRC "Kurchatov Institute", PNPI, 188300, Gatchina, Russia i Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China j Also at Frontiers Science Center for Rare Isotopes, Lanzhou University, Lanzhou 730000, People's Republic of China k Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000, People's Republic of China l Also at the Department of Mathematical Sciences, IBA, Karachi, Pakistan By analyzing 7.33 fb −1 of e + e − annihilation data collected at center-of-mass energies between 4.128 and 4.226 GeV with the BESIII detector, we report the observation of the semileptonic decay D + s → η ′ µ + νµ, with a statistical significance larger than 10σ, and the measurements of the D + s → ηµ + νµ and D + s → η ′ µ + νµ decay dynamics for the first time.The branching fractions of D + s → ηµ + νµ and D + s → η ′ µ + νµ are determined to be (2.235 ± 0.051stat ± 0.052syst)% and (0.801 ± 0.055stat ± 0.028syst)%, respectively, with precision improved by factors of 6.0 and 6.6 compared to the previous best measurements.Combined with the results for the decays D + s → ηe + νe and D + s → η ′ e + νe, the ratios of the decay widths are examined both inclusively and in several ℓ + ν ℓ four-momentum transfer ranges.No evidence for lepton flavor universality violation is found within the current statistics.The products of the hadronic form factors f η (′) +,0 (0) and the c → s Cabibbo-Kobayashi-Maskawa matrix element |Vcs| are determined.The results based on the two-parameter series expansion are f η +,0 (0)|Vcs| = 0.452 ± 0.010stat ± 0.007syst and f η ′ +,0 (0)|Vcs| = 0.504 ± 0.037stat ± 0.012syst, which help to constrain present models on f η (′) +,0 (0).The forward-backward asymmetries are determined to be A η FB = −0.059±0.031stat±0.005syst and A η ′ FB = −0.064±0.079stat±0.006syst for the first time, which are consistent with the theoretical calculation.
The couplings between the three families of leptons and the gauge bosons are expected to be equal in the standard model (SM).This property is known as lepton flavor universality (LFU).In recent years, however, hints of tensions between experimental measurements and the SM predictions were reported in the semileptonic (SL) B decays [1], the anomalous magnetic moment of the muon [2,3], and the Cabibbo angle anomaly [4,5].For example, the measured branching fraction (BF) ratios [6][7][8][9][10][11][12] deviate from the SM predictions by 3.3σ [1].Although these tensions have been explained by various theoretical models [13][14][15][16][17][18][19][20][21][22][23][24], no definite conclusion is established yet.Precision tests of LFU in different SL decays of heavy mesons provide deeper insight into these anomalies.Possible LFU in the SL D + s decays is not yet well tested [25], due to poor knowledge of the semimuonic D + s decays.Reference [26] notes that there may indeed be observable LFU violation effects in the SL decays mediated via c → sℓ + ν ℓ .
Previously, only BESIII reported the BFs of D + s → η (′) µ + ν µ [51] with large uncertainties of 20% (51%) using 0.482 fb −1 of e + e − collision data taken at a center-ofmass energy E CM = 4.009 GeV.Using 7.33 fb −1 of e + e − collision data taken at E CM between 4.128 and 4.226 GeV with the BESIII detector, we report the first observation of D + s → η ′ µ + ν µ , and the BFs of D + s → η (′) µ + ν µ are determined with improved precision by about sixfold.Based on these, we test µ-e LFU in D + s → η (′) ℓ + ν ℓ decays in the full kinematic range and several q 2 intervals.By analyzing the D + s → η (′) µ + ν µ dynamics, we determine f η (′) 0 (q 2 ) and A FB for the first time.Charge-conjugate modes are implied throughout this Letter.
A description of the design and performance of the BESIII detector can be found in Ref. [52].About 83% of the data analyzed in this Letter profits from an upgrade of the end cap time-of-flight system with multigap resistive plate chambers with a time resolution of 60 ps [53,54].Monte Carlo (MC) simulated events are generated with a geant4-based [55] simulation software, which includes the geometric description [56] , and J/ψ], q q (q = u, d, s) continuum processes, along with Bhabha scattering, µ + µ − , τ + τ − , and γγ events.The open charm processes are generated using conexc [57].The effects of ISR and final state radiation are included.Signal MC samples of the SL decays D + s → η (′) µ + ν µ are simulated with the two-parameter series expansion model [58], with parameters obtained in this work.The input cross section of e + e − → D ± s D * ∓ s is taken from Ref. [59].In the MC generation, known particle decays are generated by evtgen [60] with the BFs taken from the Particle Data Group [30], and other modes are generated using lundcharm [61]. In , and η π 0 π + π − ρ − , are called the single-tag (ST) D − s .Those in which the ST D − s , the transition γ(π 0 ) of the D * ∓ s decay, and the signal decays of D + s → η (′) µ + ν µ are simultaneously reconstructed are called double-tag (DT) events.Based on these, we determine the BF of the signal decay by where ST and N DT are the total ST and DT yields in data summing over tag mode k, and ǫ γ(π 0 )sig is the effective signal efficiency of selecting γ(π 0 )η (′) µ + ν µ in the presence of ST D − s .The ǫ γ(π 0 )sig is the averaged efficiency of ǫ γsig and ǫ π 0 sig , and estimated by k

ST
, where ǫ k ST and ǫ k DT are the ST and DT efficiencies for the k-th tag mode, respectively.
For each tag mode, the ST yield is extracted from a fit to the corresponding invariant-mass spectrum of the ST candidates.The selection criteria for all ST candidates are the same as Ref. [41], where detailed description can be found.The total ST yield is N tot ST = (817.0± 3.4 stat ) × 10 3 .
In the presence of ST D − s , we select candidates for the transition γ (π 0 ) from D * + s decay and signal D + s → η (′) µ + ν µ among the unused particles recoiling against the ST D − s .In the signal decay, the η is reconstructed via η → γγ or η → π 0 π + π − decay, and the η ′ is reconstructed via η ′ → η γγ π + π − or γπ + π − decay.Particle identification (PID) for pions combines the specific ionization information in the multilayer drift chamber and the flight time in the time-of-flight system, and PID for muons further combine the energy deposited in the electromagnetic calorimeter (EMC).Likelihoods under various particle hypotheses (L i , i = e, π, µ, and K) are calculated.Charged tracks satisfying L π > 0.001, L π > L K are assigned as pion candidates, and satisfying L µ > 0.001, L µ > L e , L µ > L K , and E EMC ∈ (0.10, 0.28) GeV are assigned as muon candidates, where E EMC is the energy deposited in the EMC of muon candidates.The selection criteria for the transition γ (π 0 ) and η (′) are the same as those in Ref. [41].The energy and momentum of the missing neutrino of the signal SL decay are derived as E νµ ≡ E CM − Σ i E i and p νµ ≡ −Σ i p i , respectively, where E i and p i are the energy and momentum of the particle i, with i running over the ST D − s , transition γ (π 0 ), η (′) , and µ + .The yield of signal events is determined by a fit to the distribution of the kinematic variable To improve the MM 2 resolution, the candidate tracks, along with the missing neutrino, are subjected to a 3-constraint kinematic fit requiring energy and momentum conservation, constraining the invariantmass of each D ± s meson to the known D ± s mass [30], and constraining the invariant-mass of the D − s γ(π 0 ) or D + s γ(π 0 ) combination to the known D * ± s mass [30].The combination with the lowest χ 2 is kept.The χ 2 for D + s → η ′ γπ + π − µ + ν µ is required to be less than 30 to further suppress the non-D ± s D * ∓ s backgrounds.
To suppress the backgrounds, the energy of any unused shower (E max γ extra ) in an event is required to be less than 0.2 GeV.The DT candidates are vetoed if they contain any additional charged tracks (N char extra ) or π 0 reconstructed by two unused photons (N π 0 extra ).To further reject the peaking backgrounds from D + s → η (′) π + and D + s → η (′) π + π 0 , the invariant-masses of η (′) µ + , M η (′) µ + , are required to be less than 1.8 GeV/c 2 for both D + s → ηµ + ν µ and D + s → η ′ µ + ν µ , and the invariant-masses of η (′) ν µ , M η (′) νµ , are required to be greater than 0.97 GeV/c 2 and 1.27 GeV/c 2 for D + s → ηµ + ν µ and D + s → η ′ µ + ν µ , respectively.After imposing all above selection criteria, the resulting MM 2 distributions of the accepted candidates for different signal decay modes are exhibited in Fig. 1.For D + s → η (′) µ + ν µ , the signal yield is extracted from a simultaneous unbinned maximum-likelihood fit to the MM 2 spectra for the two η or η ′ reconstruction modes, with BFs constrained to be the same in the fit.The signal, peaking backgrounds of  Table 1.Signal efficiencies (ǫ γ(π 0 )sig ), signal yields (NDT), and obtained BFs (Bsig).Efficiencies are averaged by ǫγsig and ǫ π 0 sig , and include the BFs of the η (′) and D * ∓ s subdecays.Numbers in the first and second parentheses are the most significant digits of the statistical and systematic uncertainties, respectively.

Decay
ηµ The systematic uncertainties in the BF measurements are listed in Table 1 of Ref. [63] and are discussed below.
The uncertainty in the ST D − s yields is studied by examining the change of the ST D − s yields by varying the matched angle for signal shape and the order of Chebychev polynomial for background shape.
The uncertainties in the tracking or PID efficiencies of π ± from the secondary η (′) decay and µ + from primary D + s decays are studied with the control samples of e + e − → K + K − π + π − and e + e − → γµ + µ − , respectively.The uncertainty of the π 0 or η reconstruction is assigned by studying the control sample of e + e − → K + K − π + π − π 0 .The uncertainty in the reconstruction of transition γ(π 0 ) from D * + s is studied with the control sample of J/ψ → π 0 π + π − [64].The uncertainty from the selection of the transition γ(π 0 ) from D * + s with the smallest |∆E| method is estimated by using the control samples of D + s → K + K − π + and D + s → ηπ + π 0 .The uncertainties due to the signal model are estimated by comparing the DT efficiencies by varying the input hadronic FFs measured by ±1σ.The uncertainties due to the M η (′) µ + and M η (′) νµ requirements are estimated by using the DT events of D + s → η (′) e + ν e .The uncertainties from the η (′) reconstruction are estimated by analyzing the control sample of J/ψ → φη (′) .The systematic uncertainties in the MM 2 fit are studied by repeating the fits with different signal and background shapes.The uncertainty in the peaking background yield is propagated by varying its size by ±1σ of the corresponding BF [30].
The uncertainties due to different multiplicities of tag environments [65] are assigned by studying of data-MC efficiency differences.The uncertainty due to the finite MC statistics, which is dominated by that of the DT efficiency, is considered as a systematic uncertainty.

The uncertainties of the E max
γ extra , N π 0 extra , and N char extra requirements are analyzed with the DT events of D + s → η (′) π + (π 0 ) and D + s → η (′) e + ν e .The uncertainties of the χ 2 requirements for D + s → η ′ γπ + π − µ + ν µ are studied with the DT events of The uncertainties due to the quoted BFs of η (′) and D * + s decays are taken from Ref. [30].The correlated and uncorrelated systematic uncertainties between the two η (′) decay modes are summarized in top and bottom section of Table 1 in Ref [63].The combined systematic uncertainties are 2.3% for D + s → ηµ + ν µ and 3.5% for D + s → η ′ µ + ν µ , taking into account correlated and uncorrelated systematic uncertainties with the method described in Ref. [66].
The decay dynamics of D + s → η (′) µ + ν µ are investigated by dividing individual candidate events into m = 8(3) intervals of q 2 and performing a least-χ 2 FF fit to the measured (∆Γ i msr ) and theoretically expected (∆Γ i th ) partial decay rates in the i-th q 2 interval.The ∆Γ i th ≡ i dΓ dq 2 dq 2 relate to the hadronic FF via [33] dΓ where G F is the Fermi coupling constant [30], The hadronic FFs f η (′) + (q 2 ) are parameterized with a two-parameter series expansion [67].We fix the pole mass at the known D * + s mass [30].The similar formulas are applied for f η (′) 0 (q 2 ) but with one-parameter series expansion due to much less contribution, and with the pole mass replaced with m D * s0 (2317) + [30].
The ∆Γ i msr are determined by ∆Γ i msr = where τ D + s is the D + s meson lifetime [30,68] and DT is the corresponding produced signal yield.The observed signal yield (N j DT ) is obtained from a similar fit of the corresponding MM 2 distribution.The signal efficiency matrix (ǫ ij ) is determined via gen is the total signal yield produced in the j-th q 2 interval, N ij rec is the number of events generated in the j-th q 2 interval but reconstructed in the i-th q 2 interval, and k sums over all tag modes.Details of q 2 divisions, the weighted signal efficiency matrices, N i DT , N i prd , and ∆Γ i msr of different q 2 intervals for D + s → ηµ + ν µ and D + s → η ′ µ + ν µ are shown in Tables 2-5 of Ref. [63], respectively.
The statistical and systematic covariance matrices are constructed as Ref. [41], which are shown in Tables 6-7 of Ref. [63].The systematic covariance matrices include with those uncertainties from the BF measurements, along with the D + s lifetime [30,68].For each signal decay, we perform a simultaneous fit on the differential decay rates measured by the two η (′) sub-decays, where the two modes are constrained to have same parameters for the hadronic FF.Figures 2(a,b) show the fits to the differential decay rates.The obtained parameters of hadronic FFs are summarized in Table 2. Figure 3 shows the projections to the extracted f η (′) +,0 (q 2 ), as well as the comparison of f η (′) +,0 (q 2 ) and various theoretical calculations.
The forward-backward asymmetry parameter A FB is defined as , where θ ℓ is ) 0.0 0.5 1.0 Fig. 3. Projections of the fits on f η (′) +,0 (q 2 ).Red circles and blue triangles with error bars are data, the yellow bands are the ±1σ limits of fitted parameters, and the curves in different colors are from various theoretical calculations described in the legend.
the angle between the momentum of the µ + in the rest frame of the W -boson and the direction to the W -boson momentum in the rest frame of D + s .The measured A FB in various q 2 intervals are shown in Figs.2(e,f).The averaged A FB are determined to be −0.059± 0.031 stat ± 0.005 syst for D + s → ηµ + ν µ and −0.064 ± 0.079 stat ± 0.006 syst for D + s → η ′ µ + ν µ , which are consistent with theoretical predictions [29,33].
In summary, we present for the first time the observation of D + s → η ′ µ + ν µ with statistical significance greater than 10σ and the studies of D + s → η (′) µ + ν µ dynamics.The BFs of D + s → ηµ + ν µ and D + s → η ′ µ + ν µ are measured with precision improved by about sixfold over the previous best measurements [30].Combining with the BESIII measurements of D + s → η (′) e + ν e [41], we calculate the R η (′)  µ/e ratios in separate q 2 intervals and in the full range.No significant evidence of LFU violation is found with current statistics.
By analyzing the dynamics of these decays, we determine the hadronic FFs of f η (′) +,0 (0), |V cs |, the shapes of f η (′) +,0 (q 2 ), and forward-backward asymmetry parameters of A FB .The obtained f η (′) 0 (q 2 ) lineshapes offer crucial data to calibrate the q 2 dependent FFs from different theories for the first time.Unlike the comparable FFs in the SL decays mediated via c → de + ν e [69], the D + s → η FF measured in this work deviate from f D→K + (0) = 0.7327 ± 0.0049 obtained via D 0 → K − µ + ν µ [70] by more than 5σ.This rules out the expectation for comparable FFs for D 0(+) → Kµ + ν µ and D + s → ηµ + ν µ [47], but supports that the spectator play important role in FFs due to effects of confinement in the considered weak transitions and the SU(4) asymmetry breaking [71].
The authors thank Prof. Haibing Fu, Prof. Shan Cheng, Prof. Xianwei Kang, and Prof. Khlopov for helpful dis-cussions.The BESIII Collaboration thanks the staff of BEPCII and the IHEP computing center APPENDICES Table 1 summarizes the sources of the systematic uncertainties in the measurements of the branching fractions of D + s → η (′) µ + ν µ .
In this table, the contributions to the systematic uncertainties listed in the top part are treated as correlated, while those in the bottom part are treated as uncorrelated.
Tables 2 and 3 give the weighted efficiency matrices averaged over all fourteen ST modes for D + s → ηµ + ν µ and D + s → η ′ µ + ν µ , respectively.Tables 4 and  5 present the numbers of the reconstructed events N i DT in data obtained from the MM 2 fits, the numbers of the produced events N i prd , and the measured partial decay rate ∆Γ i msr in different q 2 intervals for D + s → ηµ + ν µ and D + s → η ′ µ + ν µ , respectively.
Table 6 summarizes the statistical correlation matrices and relative uncertainties for the measured partial decay rates in different q 2 intervals for D + s → ηµ + ν µ .Table 7 summarizes the systematic correlation matrices and relative uncertainties for the measured partial decay rates in different q 2 intervals for D + s → η ′ µ + ν µ .
Table 1.Relative systematic uncertainties (in %) on the measurements of the branching fractions of D + s → ηµ + νµ and D + s → η ′ µ + νµ.The top and the bottom sections are correlated and uncorrelated, respectively.The uncertainty in the uncorrelated π ± tracking is obtained as the square root of the quadratic difference of the total uncertainty in the π ± tracking and the correlated portion.The last row of combined uncertainties are total uncertainties of D + s → ηµ + νµ and D + s → η ′ µ + νµ after taking into account correlated and uncorrelated systematic uncertainties.Table 2.The efficiency matrices for D + s → ηµ + νµ averaged over all fourteen ST modes, where εij represents the efficiency of events produced in the j-th q 2 interval and reconstructed in the i-th q 2 interval.

Fig. 1 .
Fig. 1.Simultaneous fits to the MM 2 distributions of the accepted candidates for D + s → η (′) µ + νµ.Points with error bars are data, solid curves are the best fits, differences between dashed and dotted curves are the D + s → η (′) π + (π 0 ) peaking backgrounds, and dotted curves are the other backgrounds.

Fig. 2 .
Fig. 2. (a,b) Fits to ∆Γ i msr .(c,d)The measured R µ/e combining the two signal channels in each q 2 interval.(e,f) Comparisons of the measured AFB and theoretical predications[33].Red circles, blue triangles, and black points with error bars are data; the error bars combine both statistical and systematic uncertainties.

d
Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany e Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education; Shanghai Key Laboratory for Particle Physics and Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240, People's Republic of China f Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of Modern Physics, Fudan University, Shanghai 200443, People's Republic of China g Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People's Republic of China h Also at School of Physics and Electronics, Hunan University, Changsha 410082, China e + e − collisions at E CM between 4.128 and 4.226 GeV, the D + and other background shapes are modeled by the individual simulated shapes taken from the inclusive MC sample.The signal and peaking background shapes are convolved with a Gaussian resolution function to account for the differences between data and simulation.The yields of the peaking backgrounds are fixed to the expectation from simulation, and the other yields are left free.The obtained signal efficiencies, signal yields, and resultant BFs are shown in Table1.The signal efficiencies have been corrected for small data-MC differences (overall factor f cor = 0.992 ∼ 1.018) in the η, π 0 reconstruction, µ + PID, requirements of E max γ, extra , N char extra , N π 0 extra , M η (′) µ + , M η (′) νµ , and χ 2 .

Table 2 .
Fitted parameters of hadronic FFs.Quantities in the first (second) parentheses are the least two significant digits of statistical (systematic) uncertainties.NDF is the number of degrees of freedom.

Table 4 .
The partial decay rates of D + s → ηµ + νµ in different q 2 intervals of data, where the uncertainties of partial decay rates are statistical only.

Table 5 .
The partial decay rates of D + s → η ′ µ + νµ in different q 2 intervals.Numbers in the parentheses are the statistical uncertainties.

Table 6 .
Statistical and systematic correlation matrices and relative uncertainties for the measured partial decay rates of D + s → ηµ + νµ in different q 2 intervals.